# STA 2102/450 Spring 2000 - S Functions for Assignment #2, Part 1 - R. M. Neal


# LOG LIKELIHOOD FUNCTION FOR THE CAUCHY MODEL.  The arguments are the location
# parameter and the vector of data points.  The constant factor of pi in the
# likelihood is omitted.

cauchy.log.likelihood <- function (theta, x)
{
  sum (-log(1+(x-theta)^2))
}


# DERIVATIVE OF LOG LIKELIHOOD FUNCTION.  Arguments as above.

cauchy.ll.deriv <- function (theta, x)
{
  sum (2*(x-theta)/(1+(x-theta)^2))
}


# SECOND DERIVATIVE OF LOG LIKELIHOOD FUNCTION.  Arguments as above.

cauchy.ll.2nd.deriv <- function (theta, x)
{
  d <- 1+(x-theta)^2
  sum (4*(x-theta)^2/d^2 - 2/d)
}


# FIND MLE FOR CAUCHY MODEL USING NEWTON-RAPHSON ITERATION.  The arguments
# are the vector of data point, an initial guess for the location parameter
# (default is the median of the data points), and the number of iterations
# of Newton-Raphson iteration to do (default is 9). 
#
# The values of theta and the log likelihood after each iteration are printed 
# so that we can see what's happening.

cauchy.mle <- function (x, theta0=median(x), n=9)
{
  theta <- theta0

  cat (0, theta, cauchy.log.likelihood(theta,x), "\n")

  for (i in 1:n)
  { theta <- theta - cauchy.ll.deriv(theta,x) / cauchy.ll.2nd.deriv(theta,x)
    cat (i, theta, cauchy.log.likelihood(theta,x), "\n")
  }
  
  theta
}


# TEST HOW WELL IT WORKS.

# A test data set.

x<-c(1.8,-3.1,2.3,1.1,11.9)

# Plot the log likelihood function for this data.

th<-seq(-15,15,by=0.1)
ll<-rep(0,length(th))
for (i in 1:length(th)) ll[i] <- cauchy.log.likelihood(th[i],x)
plot(th,ll,type="l")

# Try to find the MLE from various starting points.

cat("From median:\n\n")
mle1 <- cauchy.mle(x)

cat("\nFrom -5:\n\n")
mle2 <- cauchy.mle(x,-5)

cat("\nFrom 5:\n\n")
mle3 <- cauchy.mle(x,5)

cat("\nFrom 9:\n\n")
mle4 <- cauchy.mle(x,9)

cat("\nFrom 12:\n\n")
mle5 <- cauchy.mle(x,12)

# Show the results found above on the plot (if they fit).

abline(v=c(mle1,mle2,mle3,mle4,mle5))