cauchy.mle2 <- function(x, maxiter=20, eps=1e-8, 
                        itertrace=FALSE, maxattempts=100){

# initialize parameter vector
p <- median(x)
n <- length(x)

# calculate log-likelihood at initial parameter estimates
ll <- - n*log(pi) - sum(log(1 + (x-p)^2))

iter <- 0

while (iter < maxiter){
iter <- iter+1

# calculate shortcuts
denom <- (1 + (x-p)^2)
denom2 <- denom^2
xminusp <- x-p

# calculate second derivatives at current parameter estimates
D <- 2* sum((- denom + 2*xminusp^2)/denom2)

# calculate first derivaties at current parameter estimates
f <- 2*sum(xminusp/denom)     # first derivative of loglik wrt theta

pnew <- p - f/D

# calculate log-likelihood at new parameter estimates
llnew <- - n*log(pi) - sum(log(1 + (x-pnew)^2))

# if likelihood decreased, try random step sizes in both directions
if (llnew<ll){
attempts <- 0
while(llnew<ll && attempts<maxattempts){
attempts <- attempts+1
if (itertrace){
cat("Log-Likelihood decreased from ", ll, "to ", llnew, "\n",
"Attempting random stepsize. \n")
}
pnew <- p + (runif(1)-.5)*f/D

# calculate log-likelihood at new parameter estimates
llnew <- - n*log(pi) - sum(log(1 + (x-pnew)^2))
}
if (llnew<ll) {
cat("\n *** ERROR: Could not improve likelihood.*** \n Current estimates: \n\n")
return(list(parameter_estim = p,loglik_1stder = f, 
fisherinfo=D, loglikelihood = ll))
break
}
}

ll <- llnew

# if max(change in estimate) < eps or max(function value)<1e-10 then converged

if(max(abs(p-pnew))<eps || max(abs(f))<1e-10 ) {
cat("*** CONVERGED in ", iter, "ITERATIONS. ***\n", 
"MLE ESTIMATOR: theta = ", pnew, "\n")
return(list(parameter_estim = pnew,loglik_1stder = f, 
fisherinfo=D, loglikelihood = ll))
break
}
else {
p <- pnew
if(itertrace){
cat("*** Iter.No:", iter, " Current root = ", p,
 "\n \t\t 1st derivative of log-likelihood:", f, "\n")
}
}
}
return(cat("*** DID NOT CONVERGE IN", maxiter, "ITERATIONS. ***\n"))
}