# 610 fall 2014; compare with HW1

# mleCauchy.R
# MLE for Cauchy location:
# Draw log-likelihood for a sample from standard Cauchy			
#	with expected value superimposed.  Maybe show approximating
#	quadratic for log-likelihhod.
 
setwd("/Users/pollard/Google Drive/G_DP-courses/610/610.fall2013/610Handouts/")


# the grid points at which to evaluate G_n and G_0:
GridCauchy <- seq(-10,10,by=0.001)			

# the function G_0(t):	
makeG0 <-function(t){ 	
	integrand <- function(x) {log(dcauchy(x,loc=t))*dcauchy(x)}
	integrate(integrand,-Inf,Inf)
}
# add code to let the function accept vector arguments:
makeG0 <- Vectorize(makeG0,vec="t")

if(!exists("G0")){
	# extract the bit we need from the output:
	G0 <- unlist(makeG0(GridCauchy)[1,])	
}

# derivative of log density at theta=0;
# 	used to calculate Z_n:
ell <- function(z){
	2*z/(1+z^2)
}

# text for subtitles:
G0text <- "G_0 is dotted blue line"
Qtext <- " the quadratic approximation is in purple"
ttext <- "\nvertical lines locate each argmax"

drawCauchyMLE <- function(
	n = 100,		# sample size
	Grid = GridCauchy,
	pop=G0,
	xrange = c( min(Grid),max(Grid) ),
	showquad=F,
	showmle=T,
	showapprox=T
	)
{
	# reduce the domain for the plots
	pop <- pop[ Grid >= xrange[1] & Grid <= xrange[2] ]
	Grid <- Grid[ Grid >= xrange[1] & Grid <= xrange[2] ]
	
	# calculate G_n:
	xx <- rcauchy(n)	# the observations
	Gn <- apply(log(dcauchy( outer(xx,Grid,"-"))),2,mean)
	mle <- Grid[Gn == max(Gn)]
	
	# vertical range for the picture:
	ymax <- max(pop,Gn)+0.1
	ymin <- min(pop,Gn)- 0.1
	
	# first draw G_n:
	plot(Grid,Gn,type="l",xlab="",ylab="G_n(t)",ylim=c(ymin,ymax),xlim=xrange)
	title(paste("log-likelihood for a sample of size ",n," from the Cauchy"))
	if(showmle){
		abline(v=mle)
	}
	# then add G_0:
	lines(Grid,pop,type="l",col="blue",lty=3)
	
	
	if(!showquad){	# no quadratic to be shown:
			title(sub=G0text)
	}else{	# add approximating quadratic:
		lin <- mean(ell(xx))	# Z_n/sqrt(n)
		quad <- Gn[Grid == 0] + Grid*lin - Grid^2/4
		lines(Grid,quad,col="purple")
		subtitle <-paste(G0text,Qtext,sep=",")
		if(showapprox){
			abline(v=2*lin,col="purple")
			subtitle <-paste(subtitle,ttext,sep=",")
		}
		title(sub= subtitle)
	}
}