Estimator from subsequent Monte Carlo methods - Which one is better? 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: 


*

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

*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.
 A: 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$.
A: 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


*

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

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