monte Carlo Integration: Difference between revisions
(Initial copy) |
m (Conversion script moved page Monte Carlo Integration to monte Carlo Integration: Converting page titles to lowercase) |
||
(8 intermediate revisions by 3 users not shown) | |||
Line 14: | Line 14: | ||
:<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) = \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>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(x_i-\mu)</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 | 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 = \displaystyle\sum_{i=1}^N\mu</math> | ||
Line 32: | Line 32: | ||
====Basic Monte Carlo Integration==== | ====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 --> | <!-- EDITING, BACK OFF --> | ||
We start with a simple example: | We start with a simple example: | ||
Line 61: | Line 64: | ||
This example can be illustrated in Matlab using the code below: | This example can be illustrated in Matlab using the code below: | ||
u=rand(100,1) | u=rand(100,1); | ||
x=-log(u) | x=-log(u); | ||
h= x.* .5 | h= x.* .5; | ||
mean(h) | mean(h) | ||
%The value obtained using the Monte Carlo method | %The value obtained using the Monte Carlo method | ||
F = @ (x) sqrt (x). * exp(-x) | F = @ (x) sqrt (x). * exp(-x); | ||
quad(F,0,50) | quad(F,0,50) | ||
%The value of the real function using Matlab | %The value of the real function using Matlab | ||
Line 73: | Line 76: | ||
Find <math> I = \displaystyle\int^1_0h(x)\,dx, h(x) = x^3 </math> | Find <math> I = \displaystyle\int^1_0h(x)\,dx, h(x) = x^3 </math> | ||
# <math> I = x^4/4 = 1/4 </math> | # <math> I = x^4/4 = 1/4 </math> | ||
# <math> | # <math> w(x) = x^3*(1-0)</math> | ||
# <math> Xi \sim~Unif(0,1)</math> | # <math> Xi \sim~Unif(0,1)</math> | ||
# <math>\hat{I} = \frac{1}{N}\displaystyle\sum_{i=1}^ | # <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: | This example can be illustrated in Matlab using the code below: | ||
x = rand (1000) | x = rand (1000,1); | ||
mean(x^3) | mean(x.^3) | ||
'''Example 3''' | '''Example 3''' | ||
To estimate an infinite integral | To estimate an infinite integral | ||
such as <math> I = \displaystyle\int^\infty_5 h(x)\,dx, h(x) = 3e^{-x} </math> | 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> | # 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)}}{ | # <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> w(y) = \frac{3e^{-(\frac{1}{y}+4)}}{-y^2}(1-0)</math> | ||
# <math> Y_i \sim~Unif(0,1)</math> | # <math> Y_i \sim~Unif(0,1)</math> | ||
# <math>\hat{I} = \frac{1}{N}\displaystyle\sum_{i=1}^ | # <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 === |
Latest revision as of 09:45, 30 August 2017
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.
Frequentist vs Bayesian Example - Estimating Parameters
Estimating parameters of a univariate Gaussian:
Frequentist: [math]\displaystyle{ f(x|\theta)=\frac{1}{\sqrt{2\pi}\sigma}e^{-\frac{1}{2}*(\frac{x-\mu}{\sigma})^2} }[/math]
Bayesian: [math]\displaystyle{ f(\theta|x)=\frac{f(x|\theta)f(\theta)}{f(x)} }[/math]
Frequentist Approach
Let [math]\displaystyle{ X^N }[/math] denote [math]\displaystyle{ (x_1, x_2, ..., x_n) }[/math]. Using the Maximum Likelihood Estimation approach for estimating parameters we get:
- [math]\displaystyle{ 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]\displaystyle{ 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]\displaystyle{ \frac{dl}{d\mu} = \displaystyle\sum_{i=1}^N(\frac{x_i-\mu}{\sigma^2}) }[/math]
Setting [math]\displaystyle{ \frac{dl}{d\mu} = 0 }[/math] we get
- [math]\displaystyle{ \displaystyle\sum_{i=1}^Nx_i = \displaystyle\sum_{i=1}^N\mu }[/math]
- [math]\displaystyle{ \displaystyle\sum_{i=1}^Nx_i = N\mu \rightarrow \mu = \frac{1}{N}\displaystyle\sum_{i=1}^Nx_i }[/math]
Bayesian Approach
Assuming the prior is Gaussian:
- [math]\displaystyle{ P(\theta) = \frac{1}{\sqrt{2\pi}\tau}e^{-\frac{1}{2}(\frac{x-\mu_0}{\tau})^2} }[/math]
- [math]\displaystyle{ 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]\displaystyle{ f(\theta|x)=\frac{1}{\sqrt{2\pi}\tilde{\sigma}}e^{-\frac{1}{2}(\frac{x-\tilde{\mu}}{\tilde{\sigma}})^2} }[/math]
Where
- [math]\displaystyle{ \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{ \displaystyle\lim_{N\to\infty}\tilde{\mu} = \bar{x} }[/math]. Also note that when [math]\displaystyle{ N = 0 }[/math] we get [math]\displaystyle{ \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.
We start with a simple example:
- [math]\displaystyle{ I = \displaystyle\int_a^b h(x)\,dx }[/math]
- [math]\displaystyle{ = \displaystyle\int_a^b w(x)f(x)\,dx }[/math]
where
- [math]\displaystyle{ w(x) = h(x)(b-a) }[/math]
- [math]\displaystyle{ f(x) = \frac{1}{b-a} \rightarrow }[/math] the p.d.f. is Unif[math]\displaystyle{ (a,b) }[/math]
- [math]\displaystyle{ \hat{I} = E_f[w(x)] = \frac{1}{N}\displaystyle\sum_{i=1}^Nw(x_i) }[/math]
If [math]\displaystyle{ x_i \sim~ Unif(a,b) }[/math] then by the Law of Large Numbers [math]\displaystyle{ \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]\displaystyle{ I = \displaystyle\int^b_ah(x)\,dx }[/math]
- [math]\displaystyle{ \begin{matrix} w(x) = h(x)(b-a)\end{matrix} }[/math]
- [math]\displaystyle{ \begin{matrix} x_1, x_2, ..., x_n \sim UNIF(a,b)\end{matrix} }[/math]
- [math]\displaystyle{ \hat{I} = \frac{1}{N}\displaystyle\sum_{i=1}^Nw(x_i) }[/math]
From this we can compute other statistics, such as
- [math]\displaystyle{ SE=\frac{s}{\sqrt{N}} }[/math] where [math]\displaystyle{ s=\frac{\sum_{i=1}^{N}(Y_i-\hat{I})^2 }{N-1} }[/math] with [math]\displaystyle{ \begin{matrix}Y_i=w(i)\end{matrix} }[/math]
- [math]\displaystyle{ \begin{matrix} 1-\alpha \end{matrix} }[/math] CI's can be estimated as [math]\displaystyle{ \hat{I}\pm Z_\frac{\alpha}{2}*SE }[/math]
Example 1
Find [math]\displaystyle{ E[\sqrt{x}] }[/math] for [math]\displaystyle{ \begin{matrix} f(x) = e^{-x}\end{matrix} }[/math]
- We need to draw from [math]\displaystyle{ \begin{matrix} f(x) = e^{-x}\end{matrix} }[/math]
- [math]\displaystyle{ \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]\displaystyle{ I = \displaystyle\int^1_0h(x)\,dx, h(x) = x^3 }[/math]
- [math]\displaystyle{ I = x^4/4 = 1/4 }[/math]
- [math]\displaystyle{ w(x) = x^3*(1-0) }[/math]
- [math]\displaystyle{ Xi \sim~Unif(0,1) }[/math]
- [math]\displaystyle{ \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]\displaystyle{ I = \displaystyle\int^\infty_5 h(x)\,dx, h(x) = 3e^{-x} }[/math]
- Substitute in [math]\displaystyle{ y=\frac{1}{x-5+1} =\gt dy=\frac{1}{-(x-4)^2}dx =\gt dy=-y^2dx }[/math]
- [math]\displaystyle{ I = \displaystyle\int^1_0 \frac{3e^{-(\frac{1}{y}+4)}}{y^2}\,dy }[/math]
- [math]\displaystyle{ w(y) = \frac{3e^{-(\frac{1}{y}+4)}}{-y^2}(1-0) }[/math]
- [math]\displaystyle{ Y_i \sim~Unif(0,1) }[/math]
- [math]\displaystyle{ \hat{I} = \frac{1}{N}\displaystyle\sum_{i=1}^N(\frac{3e^{-(\frac{1}{y_i}+4)}}{y_i^2}) }[/math]