# home assignment 1, solution to simulation part
# intended solution used glm(family=quasi), but it does not work easily
# alternative solution, works but not elegant

library(stats4)

# part 1: n=10
x <- rep(c(5,10),each=5)
set.seed(250)
y <- rnorm(10,exp(0.1*x),sqrt(exp(0.1*x)))
ll <- function(beta=0.1)
         -sum(stats::dnorm(y, exp(beta*x), sqrt(exp(beta*x)), log=TRUE))
fit <- mle(ll)
ml10 <- coef(fit)
# fit <- glm(y ~ x-1, family=quasi(link="log", variance="mu"))
# does not work consistently, start values don't help
f <- function(beta) (sum(x*y)-sum(x*exp(beta*x)))^2
# note: estimating equation is: sum(x*y) = sum(x*exp(beta*x)), from (5.55) in MS
fit <- nlm(f,0)
mql10 <- fit$estimate
mql10min <- fit$minimum

for (i in 2:1000)
{ 
  y <- rnorm(10,exp(0.1*x),sqrt(exp(0.1*x)))
  ml10 <- c(ml10, coef(mle(ll)))
  fit <- nlm(f,0)
  mql10 <- c(mql10, fit$estimate)
  mql10min <- c(mql10min, fit$minimum)
}
names(ml10) <- NULL
mean(ml10)
sd(ml10)
mean(mql10[mql10min<0.01])
sd(mql10[mql10min<0.01])

# part 2: n=50
x <- rep(c(5,10),each=25)
set.seed(260)
y <- rnorm(50,exp(0.1*x),sqrt(exp(0.1*x)))
fit <- mle(ll)
ml50 <- coef(fit)
fit <- nlm(f,0)
mql50 <- fit$estimate
mql50min <- fit$minimum

for (i in 2:1000)
{ 
  y <- rnorm(50,exp(0.1*x),sqrt(exp(0.1*x)))
  ml50 <- c(ml50, coef(mle(ll)))
  fit <- nlm(f,0)
  mql50 <- c(mql50, fit$estimate)
  mql50min <- c(mql50min, fit$minimum)
}
names(ml50) <- NULL
mean(ml50)
sd(ml50)
mean(mql50[mql50min<0.01])
sd(mql50[mql50min<0.01])