So I used the rgl package and created my own likelihood function to output the log likelihood of a sample from a normal distribution. I was doing this really just to learn how to program this myself so I could better understand how likelihood works and also how MLE works. Anyways, I noticed something particularly strange and I wanted to know if someone knew the answer here. When I plot the graph, it comes out in a folded curve shape, but I supposed I was expecting more of a cone type shape. Basically, what im curious about is why when the plot peaks at the sigma^2 value (on this axis, there is a good decline on both sides of the peak), the mu value stays roughly the same? It's as if once the sigma^2 parameter has reached the optimal level, the differences in likelihood between mu values are pretty small. For example, when I check the variance of the likelihoods of the maximum point of sigma (keeping it constant), it's 11.5. In contrast, when I check the variance of the mu's across that same point, the variance is 23402. Since I can't yet post images since I don't have enough reputation, I will just post my R-code that produces the graph.
#Define LL function
LL <- function(X, theta)
{
mu <- theta[1]
sigma2 <- theta[2]
log.likelihood <- 0
n <- length(X)
for (i in 1:length(X))
{
log.likelihood <- log.likelihood - (((X[i]-mu)^2)/(2*sigma2)) -
log(sqrt(2*pi*sigma2))
}
return(log.likelihood)
}
#Parameters
Mu <- 100
Sigma2 <- 50
#Sample
N <- 100
set.seed(1)
IQs <- rnorm(N, mean=Mu, sd=sqrt(Sigma2))
#Possible values to test
x <- posMu <- seq(80, 120, length.out=200)
y <- posSig <- seq(20, 60, length.out=200)
#x1 <- sort(x, decreasing=T)
#Produce LLs for plotting
LLlist <- NULL
for (m in 1:length(posMu)){
LLs <- NULL
for(s in 1:length(posSig)){
posTheta <- cbind(posMu[m],posSig[s])
LLs <- c(LLs, LL(IQs,posTheta))
}
LLlist <- cbind(LLlist,LLs, deparse.level=0)
}
z <- LLlist
#Find the approximate MLE
mLL <- which(LLlist == max(LLlist), arr.ind=TRUE)
cbind(posMu[mLL[2]],posSig[mLL[1]],LLlist[mLL])
#Graph the LLs
library(rgl)
open3d()
plot3d(mean(x),mean(y),mean(z), xlab="Mu", ylab="Sigma2", zlab="log L", xlim=c(min(x),max(x)), ylim=c(min(y),max(y)), zlim=c(min(z),max(z)))
surface3d(x, y, z, color=rainbow(length(x)))
So, is my code just wrong? Or is this what a LL curve should look like? If so, why is it that sigma^2 seems to show a clear curve and height whereas mu hardly differs at the maximum? Thanks in advance!