# Daniel Bernoulli's MLE

# 1D. First, his error distribution with r=1 (r FIXED)

# 1/2 circle error distribution,
# centered on theta, with radius r=1

dx=0.01

Y = seq(-1,2,dx)

pdf = function(x,theta,r) {
 if( abs(x-theta) <= r) density=( 2/(pi*(r^2)) ) *
                          sqrt(r^2-(x-theta)^2 )
 if( abs(x-theta) > r)  density=0
 return(density)
}  

PDF = unlist( lapply(y, pdf, theta=0, r=2 ))
sum(PDF*dx)

observed.y = c(0,0.2,1) # **********
n = length(observed.y)

R=1

par(mfrow=c(3,3),mar = c(4,4,3,1))

THETA = seq(0.05,0.85,0.1)
for (parameter in THETA ) {
 MAIN=""
 if(parameter==THETA[2]) MAIN="Bernoulli (<=1778), R(2016)"
 PDF = unlist( lapply(Y, pdf, theta= parameter, r=R ))
 plot(Y, PDF,type="l",
   xlim=c( min(THETA)-R, max(THETA)+R), ylim=c(0,.8),
   main=MAIN )
 points(observed.y,rep(0,n),
  pch=19)
 points(parameter,0,pch=18, col="red")
 text(parameter,0,toString(parameter), 
  col="red", adj=c(0.5,-0.5))
 if(parameter==THETA[1]) text(1.6,0.5,
 "log-likelihood
 contributions", adj=c(1,0))
 if(parameter==THETA[9]) {
 	text(-0.9,0.1,"observed
data", adj=c(0,0.5))
 	points(-1.0,0.1,pch=19)
 }  
 if(parameter==mean(THETA)) text(parameter,0.1,"theta", 
  col="red", adj=c(0.5,-0.5))
 PDF = unlist(lapply(observed.y,pdf,theta=parameter,r=R))
 LOGLik = round(log(PDF),2)
 overall = logLik = round(sum(log(PDF)),2)
 ADJ = 1.1*(observed.y < parameter)-
       0.1*(observed.y >= parameter)
 for(i in 1:n) {
 	segments(observed.y[i],0,
                        observed.y[i],PDF[i] )
    text(observed.y[i],PDF[i], 
         toString(LOGLik[i]), adj=c(ADJ[i],0) )
 }
 text(parameter,0.8,
    paste("LOGLIK = ", 
         toString(overall)), adj=c(0.5,1), cex=1.5 ) 
}


###################################################

# 2D. His error distribution with r ESTIMATED 

####################   estimating both r and theta

observed.y = c(0,0.2,1)


#  good to work with the log(r) scale to avoid negative r

# dy assumed to be small, so that probability of
# observing each reported y value is close to pdf(y)*dy

# the log of the probability = log(pdf) + log(dy)

# dy common to all observations, and does not involve
# the parameters, so we focus on log density itself 

logLik2D = function(par) {
 theta=par[1] 
 r = exp(par[2]) 
 densities =( 2/(pi*(r^2)) ) *
          sqrt(r^2-(observed.y-theta)^2 )
  return( sum( log(densities) ) )
} 

logLik2D(   c(1,log(1)) )

ml.fit = optim(c(0.4,1), logLik2D,
               hessian=T,method="BFGS",
               control=list(fnscale=-1,
                   trace=5,
                   parscale=c(0.1,0.1))); 

ml.fit

exp(ml.fit$par[2]) 

# jh 2016.10.17