<p>So someone at mathstackexchange posted the answer I was satisfied with. Here is the <a href="https://math.stackexchange.com/questions/2823164/variance-of-estimator-monte-carlo-integration/2823170#2823170">link</a></p>

<p>And here is the quoted answer.</p>

<blockquote>
  <p>This computation really has nothing to do with the Monte Carlo
  integration context. It is a general fact that if $X_i$ are iid random
  variables with variance $\sigma^2$ then $\frac{1}{n} \sum_{i=1}^n X_i$
  has variance $\sigma^2/n$. This is because of two related facts.
  First, if you sum two independent random variables, you add their
  variances. This is because:</p>
  
  <p>$$\sigma_{X+Y}^2=E[(X+Y-\mu_X-\mu_Y)^2] \\=E[((X-\mu_X)+(Y-\mu_Y))^2]
&gt; \\
&gt; =E[(X-\mu_X)^2]+E[(Y-\mu_Y)^2]+2E[(X-\mu_X)(Y-\mu_Y)] \\=\sigma_X^2+\sigma_Y^2.$$</p>
  
  <p>In the last step you use that independent random variables are
  uncorrelated. This follows from the fact that $E[XY]=E[X]E[Y]$ if $X$
  and $Y$ are independent.</p>
  
  <p>Second, if you multiply a random variable by a constant, you multiply
  its variance by the square of that constant. This is obvious from the
  definition and linearity of expectation.</p>
  
  <p>So in the context of $\frac{1}{n} \sum_{i=1}^n X_i$, the summation
  multiplies the variance by $n$ while the division by $n$ multiplies
  the variance by $1/n^2$, giving the result.</p>
</blockquote>

<p>So if we apply this $E[XY]=E[X]E[Y]$, we get the final result.</p>

CC BY-SA 4.0
2018-06-18T00:47:56.717Z
2018-06-18T00:47:56.717Z
gallickgunner
<p>Quoting from "Advanced global illumination" book:</p>

<p><a href="https://i.sstatic.net/8ovuh.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/8ovuh.png" alt="enter image description here"></a></p>

<p>The Monte Carlo estimator is defined as:</p>

<p><a href="https://i.sstatic.net/KSMPl.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/KSMPl.png" alt="enter image description here"></a></p>

<p>And its variance according to the above definition(continuous random variable) would be:</p>

<p><a href="https://i.sstatic.net/UhgvP.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/UhgvP.png" alt="enter image description here"></a></p>

<p><strong>Edit:</strong></p>

<p>An example for how the variance of an estimator G is calculated:</p>

<p>(note how the summation is written as N times the integrand and is cancelled out later.)</p>

<p><a href="https://i.sstatic.net/eLSml.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/eLSml.png" alt="enter image description here"></a></p>

CC BY-SA 4.0
2018-06-13T11:22:55.973Z
Hmm how did you remove the summation, I would appreciate though if you will tell what's wrong in my calculation rather than proving from a different way.
gallickgunner
CC BY-SA 4.0
2018-06-13T12:20:32.980Z
the summation is written as N times the f(x)/p(x) which cancels out one of the N factor in variance. I think you missed this on your last proof. Have look at the edit.
ali
CC BY-SA 4.0
2018-06-13T19:57:52.220Z
Yes I know that, but it can only be done if we are taking the  $\sum  E[f(x_i)] $. That is when we are taking the expected value inside the summation, that reduces $f(x_i)$ to $f(x)$ and then we can write it as $N*f(x)$. However in my proof I am having problem taking the Expected value of the whole square of the summation. You can't take the E operator inside of the whole square term. example,  $E [ (  \sum f(x_i) )^2]$
gallickgunner
CC BY-SA 4.0
2018-06-13T08:00:33.657Z
2018-06-13T15:54:42.260Z
2018-06-13T15:54:42.260Z
ali
<p>So I was reading this paper by Lafortune, "Mathematical Models and Monte Carlo algorithms" and in it he writes.</p>

<p>We have a function or integrand I we want to estimate given as,</p>

<p>$I = \int f(x) dx$</p>

<p>We then have a primary estimator for this as,</p>

<p>$\hat{I_p} = f(\xi)$</p>

<p>where $\xi$ is a random number generated in the interval of the integrand $I$.</p>

<p>The secondary estimator is defined as,</p>

<p>$\hat{I_{sec}} = \frac{1}{N} \sum f(\xi_i)$</p>

<p>Then there is an explanation that goes to show how the expected value of the secondary estimator is equal to the function/integrand I. i.e,</p>

<p>$E(\hat{I_{sec}}) = I$</p>

<p>All fair. The problem comes when he tries to find the variance of this secondary estimator as follows. Im posting a screenshot of the paper here since it's getting a little complex.</p>

<p><a href="https://i.sstatic.net/56MPu.png" rel="noreferrer"><img src="https://i.sstatic.net/56MPu.png" alt="enter image description here"></a></p>

<p>I don't understand how did step 2 follow from step 1 of Equation 3.5 in the variance calculation. Note that he has assumed $PDF = 1$ for now</p>

