# home assignment 1, program for simulation part (see 2004 version for details)
# 10/2009: revised version with more elaborate specification of minimization

library(stats4)

# part 1: n=10, beta=0.1
x <- rep(c(5,10),each=5)
set.seed(250)
ll <- function(beta=0.1)
         -sum(stats::dnorm(y, exp(beta*x), sqrt(exp(beta*x)), log=TRUE))
ml10<- numeric()
f <- function(beta) log(1+(sum(x*y)-sum(x*exp(beta*x)))^2)
# fit <- glm(y ~ x-1, family=quasi(link="log", variance="mu"))
# does not work consistently, start values don't help
# old 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
mql10<- numeric()
mql10min<- numeric()

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

# part 2: n=50, beta=0.1
x <- rep(c(5,10),each=25)
set.seed(260)
ml50<- numeric()
mql50<- numeric()
mql50min<- numeric()

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


# part 3: n=10, beta=1
x <- rep(c(5,10),each=5)
set.seed(250)
ml10<- numeric()
mql10<- numeric()
mql10min<- numeric()

for (i in 1:1000)
{ 
  y <- rnorm(10,exp(1*x),sqrt(exp(1*x)))
  ml10 <- c(ml10, coef(mle(ll))[[1]])
  fit <- nlm(f,.95,typsize=c(1),fscale=1)
  mql10 <- c(mql10, fit$estimate)
  mql10min <- c(mql10min, fit$minimum)
}
mean(ml10)
sd(ml10)
mean(mql10[mql10min<4])
sd(mql10[mql10min<4])

# part 4: n=50, beta=1
x <- rep(c(5,10),each=25)
set.seed(260)
ml50<- numeric()
mql50<- numeric()
mql50min<- numeric()

for (i in 1:1000)
{ 
  y <- rnorm(50,exp(1*x),sqrt(exp(1*x)))
  ml50 <- c(ml50, coef(mle(ll))[[1]])
  fit <- nlm(f,.95,typsize=c(1),fscale=1)
  mql50 <- c(mql50, fit$estimate)
  mql50min <- c(mql50min, fit$minimum)
}
mean(ml50)
sd(ml50)
mean(mql50[mql50min<8])
sd(mql50[mql50min<8])


# for beta=10, all mle values are equal to 10, and the mql estimation
# is essentially impossible by the procedure used above because the 
# target function is incredibly peaked around the true value (10). Basically,
# the observations are out of any sensible scale.