# 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