<p>I was trying to calculate the variance using a more standard approach and ignoring the multiple integrals. We know the variance of an estimator $\hat{\theta}$ is given as</p>

<p>$Var(\hat{\theta}) = E[(\hat{\theta} - E(\hat{\theta}))^2]$</p>

<p>If we apply this to above we get</p>

<p>$Var(\hat{I_{sec}}) =  E[(\hat{I_{sec}} - E(\hat{I_{sec}}))^2] $</p>

<p>$\hspace{20mm}      = E[(\hat{I_{sec}} - I)^2]$</p>

<p>$\hspace{20mm}      = E[(\hat{I_{sec}})^2 - 2* I *(\hat{I_{sec}}) + I^2 ] $</p>

<p>$\hspace{20mm}      = E[(\hat{I_{sec}})^2] - 2* I *E[\hat{I_{sec}}] + I^2 $</p>

<p>$\hspace{20mm}      = E[(\hat{I_{sec}})^2] - 2* I^2 + I^2$</p>

<p>$\hspace{20mm}      = E[(\hat{I_{sec}})^2] - I^2 $</p>

<p>$\hspace{20mm}      = E[(\frac{1}{N} \sum\limits_{i=1}^{N} f(\xi_i))^2] - I^2$</p>

<p>Now I am stuck here, I could do this,</p>

<p>$\hspace{20mm}      = E[\frac{1}{N^2} \sum\limits_{i=1}^{N} f^2(\xi_i) + \frac{1}{N^2}\sum\limits_{i\neq j}^{N} f(\xi_i) f(\xi_j) ] - I^2  \hspace{10mm} \because 
[\sum\limits_{i=1}^{N} f(x_i) ]^2 =  \sum\limits_{i=1}^{N} f^2(x_i) + \sum\limits_{i\neq j}^{N} f(\xi_i) f(\xi_j)$</p>

<p>$\hspace{20mm}      = \frac{1}{N^2} \sum\limits_{i=1}^{N} E[f^2(\xi_i)] + \frac{1}{N^2}\sum\limits_{i\neq j}^{N} E[f(\xi_i) f(\xi_j) ] - I^2$</p>

<p>$\hspace{20mm}      = \frac{1}{N^2} \sum\limits_{i=1}^{N} \int f^2(x)dx + \frac{1}{N^2}\sum\limits_{i\neq j}^{N} E[f(\xi_i) f(\xi_j) ] - I^2$</p>

<p>$\hspace{20mm} \because E[f(X)] = \int f(X) p(X) dx$</p>

<p>Note $p(X) = 1$,</p>

<p>$\hspace{20mm}      = \frac{1}{N} \int f^2(x)dx + \frac{1}{N^2}\sum\limits_{i\neq j}^{N} E[f(\xi_i) f(\xi_j) ] - I^2$</p>

<p>Don't know what to do next from here, any tips? Or did I do something wrong?</p>

<monte-carlo>
Variance of estimator (Monte Carlo Integration)
gallickgunner
# Variance of estimator (Monte Carlo Integration)

(Asked by: username_0 on 2018-06-10)

So I was reading this paper by Lafortune, "Mathematical Models and Monte Carlo algorithms" and in it he writes.

We have a function or integrand I we want to estimate given as,

$I = \int f(x) dx$

We then have a primary estimator for this as,

$\hat{I\_p} = f(\xi)$

where $\xi$ is a random number generated in the interval of the integrand $I$.

The secondary estimator is defined as,

$\hat{I\_{sec}} = \frac{1}{N} \sum f(\xi\_i)$

Then there is an explanation that goes to show how the expected value of the secondary estimator is equal to the function/integrand I. i.e,

$E(\hat{I\_{sec}}) = I$

All fair. The problem comes when he tries to find the variance of this secondary estimator as follows. Im posting a screenshot of the paper here since it's getting a little complex.

