# Double Monte Carlo Integration of mixed normal distribution

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$.