Estimator from subsequent Monte Carlo methods - Which one is better?

Question

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.

Answer 1

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

Answer 2

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