[![enter image description here](https://i.sstatic.net/56MPu.png)](https://i.sstatic.net/56MPu.png)

I don't understand how did step 2 follow from step 1 of Equation 3.5 in the variance calculation. Note that he has assumed $PDF = 1$ for now

I was trying to calculate the variance using a more standard approach and ignoring the multiple integrals. We know the variance of an estimator $\hat{\theta}$ is given as

$Var(\hat{\theta}) = E[(\hat{\theta} - E(\hat{\theta}))^2]$

If we apply this to above we get

$Var(\hat{I\_{sec}}) = E[(\hat{I\_{sec}} - E(\hat{I\_{sec}}))^2] $

$\hspace{20mm} = E[(\hat{I\_{sec}} - I)^2]$

$\hspace{20mm} = E[(\hat{I\_{sec}})^2 - 2\* I \*(\hat{I\_{sec}}) + I^2 ] $

$\hspace{20mm} = E[(\hat{I\_{sec}})^2] - 2\* I \*E[\hat{I\_{sec}}] + I^2 $

$\hspace{20mm} = E[(\hat{I\_{sec}})^2] - 2\* I^2 + I^2$

$\hspace{20mm} = E[(\hat{I\_{sec}})^2] - I^2 $

$\hspace{20mm} = E[(\frac{1}{N} \sum\limits\_{i=1}^{N} f(\xi\_i))^2] - I^2$

Now I am stuck here, I could do this,

$\hspace{20mm} = E[\frac{1}{N^2} \sum\limits\_{i=1}^{N} f^2(\xi\_i) + \frac{1}{N^2}\sum\limits\_{i\neq j}^{N} f(\xi\_i) f(\xi\_j) ] - I^2 \hspace{10mm} \because
[\sum\limits\_{i=1}^{N} f(x\_i) ]^2 = \sum\limits\_{i=1}^{N} f^2(x\_i) + \sum\limits\_{i\neq j}^{N} f(\xi\_i) f(\xi\_j)$

$\hspace{20mm} = \frac{1}{N^2} \sum\limits\_{i=1}^{N} E[f^2(\xi\_i)] + \frac{1}{N^2}\sum\limits\_{i\neq j}^{N} E[f(\xi\_i) f(\xi\_j) ] - I^2$

$\hspace{20mm} = \frac{1}{N^2} \sum\limits\_{i=1}^{N} \int f^2(x)dx + \frac{1}{N^2}\sum\limits\_{i\neq j}^{N} E[f(\xi\_i) f(\xi\_j) ] - I^2$

$\hspace{20mm} \because E[f(X)] = \int f(X) p(X) dx$

Note $p(X) = 1$,

$\hspace{20mm} = \frac{1}{N} \int f^2(x)dx + \frac{1}{N^2}\sum\limits\_{i\neq j}^{N} E[f(\xi\_i) f(\xi\_j) ] - I^2$

Don't know what to do next from here, any tips? Or did I do something wrong?

**Tags:** monte-carlo


---

## 2 answers

### Answer by username_0 on 2018-06-18. Score: 1 | Accepted answer

So someone at mathstackexchange posted the answer I was satisfied with. Here is the [link](https://math.stackexchange.com/questions/2823164/variance-of-estimator-monte-carlo-integration/2823170#2823170)

And here is the quoted answer.

> This computation really has nothing to do with the Monte Carlo
> integration context. It is a general fact that if $X\_i$ are iid random
> variables with variance $\sigma^2$ then $\frac{1}{n} \sum\_{i=1}^n X\_i$
> has variance $\sigma^2/n$. This is because of two related facts.
> First, if you sum two independent random variables, you add their
> variances. This is because:
>
> $$\sigma\_{X+Y}^2=E[(X+Y-\mu\_X-\mu\_Y)^2] \\=E[((X-\mu\_X)+(Y-\mu\_Y))^2]
> > \\
> > =E[(X-\mu\_X)^2]+E[(Y-\mu\_Y)^2]+2E[(X-\mu\_X)(Y-\mu\_Y)] \\=\sigma\_X^2+\sigma\_Y^2.$$
>
> In the last step you use that independent random variables are
> uncorrelated. This follows from the fact that $E[XY]=E[X]E[Y]$ if $X$
> and $Y$ are independent.
>
> Second, if you multiply a random variable by a constant, you multiply
> its variance by the square of that constant. This is obvious from the
> definition and linearity of expectation.
>
> So in the context of $\frac{1}{n} \sum\_{i=1}^n X\_i$, the summation
> multiplies the variance by $n$ while the division by $n$ multiplies
> the variance by $1/n^2$, giving the result.

So if we apply this $E[XY]=E[X]E[Y]$, we get the final result.


---

### Answer by username_1 on 2018-06-13. Score: 1

Quoting from "Advanced global illumination" book:

[![enter image description here](https://i.sstatic.net/8ovuh.png)](https://i.sstatic.net/8ovuh.png)

The Monte Carlo estimator is defined as:

[![enter image description here](https://i.sstatic.net/KSMPl.png)](https://i.sstatic.net/KSMPl.png)

And its variance according to the above definition(continuous random variable) would be:

[![enter image description here](https://i.sstatic.net/UhgvP.png)](https://i.sstatic.net/UhgvP.png)

**Edit:**

An example for how the variance of an estimator G is calculated:

(note how the summation is written as N times the integrand and is cancelled out later.)

[![enter image description here](https://i.sstatic.net/eLSml.png)](https://i.sstatic.net/eLSml.png)

> **Comments:**
> **username_0**: Hmm how did you remove the summation, I would appreciate though if you will tell what's wrong in my calculation rather than proving from a different way.  
> **username_1**: the summation is written as N times the f(x)/p(x) which cancels out one of the N factor in variance. I think you missed this on your last proof. Have look at the edit.  
> **username_0**: Yes I know that, but it can only be done if we are taking the  $\sum  E[f(x_i)] $. That is when we are taking the expected value inside the summation, that reduces $f(x_i)$ to $f(x)$ and then we can write it as $N*f(x)$. However in my proof I am having problem taking the Expected value of the whole square of the summation. You can't take the E operator inside of the whole square term. example,  $E [ (  \sum f(x_i) )^2]$  

---

