# dowloaded from https://rpubs.com/Koba/MLE-Normal
# 1/27/2022
## the name of the tab-delimited text file containing the data
file_name <- "Ant_Nest_Data.txt"

## PROGRAM CODE ##

## load the ant nest density data from the text file
nd <- read.table(file_name, header=T)
forest <- nd$NestsPerQuadrat[nd$Habitat=="Forest"]


## Example with mean 10 and variance 25
x <- rnorm(1000,10,sqrt(25))
mo <- mean(x)
vo <- var(x)
m <- seq(9,11,.01)
v <- seq(20,32,length.out=length(m))
#######################################

## Example with class data 
# x <- forest
# mo <- mean(x)
# vo <- var(x)
# m <- seq(2,10,0.1)
# v <- seq(1,40,length.out=length(m))
#######################################

normalF <- function(parvec) {
  # Log of likelihood of a normal distribution
  # parvec[1] - mean
  # parvec[2] - standard deviation
  # x - set of observations. Should be initialized before MLE
  sum ( -0.5* log(parvec[2]) - 0.5*(x - parvec[1])^2/parvec[2] )
}

normalF(c(mean(x),sqrt(var(x)))) # log likelihood function value for given x and mu=sd=1 


lik <- array(0,c(length(m),length(v)))
likp <- array(0,c(length(m)*length(v),3))
lik <- as.matrix(lik)
s <- 0
for (i in 1:length(m)){
  for (j in 1:length(v)) {
    s <- s+1
    lik[i,j] <-normalF(c(m[i],v[j])) # log likelihood function value  
    likp[s,2] <- m[i]
    likp[s,1] <- v[j]
    likp[s,3] <-normalF(c(m[i],v[j])) # log likelihood function value  
  }
}


par(mfrow=c(1,1))
#image(lik)
likp[,3] <- likp[,3]-min(likp[,3])
likp <- as.data.frame(likp)
print(paste(c("Obs",round(vo,1),"MxLk", round(likp[,1][which(likp[,3]==max(likp[,3]))],1))))
print(paste(c("Obs",round(mo,1),"MxLk", round(likp[,2][which(likp[,3]==max(likp[,3]))],1))))

########################
library(plotly)      
fig <- plot_ly(x=unique(likp[,1]),y=unique(likp[,2]),z=lik) %>% add_surface()

fig