# Solution to STA 410/2102 Assignment #2, Spring 2003 - Commands for tests done.


# Small data set for testing.

testx <- c(1.56,5.20,1.49,7.52,3.43,5.96,5.16,7.62,3.93,4.22)


# PART I. 

cat("\nTESTS FOR PART I\n\n")

# Try estimating mu with omega fixed at 1, from various starting points.

print (logistic.mle1(testx,10,debug=T))  # Default initial value of median

for (mu0 in c(0.1,1,3,10))
{ print (logistic.mle1(testx,10,mu0,debug=T))
}


# PART II. 

cat("\nTESTS FOR PART II\n\n")

# Try estimating mu and omega, from various starting points.

print (logistic.mle2(testx,10,debug=T))  # Default initial values

for (mu0 in c(0.1,1,3,10))
{ for (omega0 in c(0.1,0.3,1,3,10))
  { print (logistic.mle2(testx,10,mu0,omega0,debug=T))
  }
}


# PART III.

cat("\nTESTS FOR PART III\n\n")

# Try estimating mu with omega fixed at 1, from various starting points,
# using nlm.

cat ("\nInitial value:",median(testx),"\n\n")
print (nlm (function (mu) -logistic.log.lik(testx,mu), median(testx)))

for (mu0 in c(0.1,1,3,10))
{ cat ("\nInitial value:",mu0,"\n\n")
  print (nlm (function (mu) -logistic.log.lik(testx,mu), mu0))
}

# Try estimating mu and omega, from various starting points, using nlm.

cat ("\nInitial values:",median(testx),
      (quantile(testx,0.75)-quantile(testx,0.25))/log(9),"\n\n")
print (nlm (function (v) -logistic.log.lik(testx,v[1],v[2]), 
       c (median(testx), (quantile(testx,0.75)-quantile(testx,0.25))/log(9))))

for (mu0 in c(0.1,1,3,10))
{ for (omega0 in c(0.1,0.3,1,3,10))
  { cat ("\nInitial values:",mu0,omega0,"\n\n")
    print (nlm (function (v) -logistic.log.lik(testx,v[1],v[2]), c(mu0,omega0)))
  }
}

# Compare times on a large data set.

set.seed(1)   # Set random number seed so results will be reproducible.

u <- runif(10000)
x <- 4 + 3*log(u/(1-u))

cat("\nTime for Newton-Raphson\n\n")
print(system.time(print (logistic.mle2(x,2))))

cat("\nTime for nlm\n\n")
print(system.time(print (nlm (function (v) -logistic.log.lik(x,v[1],v[2]), 
                   c (median(x), (quantile(x,0.75)-quantile(x,0.25))/log(9))))))