Monte Carlo Methods

A Monte Carlo method is a computational / numerical method that uses random numbers to compute / estimate a quantity of interest. The quantity of interests may be the mean of a random variable , functions of several means , distributions of random variables or high dimensional integrals.

Basically, Monte Carlo methods may be grouped into two types:

The direct/simple/classical Monte Carlo methods involves generating identical and independent sequences of random samples. And the other, in a sense, involves generating a sequence of random samples, which are not independent, and is the Markov Chain Monte Carlo methods.

Monte Carlo Integration

Let f(x) be a function of x and suppose we are interested in computation of the integral:

I= \int_0^1 f(x) dx

We can write the integral as,

I=\int_0^1 f(x) p(x) dx =E(f(X))\\
\textit{; where p(x) is the pdf of a r.v. X $\sim$ Unif(0,1)} 

Now suppose that x_1,x_2,\dots,x_n are independent random samples drawn from Unit(0,1), then by the law of large number we have,

\frac{1}{n} \sum_{i=1}^n f(x_i) \rightarrow E(f(X)) 

Thus an estimator of I may be:

\hat{I}=\frac{1}{n} \sum_{i=1}^n f(x_i)

On a more general note, if a < b < \infty then,

I= \int_a^b f(x) dx \\
\textit{Taking $y=\frac{x-a}{b-a}$}\\
I= \int_0^1 (b-a) f \left(\frac{a+(b-a)y}{b} \right) dy  \\
=  \int_0^1 h(y) dy \quad dy \\ \text{;where } h(y)=f \left(\frac{a+(b-a)y}{b} \right)  \\
= E\left[ h(Y) | Y \sim Unif(0,1) \right]

And when b=\infty ,

I= \int_a^\infty f(x) dx \\ 
\text{Taking $y=\frac{1}{x+1}$} \\
I= - \int_1^0  f\left( \frac{1}{y} -1 \right) \frac{dy}{y^2} \\
= \int_0^1  h(y) dy \\
\text{; where } h(y) =  f\left( \frac{1}{y} -1 \right)/y^2 \\
=E(h(Y) | Y \sim Unif(0,1))

STB – 2019

  • Consider a count variable X following a Poisson distribution with parameter θ > 0, where zero count (i.e., X = 0) is not observable. We have n observations X1, . . . , Xn from this distribution. Let  denote the sample mean.

a)  Derive the quantity for which  is an unbiased estimator. 

b)  Suppose that the observed value of  is strictly greater than 1. Show that the likelihood function of θ has a unique maximiser.

The pmf of a Poisson distribution with parameter   is given by:

We know that X follows a Poisson distribution with parameter \theta   and X=0 is not observable. Under such condition the probability mass function(pmf) of X is given by:

g(x | \theta ) =Pr⁡[ X=x | X≠0] \\
\hskip{1.8cm}=Pr⁡[X=x]/Pr⁡[X≠0] \\
\hskip{0.8cm}=f(x)/(1-f(0) ) \\
\hskip{2.8cm}=\frac{θ^x e^{-θ}} {x!(1-e^{-θ} )} \quad ,x=1,2,3,…

(a) Now,

E (X)= ∑_{x=1}^∞ E\left[x g(x)\right]= ∑_{x=1}^∞E\left[x \frac{θ^x e^{-θ}} {x!(1-e^{-θ} )} \right]= ∑_{x=1}^∞\frac{θ^x e^{-θ}} {(x-1)!(1-e^{-θ} )} \\
= \frac{\theta e^{-θ}} {(1-e^{-θ} )}∑_{x=1}^∞\frac{θ^{x-1}} {(x-1)!} = \frac{\theta e^{-θ}} {(1-e^{-θ} )}∑_{x=0}^∞\frac{θ^{x}} {x!} =  \frac{\theta} {(1-e^{-θ} )}

Also, as  are random samples drawn from the distribution so,

E(\bar{X})= \frac{1}{n} \sum_{i=1}^n E(X_i) =  \frac{\theta} {(1-e^{-θ} )} \quad
\left[ i.e. \bar{X} \text{is an u.e. of }   \frac{\theta} {(1-e^{-θ} )} \right]

(b) The Likelihood function of x1,x2,…,xn is given by:

L(θ)= ∏_{x=1}^ng(x_i|\theta ) \\
\hskip{2.4cm} =∏_{x=1}^n
((1-e^{-θ} )^{-1} \frac{ θ^{x_i} e^{-θ}}{(x_i !)} \\
\hskip{1.8cm}=\frac{1}{(e^θ-1)^n}   \frac{θ^{\sum x_i }}{∏_{i=1}^{n}x_i !}

Taking ‘ln’ on both sides we get,

l(θ)=ln(L(θ))=-n log⁡(e^θ-1)+n \bar{x}  ln⁡(θ)-ln⁡(∏x_i !)

Differentiation wrt  \theta would give,

\frac{d}{d\theta}l(\theta) =  \frac{n}{1-e^{-\theta}} + \frac{n\bar{x}}{\theta} \\
\hskip{.5cm} \frac{d^2}{d\theta^2}l(\theta) = - \frac{n}{(1-e^{-\theta})^2}-\frac{n\bar{x}}{\theta^2}


\left[\frac{d}{d\theta}l(\theta) \right]_{\theta=\hat{\theta}} =0 \Rightarrow \frac{\hat{\theta}}{1-e^{-\theta}} = \bar{x} \\
 \left[\frac{d^2}{d\theta^2}l(\theta) \right]_{\theta=\hat{\theta}}= -\frac{n\bar{x}^2}{\hat{\theta}^2}-\frac{n\bar{x}}{\hat{\theta}^2}<0

Thus, the likelihood function of θ has a unique maximiser.

Ratio Estimation


In Survey Sampling, we often use information on some auxiliary variable, to improve our estimator of the finite population parameter, by giving our estimator protection against selection of bad sample. One such estimator is the ratio estimator introduced as follows:

In practice, we are often interested to estimate the ratio of the type:

R= \frac{\bar{Y}}{\bar{X}} = \frac{\sum Y_i}{\sum X_i} 

For example, in different socio-economic survey, we may be interested in per-capita expenditure on food-items, infant mortality rate, literacy rate, etc. So, the estimation of R itself is of interest to us and beside that, we can get an improved estimate of Y ̄, as follows:


Population Size N
Population U=(U_1,U_2,\dots,U_N)
Study variable Y=(Y_1,Y_2,\dots,Y_N)
Auxilliary Variable X=(X_1,X_2,\dots,X_N)
Mean of Y(study variable) \bar{Y} = \frac{1}{N} \sum Y_\alpha
Mean of X( auxilliary variable) \bar{X} = \frac{1}{N} \sum X_\alpha
Sample size n
(drawn by SRSWOR from U)
Sample Mean of Y(study variable) \bar{y} = \frac{1}{n} \sum y_i
Sample Mean of X( auxilliary variable) \bar{x} = \frac{1}{n} \sum x_i

Since,  \bar{Y} = R \bar{X}  , where R is unknown but  \bar{X}  is known, we can take

\hat{\bar{Y}}_R = \hat{R} \bar{X} = \frac{\bar{y}}{\bar{x}} \bar{X}

as an estimate of  \bar{Y}  and this estimator is called the ratio estimator of  \bar{Y}  .


\hat{\bar{Y}_R} is not an unbiased estimator of \hat{\bar{Y}} and its approximate bias is given by:

B(\hat{\bar{Y_R}}) = \frac{1}{\bar{X}} \frac{1-f}{n} \left( RS_X^2 - S_{XY} \right) 

;where \quad S_X^2= \frac{1}{N-1} \sum_{i=1}^{N} \left( X_i - \bar{X} \right)^2 , \\
\quad \quad  S_{XY}= \frac{1}{N-1} \sum_{i=1}^{N} \left( X_i - \bar{X} \right) \left(Y_i - \bar{Y} \right)


\frac{\left| B( \hat{\bar{Y}}_R) \right|}{ \sigma_{\hat{\bar{Y}}_R}} \leq \left| C.V.(\bar{x}) \right|


Mean square error of \hat{\bar{Y}_R} is given by:

MSE(\hat{\bar{Y}}_R) = \frac{1-f}{n} \left[ S_Y^2 + R^2 S_X^2 - 2 R S_{XY}   \right]


In SRSWOR, for large n, an approximation of the variance of \hat{\bar{Y}_R} is given by:

V(\hat{\bar{Y}}_R) = \frac{1-f}{n}  \frac{1}{N-1} \sum_{i=1}^{N} U_i^2 \\
\textit{; where} \quad U_i= Y_i-RX_i, \forall i=1(1)N