[Math] Monte Carlo integration, expected value of the sample mean and expected value of f(x)

monte carloprobability

I am still progressing in my learning of probability and monte carlo method. I understand a basic MC estimator can be written as:

$$\bar x = { 1 \over N } \sum_{i=1}^N f(x_i) \approx E[f(x)]$$

I understand this can be seen as an approximation of $E[f(x)]$. Where it assumes that $x$ and $x_i$ are distributed according to the same probability distribution function $p$. I understand these things so far. Now what I don't understand in the book I am reading is the following thing: it computes the expected value of $\bar x$ and I am not sure how to interpret that? This is confusing because $\bar x$ is already an approximation of the mean $E[f(x)]$. So I am not too sure what the "mean of the approximation of the mean $E[f(x)]$" means. It's written in the book that:

$$E[\bar x] = E[f(x)]$$

The way I interpret this is the following and I'd like to know if this is accurate. The sample mean $\bar x$ is just an approximation. Thus if we compute the mean of this approximation (which comes back to say that you compute an infinity of these sample mean and take the average value), then the result is actually identical to the mean $E[f(x)]$. Wouldn't that be the same though than computing the sample mean $\bar x$ where the number of samples N goes to infinity?

Where I am also confused is that in the book it says that:

$$E[\bar x] = {1 \over N} \sum_{j=1}^N E[f(x_j)]$$

I just don't get where the $E[f(x_j)]$ on the right inside comes from. Why would computing the mean of $\bar x$ implies a sum of mean $E[f(x_j)]$? It's probably simple but I can't make sense of what's going on really.

Thank you.

Best Answer

It is true that if $Y_i$ have finite means, then $E[\sum a_iY_i] = \sum a_i E[Y_i]$. This is known as linearity of expectation. In your case, the author took the definition of $\overline{x}$ $$ \overline{x} = \frac{1}{N}\sum_{i=1}^N f(x_i) $$ and took the expectation of both sides, so $$ E[\overline{x}] = E[\frac{1}{N}\sum_{i=1}^N f(x_i)] = \frac{1}{N}\sum_{i=1}^N E[f(x_i)]. $$

The reason the author does this is because, although $\overline{x} \approx E[f(x)]$ is intuitively clear, this statement has no mathematical meaning. However, the formula for $E[\overline{x}]$ above provides useful information on how good the approximation is.

Related Solutions

[Math] Monte Carlo estimator

It appears like you are trying to do something like this (see slide 15). In that case, what the monte carlo estimator $\langle I \rangle$ is doing does not make sense if you think of it as calculating the $E[f(x)]$. What that estimator is doing is calculating $\int_{-\infty}^{\infty} f(x) dx$ using a Monte Carlo method that samples the $x_i$ not uniformly but according to the density function $p(x_i)$, where $p(x_i)$ is chosen to sample the places where $f(x_i)$ is large more often than where it is small, thereby making better use of limited computer resources - this is called importance sampling (as copper.hat pointed out). To calculate the expected value of $f(x_i)$ you first need to specify the distribution of the $x_i$ and then generate Monte Carlo values of $f(x_i)$ with x's drawn from that distribution.

To summarize, the Monte Carlo estimator for the average of a random function and for the integral of a deterministic function are completely different things. You can't necessarily mix them. I think your confusion is thinking $\langle I \rangle$ is for $E[f(x_i)]$, which is it not. Your intuition on how to calculate $E[f(x_i)]$ is correct.

[Math] Integration using Monte Carlo Method

Recall that if $Y$ is a random variable with density $g_Y$ and $h$ is a bounded measurable function, then $$\mathbb E[h(Y)] = \int_{\mathbb R} h(y)g_Y(y)\,\mathsf dy. $$ Moreover, if $Y\sim\mathcal U(0,1)$, then $a+(b-a)U\sim\mathcal U(a,b)$. So applying the change of variables $x=a+(b-a)u$ (with $a=0$, $b=\pi$) to the given integral, we have $$I = \int_0^1 \frac{\pi}{\sqrt{2\pi}} e^{-\frac12\sin^2 (\pi u) }\,\mathsf du=\int_0^1 h(u)\,\mathsf du, $$ with $h(u)=\sqrt{\frac\pi 2} e^{-\frac12\sin^2 (\pi u) }$. It follows then that $I=\mathbb E[h(U)]$ with $U\sim\mathcal U(0,1)$. Let $U_i$ be i.i.d. $\mathcal U(0,1)$ random variables and set $X_i=h(U_i)$, then for each positive integer $n$ we have the point estimate $$\newcommand{\overbar}[1]{\mkern 1.75mu\overline{\mkern-1.75mu#1\mkern-1.75mu}\mkern 1.75mu} \widehat{I_n} =: \overbar X_n= \frac1n \sum_{i=1}^n X_i$$ and the approximate $1-\alpha$ confidence interval $$\overbar X_n\pm t_{n-1,\alpha/2}\frac{S_n}{\sqrt n}, $$ where $$S_n = \sqrt{\frac1{n-1}\sum_{i=1}^n \left(X_i-\overbar X_n\right)^2} $$ is the sample standard deviation.

Here is some $\texttt R$ code to estimate an integral using the Monte Carlo method:

# Define "h" function
hh <-function(u) {
  return(sqrt(0.5*pi) * exp(-0.5 * sin(pi*u)^2))
}

n <- 1000 # Number of trials
alpha <- 0.05 # Confidence level
U <- runif(n) # Generate U(0,1) variates
X <- hh(U) # Compute X_i's
Xbar <- mean(X) # Compute sample mean
Sn <- sqrt(1/(n-1) * sum((X-Xbar)^2)) # Compute sample stdev
CI <- (Xbar + (c(-1,1) * (qt(1-(0.5*alpha), n-1) * Sn/sqrt(n)))) # CI bounds

# Print results
cat(sprintf("Point estimate: %f\n", Xbar))
cat(sprintf("Confidence interval: (%f, %f)\n", CI[1], CI[2]))

For reference, the value of the integral (as computed by Mathematica) is $$e^{-\frac14}\sqrt{d\frac{\pi }{2}} I_0\left(\frac{1}{4}\right) \approx 0.991393, $$ where $I_\cdot(\cdot)$ denotes the modified Bessel function of the first kind, i.e. $$I_0\left(\frac14\right) = \frac1\pi\int_0^\pi e^{\frac14\cos\theta}\,\mathsf d\theta. $$

Related Question

[Math] Monte Carlo Integration – Metropolis Algorithm using MATLAB