I was wondering about how to integrate via Monte Carlo Methods. Suppose that $Y_{1:n}$ is a sample of $n$ values, there are simulated from a mixture of normal distributions, $Y_{1:n} \sim p_{0} \mathcal{N}(\mu_{0},1)+ (1-p_{0}) \mathcal{N}(\lambda_{0}, 1)$ where $p_{0},\mu_{0}, \lambda_{0}$ are known. Assume that I put a prior on the parameters and assume that $Y_{1:n}$ are simulated from a $p \mathcal{N}(\mu,1)+ (1-p) \mathcal{N}(\lambda, 1)$ model, where, $\mu,\lambda \sim \mathcal{N}(0,1)$, and $p \sim U(0,1)$. I am trying to find a Monte Carlo way of computing $\int\limits_{[0,1] \times \mathbb{R}^{2}} \prod\limits_{i=1}^{n} \big[p \phi(y_{i}-\mu)+(1-p)\phi(y_{i}-\lambda)\big] \phi(\mu)\phi(\lambda) 1 d \mu d\lambda dp$.


Simulate from the prior, calculate the likelihood for a fixed dataset (your expression is a function of chosen $y_1,\ldots,y_n$), and then store the resulting number. Do that over and over again and calculate the average of your numbers.

Pick a number $B$ for how many times you want to simulate. For $k=1,\ldots,B$,

  1. draw $\mu^k, \lambda^k, p^k$
  2. evaluate and store $g(\mu^k, \lambda^k, p^k) = \prod\limits_{i=1}^{n} \big[p^k \phi(y_{i}-\mu^k)+(1-p^k)\phi(y_{i}-\lambda^k)\big]$
  3. report $B^{-1}\sum_{k=1}^Bg(\mu^k, \lambda^k, p^k)$.

This works for any choice of $y_1,\ldots,y_n$.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.