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

Question

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!

Answer 1

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]]

Answer 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)))}