# How do you perform a two-dimensional grid-search for the MLE in R?

I'm currently working on parameter estimation for a mixture of two normal distributions. I have two parameters I am trying to estimate using the maximum likelihood method: the probability and the population mean. In the equation above, $$f_1$$ follows a normal distribution with mean 0 and standard deviation 1, while $$f_2$$ follows a normal distribution with mean $$\mu$$ and standard deviation 1.

I have already defined the log-likelihood function and also generated a vector $$x$$ of 100 observations setting $$\pi$$ equal to 0.25, and $$\mu$$ equal to 1. I would like to find a combination of these two parameters that generates the highest log-likelihood, but I am not sure how to perform an exhaustive search in two dimensions. (I will later implement an optimizing function, but I would like to perform a grid search first).

Here is the code I have thus far:

loglike = function(pi, mu, x) {

llh = NULL # initializing an empty list

for (i in x) {
f1 = dnorm(i, 0, 1)
f2 = dnorm(i, mu, 1)
llh = c(llh, log(pi*f1 + (1-pi)*f2)) # appending the values to the list
}
return(sum(llh))
}

set.seed(4185)

group = sample(c(1,2), 100, replace=TRUE, prob=c(.25,.75))
rn1 = rnorm(100,0,1)
rn2 = rnorm(100,1,1)
x = rn1*(group==1) + rn2*(group==2)

piSeq = seq(0, 1, .01)
muSeq = seq(min(x), max(x), .01)


Any help or advice is appreciated. Thanks!

• Notice the slight modification in the code (I had the rows and columns switched for ll_save – dlnB May 1 '20 at 19:59
• @dlnB Thanks! It works. – Aspire May 1 '20 at 20:19

Here is an algorithm for an exhaustive search (i.e. searching the full grid rather than searching until you have a local max. You need to choose a step size for the $$\pi$$ loop and for the $$\mu$$ loop (it looks like you chose .01 in your code, so I will use that here).

mustep<-.01
pistep<-.01

mu_seq<-seq(min(x), max(x),mustep)
pi_seq<-seq(0,1,pistep)

ll_save<-matrix(, nrow = length(pi_seq), ncol = length(mu_seq))

for (i in pi_seq){
for (j in mu_seq){

loglike(i,j,x)

}
}

indices<-which(ll_save == max(ll_save), arr.ind=TRUE)

pi<-pi_seq[indices[1]]
mu<- mu_seq[indices[2]]


You can create the mixture without the need of discarding most of the random values generated in the Gaussians. You can append two random Gaussians, in your experiment you don't need them to be mixed, so you can append them, each of them in the proportion you want.

You can also obtain the likelihood or loglikelihood without the loop. You can apply the dnorm function to vectors.

You can try with this, where I try 1000 values of the mean from 0 to 10. My code by the way is not the best.

mus=c()
loglikelihood=c()
for (i in 1:1000){
mus[i]=i/100
loglikelihood[i]=sum(log(0.25*dnorm(x,0,1)+0.75*dnorm(x,i/100,1)))}

• Sorry I focused on using dnorm with vectors. You can easily apply the idea with different values of the proportions, for example repeting the code for different N values of the proportion, you can create an 1000*N matrix of loglikelihoods. Of course, you can also change 1000 and use the number you prefer. – javierazcoiti Apr 30 '20 at 16:06