2015年4月10日星期五

Use Newton-Raphson Method to calculate MLE of a Cauchy distribtuion


mlecauchy=function(x,toler=.001){      #x is a vector here
startvalue=median(x)
n=length(x);
thetahatcurr=startvalue;
# Compute first deriviative of log likelihood
firstderivll=2*sum((x-thetahatcurr)/(1+(x-thetahatcurr)^2))
# Continue Newton’s method until the first derivative
# of the likelihood is within toler of 0.001
while(abs(firstderivll)>toler){
# Compute second derivative of log likelihood
secondderivll=2*sum(((x-thetahatcurr)^2-1)/(1+(x-thetahatcurr)^2)^2);
# Newton’s method update of estimate of theta
thetahatnew=thetahatcurr-firstderivll/secondderivll;
thetahatcurr=thetahatnew;
# Compute first derivative of log likelihood
firstderivll=2*sum((x-thetahatcurr)/(1+(x-thetahatcurr)^2))
}
list(thetahat=thetahatcurr);
}
x<-c(-1.94,0.59,-5.98,-0.08,-0.77)
mlecauchy(x,0.0001)

###########result########'
####   $thetahat
####  [1] -0.5343968

#####use the R biuld in funcion ######

optimize(function(theta) -sum(dcauchy(x, location=theta, log=TRUE)), c(-100,100)) #we use negative sign here,
                                                                                  #by default it return a minimum
#optimize(function(theta) sum(dcauchy(x, location=theta, log=TRUE)),maximum=T, c(-100,100))

#########Result###############

###  $minimum
###  [1] -0.5343902
###   $objective
###  [1] 11.2959


Almost the same for the two methods.