# Monte Carlo integration and mixture distribution

How could I estimate an integral using Monte Carlo method when I have a mixture distribution? For example I want to estimate the below integral:

$$I=\int_{0}^{1}f(x)dx$$

And my distribution is:

$$p(x)=\sum_{k=1}^{n}\lambda_kp_k(x)$$

And $$\langle I^N\rangle= \frac{1}{N}\sum_{i=1}^{N}\frac{f(x_i)}{p(x_i)}$$ is a Monte Carlo estimate of $$I$$.

Is $$p(x_i)$$ the mixture density or conditional density?

Both solutions work in the sense that, if $$Z$$ denotes the index of the component of the mixture used to produce $$X$$, there is equality between the expected weights: $$\mathbb E^X\left[\dfrac{f(X)}{p(X)}\right] = \mathbb E^{X|Z}\left[\dfrac{f(X)}{p(X|Z)}\right] = I$$

Here is a quick illustration of the fact:

#target function
f <- function(x) sin(3*x)
#mixture
p <- function(x) .3*dbeta(x,1,2)+.7*dbeta(x,2,4)
#complete sample
N=1e6
z=sample(0:1,N,replace = TRUE,prob = c(3,7))
y=rbeta(N,shape1 = 1+z,shape2 = 2+2*z)
w1=f(y)
w2=w1/p(y)
w1=w1/dbeta(y,shape1 = 1+z,shape2 = 2+2*z)
#estimators
est1=cumsum(w1)/(1:N)
est2=cumsum(w2)/(1:N)


with both estimators having similar behaviour in that simple example (blue for the unconditional version and orange for the conditional version): #plots
plot(log(1:N),est1,ty="l",ylim=c(.3,.8),col="midnightblue",xlab="log(n)",ylab="")
lines(log(1:N),est2,col="orange2",lty=2)
abline(h=integrate(f,0,1)\$value,col="tomato",lty=3)


I think you should use following steps

lets we want to generate $$m$$ samples to use in Monte-Carlo estimator

1) Generate $$m$$ sample from following distribution (with replacement )

$$\begin{eqnarray} \begin{array}{c|ccc} k & 1 & 2 & \cdots & n \\ \hline probability & \frac{\lambda_1}{\sum \lambda_i} & \frac{\lambda_2}{\sum \lambda_i} & \cdots & \frac{\lambda_n}{\sum \lambda_i} \\ \end{array} \end{eqnarray}$$

and save it(in vector named $$K$$,In standard notation $$\sum \lambda_i=1$$ So you can just ignore $$\sum \lambda_i$$).

2) Now for each $$K=k$$ generate sample from $$\mathbb P_k(x)$$

lets do it with a simple example

$$X\sim .2 Normal(0,1)+.8 Normal(5,1)$$

lets we want to estimate $$E(X)$$ (that exact value is $$4$$ and estimate should be near $$4$$). I do it in R software.

  > m<-100000
> k.vector<-sample(1:2,size=m,p=c(.2,.8)/sum(c(.2,.8)),replace=T)
> l1<-length(which(k.vector==1))
> l2<-m-l1
> random.g<-c(rnorm(l1,0,1),rnorm(l2,5,1))
> mean(random.g)
 4.003443


This algorithm is a simple method based on conditional porbability

$$f_X(x)=\sum \mathbb P(X=x|K=k) P(K=k)$$ so in first step we produce one $$k$$ and then generate one $$x$$ form $$\mathbb P(X|K=k)=\mathbb P_k(x)$$. Next use this generated sample to estimate what you want.

• What I asked was vague so I edited my post. Apr 9, 2020 at 10:54