###*********************###
###*** MLE_n_LogLH.r ***###
###*********************###
rm(list=ls())

## variables #################################################################
# (also, need to modify 1:mtext line, 2: LogLH function, for different distribution)

# parameter
iParam="ƒÎ"
# minimum and maximum interval values
iMin=0.4
iMax=0.8
# n and x
iN = 1000
iX = 620

### Functions #################################################################

# Log Likelihood Function of theta
FuncLogLH = function(ValuesPassed){
    log((ValuesPassed^iX)*(1-ValuesPassed)^(iN-iX))
}

# Since I couldn't find any command to omit 'Inf' in data, so first I change to 'NA', 
# then later use "na.omit". (this is for the plot below)
ChangeInfToNA = function(iData){
   iDataBKUP = iNoInfData = numeric(length(iData))
   iDataBKUP = iData
   for (i in 1:(length(iData))){
      {
      if (iDataBKUP[i]==Inf || iDataBKUP[i]==-Inf ){
        iNoInfData[i] = NA
      }#if
      else iNoInfData[i] = iDataBKUP[i]
      }#if-else
   }#for
   return(iNoInfData)
}#function

## MAXIMUM LIKELIHOOD ESTIMATE -"optimize"- ###################################
iThetaHatOptimized = optimize(FuncLogLH, c(iMin,iMax), maximum=TRUE)

## MAXIMUM LIKELIHOOD ESTIMATE ################################################
iThetaHatByMLE = iThetaHatOptimized$maximum

## PLOT TO SEE HOW THE LOG LOKELIHOOD OF THETA LOOKS LIKE #####################
par(omi=c(0,0,0.3,0))
# x-axis
iTheta = seq(iMin,iMax,0.01)
# y-axis
iLogLH = FuncLogLH(iTheta)
# find min and max y-axis values, except "(-)Inf"
iYMin = min(na.omit(ChangeInfToNA((FuncLogLH(iTheta)))))
iYMax = max(na.omit(ChangeInfToNA((FuncLogLH(iTheta)))))
# ...and plot.
plot(iTheta, iLogLH, type="l", xlab=bquote(.(iParam)), 
     ylab=bquote(paste("logLH(",.(iParam),")")), 
     main=bquote(paste("Log Likelihood(",.(iParam),"), and Maximum Likelihood Estimate; ",hat(.(iParam)))))
arrows(iThetaHatByMLE,iYMax,iThetaHatByMLE,iYMin,col="red")
text(iThetaHatByMLE+0.05,iYMin+((iYMax-iYMin)/5),label=bquote(paste(hat(.(iParam))," = ",.(round(iThetaHatByMLE,4)))),
     pos=4,col="red")
mtext(side=3, line=-1.0, outer=T, bquote(paste("X~Bin(n, p): n=",.(iN),", X=",.(iX))), cex=1.3)
par(omi=c(0,0,0,0))

### From the handout ###
iPie = seq(0.001,0.999,0.001)
iLike = dbinom(620,1000,iPie)
X11()
plot(iLike, type="l", col="blue", xlab="n",ylab="Density", main="Numerical Maximization")
text(0,0.015,pos=4,label=bquote(paste("pie[like==max(like)] is ",.(iPie[iLike==max(iLike)]))))
########################