I have a likelihood function that contains an integral of a latent parameter. I would like to numerically integrate it using Monte Carlo, as in, $L = \prod_{i=1}^N \int f(X, \tilde{\theta}; \beta) d F_\theta (\tilde{\theta})$, with other parameters of interest $\beta$. I assume that this latent parameter follows a mixture of two Normal distributions with mean zero and unknown standard deviations, $\theta \sim p N(0,\sigma_1) + (1-p) N(0,\sigma_2)$.
I know how to employ this when I have a non-mixture distribution. If $\theta \sim N(0,\sigma)$ where $\sigma$ is unknown, I would sample $mc \sim N(0,1)$ $M$ times, use it as an input of the likelihood function, and transform it within it by multiplying $mc_i$ with the parameter $\hat{\sigma}$. The likelihood would be simply $L= \prod_{i=1}^N M^{-1} \sum_{k=1}^M f(X_i, mc_k * \hat{\sigma}; \hat{\beta})$, and I minimize it with respect to $(\sigma, \beta)$.
However, in the case of a mixture, I'm not sure how to implement the monte carlo integration. I know that sampling from a mixture of two normals involve the repetition of the following steps:
1) Generate a r.v. $u \sim U(0,1)$
2) If $u < p$, then generate a sample from $N(0,\sigma_1)$, if $u \geq p$, then from $N(0,\sigma_2)$, where $p$ is the mixture probability
So, mirroring the non-mixture case, what I did was to sample $mc \sim N(0,1)$ and $u \sim U(0,1)$, and use them as inputs of the likelihood function. Within the likelihood, I need to transform $mc_i$ from $N(0,1)$ to the mixture distribution. My question is in how to do this effectively.
If I transform it using the algorithm detailed above, whenever I use gradient-based methods, the parameter for the mixture probability, $p$, never changes from its initial guess. On the other hand, using Nelder-Mead lets $p$ change and arrive at a solution.
Note that I'm not trying to write down the likelihood for the mixture. I do not have data that comes from a mixture of normals, I want to integrate out a latent variable that I assume comes from a mixture.
I was wondering if I'm sampling it correctly. I know it is the correct way if I simply want to generate some fake data, but there must be another way of writing it down within the likelihood.
function likelihood(param, MC, data)
baseprob = MC[:,2] # Uniform(0,1)
draws = MC[:,1] # Normal(0,1)
mix = 1 / (1 + exp(-param[8])) # bound guess to lie in (0,1)
selection = baseprob .> mix # boolean array indicating which dist to sample
draws = @. selection * draws * param[6] + (1 - selection) * draws * param[7]
(... write down the likelihood and integrate it using draws...)
end
I've read about the E-M algorithms for estimating the parameters of mixture models, but I can't see how they are applicable here. Could someone point me in the right direction?
