Homework 7  Due March 28


Suppose X1, X2, ... Xn  are independent random variables.

X_i = Bernoulli( p_i = 0.4 ) for odd i
X_i = Bernoulli( p_i = 0.6 ) for even i

Show that

         exp{ 1/(n sqrt n) [\sum_{i=1}^n ( X_i - p_i )]^2 }  

converge to a limit in probability as n go to infinity. (what is the limit?)

==============================================================================


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

# initialize parameter vector
f <- p <- c(NA,NA)		
p[1] <- median(x)
p[2] <- IQR(x)/2

n <- length(x)

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

# initialize matrix of second derivatives
D = matrix(ncol=2, nrow=2)

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

	# calculate shortcut for denominator
	denom = (p[2]^2 + (x-p[1])^2)
	denom2 = denom^2

	# calculate shortcut for x-p[1]
	x_p1 = x-p[1]

	
	# calculate second derivatives at current parameter estimates
	D[1,1] = 2* sum((-p[2]^2+x_p1^2)/denom2)
	D[1,2] = -4*p[2]* sum(x_p1/denom2)
	D[2,1] = D[1,2]
	D[2,2] = -n/p[2]^2 - D[1,1]

	# calculate first derivaties at current parameter estimates
	f[1] = 2*sum(x_p1/denom)     # first derivative of loglik wrt theta
	f[2] = n/p[2] - 2*sum(p[2]/denom) # 1st deriv of loglik wrt gamma

	pnew = p - solve(D,f)
	
	# if estimate for gamma is negative, replace with small value
	if (pnew[2]<=0) pnew[2] = 1e-5
		

	# calculate log-likelihood at new parameter estimates
	llnew = n*log(pnew[2]) - n*log(pi) - sum(log(pnew[2]^2 + (x-pnew[1])^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)*2*solve(D,f)
		
	# if estimate for gamma is negative, replace with small value
			if (pnew[2]<=0) pnew[2] = 1e-5
		
	# calculate log-likelihood at new parameter estimates
			llnew = n*log(pnew[2]) - n*log(pi) - 
				sum(log(pnew[2]^2 + (x-pnew[1])^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))
		}
	}
	
	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 = ", p[1], " gamma = ", p[2], "\n")
		return(list(parameter_estim = p,loglik_1stder = f, 
					fisherinfo=D, loglikelihood = ll))
		}
	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"))
}