View source for monte Carlo Integration

===Monte Carlo Integration - May 26, 2009===
Today's lecture completes the discussion on the Frequentists and Bayesian schools of thought and introduces '''Basic Monte Carlo Integration'''.<br><br>

====Frequentist vs Bayesian Example - Estimating Parameters====

Estimating parameters of a univariate Gaussian:

Frequentist: <math>f(x|\theta)=\frac{1}{\sqrt{2\pi}\sigma}e^{-\frac{1}{2}*(\frac{x-\mu}{\sigma})^2}</math><br>
Bayesian: <math>f(\theta|x)=\frac{f(x|\theta)f(\theta)}{f(x)}</math>

=====Frequentist Approach=====

Let <math>X^N</math> denote <math>(x_1, x_2, ..., x_n)</math>. Using the Maximum Likelihood Estimation approach for estimating parameters we get:<br>
:<math>L(X^N; \theta) = \prod_{i=1}^N \frac{1}{\sqrt{2\pi}\sigma}e^{-\frac{1}{2}(\frac{x_i- \mu} {\sigma})^2}</math>
:<math>l(X^N; \theta) = \sum_{i=1}^N -\frac{1}{2}log (2\pi) - log(\sigma) - \frac{1}{2} \left(\frac{x_i- \mu}{\sigma}\right)^2 </math>
:<math>\frac{dl}{d\mu} = \displaystyle\sum_{i=1}^N(\frac{x_i-\mu}{\sigma^2})</math>
Setting <math>\frac{dl}{d\mu} = 0</math> we get
:<math>\displaystyle\sum_{i=1}^Nx_i = \displaystyle\sum_{i=1}^N\mu</math>
:<math>\displaystyle\sum_{i=1}^Nx_i = N\mu \rightarrow \mu = \frac{1}{N}\displaystyle\sum_{i=1}^Nx_i</math><br>

=====Bayesian Approach=====

Assuming the prior is Gaussian:
:<math>P(\theta) = \frac{1}{\sqrt{2\pi}\tau}e^{-\frac{1}{2}(\frac{x-\mu_0}{\tau})^2}</math>
:<math>f(\theta|x) \propto \frac{1}{\sqrt{2\pi}\sigma}e^{-\frac{1}{2}(\frac{x_i-\mu}{\sigma})^2} * \frac{1}{\sqrt{2\pi}\tau}e^{-\frac{1}{2}(\frac{x-\mu_0}{\tau})^2}</math>
By completing the square we conclude that the posterior is Gaussian as well:
:<math>f(\theta|x)=\frac{1}{\sqrt{2\pi}\tilde{\sigma}}e^{-\frac{1}{2}(\frac{x-\tilde{\mu}}{\tilde{\sigma}})^2}</math>
Where
:<math>\tilde{\mu} = \frac{\frac{N}{\sigma^2}}{{\frac{N}{\sigma^2}}+\frac{1}{\tau^2}}\bar{x} + \frac{\frac{1}{\tau^2}}{{\frac{N}{\sigma^2}}+\frac{1}{\tau^2}}\mu_0</math>
The expectation from the posterior is different from the MLE method.
Note that <math>\displaystyle\lim_{N\to\infty}\tilde{\mu} = \bar{x}</math>. Also note that when <math>N = 0</math> we get <math>\tilde{\mu} = \mu_0</math>.

====Basic Monte Carlo Integration====

Although it is almost impossible to find a precise definition of "Monte Carlo Method", the method is widely used and has numerous descriptions in articles and monographs. As an interesting fact, the term '''Monte Carlo''' is claimed to have been first used by Ulam and von Neumann as a Los Alamos code word for the stochastic simulations they applied to building better atomic bombs.

<!-- EDITING, BACK OFF -->
We start with a simple example:
:<math>I = \displaystyle\int_a^b h(x)\,dx</math>
::<math>  = \displaystyle\int_a^b w(x)f(x)\,dx</math>
where
:<math>w(x) = h(x)(b-a)</math>
:<math>f(x) = \frac{1}{b-a} \rightarrow</math> the p.d.f. is Unif<math>(a,b)</math>
:<math>\hat{I} = E_f[w(x)] = \frac{1}{N}\displaystyle\sum_{i=1}^Nw(x_i)</math>
If <math>x_i \sim~ Unif(a,b)</math> then by the '''Law of Large Numbers''' <math>\frac{1}{N}\displaystyle\sum_{i=1}^Nw(x_i) \rightarrow \displaystyle\int w(x)f(x)\,dx = E_f[w(x)]</math>

=====Process=====
Given <math>I = \displaystyle\int^b_ah(x)\,dx</math>
# <math>\begin{matrix} w(x) = h(x)(b-a)\end{matrix}</math>
# <math> \begin{matrix} x_1, x_2, ..., x_n \sim UNIF(a,b)\end{matrix}</math>
# <math>\hat{I} = \frac{1}{N}\displaystyle\sum_{i=1}^Nw(x_i)</math>

From this we can compute other statistics, such as
# <math> SE=\frac{s}{\sqrt{N}}</math> where  <math>s=\frac{\sum_{i=1}^{N}(Y_i-\hat{I})^2 }{N-1} </math> with <math> \begin{matrix}Y_i=w(i)\end{matrix}</math>
# <math>\begin{matrix} 1-\alpha \end{matrix}</math> CI's can be estimated as <math> \hat{I}\pm Z_\frac{\alpha}{2}*SE</math>

'''Example 1'''

Find <math> E[\sqrt{x}]</math>  for  <math>\begin{matrix} f(x) = e^{-x}\end{matrix} </math>

# We need to draw from <math>\begin{matrix} f(x) = e^{-x}\end{matrix} </math>
# <math>\hat{I} = \frac{1}{N}\displaystyle\sum_{i=1}^Nw(x_i)</math>

This example can be illustrated in Matlab using the code below:

 u=rand(100,1);
 x=-log(u);
 h= x.* .5;
 mean(h)
 %The value obtained using the Monte Carlo method
 F = @ (x) sqrt (x). * exp(-x);
 quad(F,0,50)
 %The value of the real function using Matlab

'''Example 2'''
Find <math> I = \displaystyle\int^1_0h(x)\,dx, h(x) = x^3 </math>
# <math> I = x^4/4 = 1/4 </math>
# <math> w(x) = x^3*(1-0)</math>
# <math> Xi \sim~Unif(0,1)</math>
# <math>\hat{I} = \frac{1}{N}\displaystyle\sum_{i=1}^Nx_i^3</math>

This example can be illustrated in Matlab using the code below:
 x = rand (1000,1);
 mean(x.^3)

'''Example 3'''
To estimate an infinite integral
such as <math> I = \displaystyle\int^\infty_5 h(x)\,dx, h(x) = 3e^{-x} </math>
# Substitute in <math> y=\frac{1}{x-5+1} => dy=\frac{1}{-(x-4)^2}dx => dy=-y^2dx </math>
# <math> I = \displaystyle\int^1_0 \frac{3e^{-(\frac{1}{y}+4)}}{y^2}\,dy </math>
# <math> w(y) = \frac{3e^{-(\frac{1}{y}+4)}}{-y^2}(1-0)</math>
# <math> Y_i \sim~Unif(0,1)</math>
# <math>\hat{I} = \frac{1}{N}\displaystyle\sum_{i=1}^N(\frac{3e^{-(\frac{1}{y_i}+4)}}{y_i^2})</math>

===[[Importance Sampling and Monte Carlo Simulation]] - May 28, 2009 ===