I want to approximate a probability via Monte Carlo simulation.

This probability depends upon some randomness that happens in a system. To be more specific the system can have two states. One with 20% probability and one with 80%.

What I am currently doing is the following simulation:

  • Generate a random state of the system (20% to be state 1, 80% to be state 2).
  • Check if the event happens.
  • Update the frequency and estimated probability accordingly.

This seems to be working fine. For example I get good results for N=1000 trials.

Now let's assume that for some reason I don't have access to the above method. Instead I want to do either of the following:

  1. Run 500 tests for state 2 of the system and get a probability estimate P2. Run 500 tests for state 1 and get another probability estimate P1. Get an estimator of the whole system's probability as: P = 0.2*P1 + 0.8*P2.
  2. Alternatively, run 800 tests for state 2 and get a P2. Then run 200 tests for state 1 and get a P1. Then calculate the whole system's probability as: P = 0.5*P1 + 0.5*P2.

My question is which of the two methods is a better or worse alternative to the main method (which if I am not mistaken is the proper Monte Carlo estimation). I try to keep the number of runs the same (1000 in total) so as to not increase the computational load.


2 Answers 2


First thing : the second method doesn't give you an unbiased estimate of the probability of failure of the system.

Let's say the probability of system's failure in state $1$ is $p_1$ and in state $2$ is $p_2$. Then the probability of failure of whole system is $p = 0.2 p_1 + 0.8 p_2$. If you have an unbiased estimate of $p_1$ (say $E_1$) and an unbiased estimate of $p_2$ (say $E_2$) than the only convex combination of them that will give you an unbiased estimate of $p$ is $0.2 E_1 + 0.8 E_2$.

Now you can think, how many samples should you draw from state $1$ and $2$. One possible answer is to choose the number that will minimize the variance of your estimator. You can sample $A$ samples from state $1$, and get $a$ failures. Then you can sample $B = 1000-A$ samples from state $2$ and get $b$ failures. The unbiased estimator you get is : $$ 0.2 \frac{a}{A} + 0.8 \frac{b}{B} $$ The variance of this estimator is equal to: $$ 0.04 \frac{p_1(1-p_1)}{A} + 0.64 \frac{p_2(1-p_2)}{B} $$ So the optimal numbers $A$, $B$ depend on $p_1, p_2$. If for example you assume that they are equal, then to minimize variance of the estimator you want to have $B = 16A$, which would be best approximated by $A=59$.

  • $\begingroup$ Thank you.. I don't fully understand how the variance is calculated? Is it p1*p2 on the nominator? $\endgroup$
    – KPr
    Sep 20, 2016 at 12:16

Both your solutions run stratified sampling, meaning you remove the variance due to the random choice of state. Your resulting estimator is a weighted average of both estimators, with variance

  1. $0.2^2 \frac{P_1(1-P_1)}{500}+0.8^2\frac{P_2(1-P_2)}{500}=0.8\,10^{-4}\,P_1(1-P_1)+1.28\,10^{-3}\,P_2(1-P_2)$
  2. $0.5^2\left[\frac{P_1(1-P_1)}{200}+\frac{P_2(1-P_2)}{800}\right]=1.25\,10^{-3}\,P_1(1-P_1)+3.125\,10^ {-4}\,P_2(1-P_2)$

Hence it all depends on the relative values of $P_1$ and $P_2$.

By comparison, your original estimator can be written as $$10^{-3}[\hat{N}\hat{P_1}+(10^3-\hat{N})\hat{P_2}]$$ with $\hat{N}$ a Binomial Bin$(10^3,.2)$, $\hat{N}\hat{P_1}\sim\text{Bin}(\hat{N},P_1)$ and $(10^3-\hat{N})\hat{P_2}\sim\text{Bin}(10^3-\hat{N},P_2)$. Its variance is equal to $$10^{-6}\left[.2\times.8\,10^3\,[P_1+P_2]+\mathbb{E}\left\{\hat{N}^2\frac{P_1(1-P_1)}{\hat{N}}+(10^3-\hat{N})^2\frac{P_2(1-P_2)}{(10^3-\hat{N})}\right\}\right]$$by the decomposition $$\text{var}(\delta)=\text{var}\mathbb{E}[\delta|\hat{N}]+\mathbb{E}[\text{var}(\delta|\hat{N})]$$that is $$.16\,10^{-3}[P_1+P_2]+.2\,10^{-3}\,P_1(1-P_1)+.8\,10^{-3}\,P_2(1-P_2)$$

  • $\begingroup$ Why is there Pi*(1-Pi) in the nominators, I am not getting this part? $\endgroup$
    – KPr
    Sep 21, 2016 at 18:51
  • $\begingroup$ Because the variance of a Binomial $B(N,p)$ is $Np(1-p)$... $\endgroup$
    – Xi'an
    Sep 21, 2016 at 19:05

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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