Monte-Carlo
  • Introduction
  • AIS
  • Generate R.V.
    • Special Distribution
    • Copulas
    • Minimum of Two Exponential
  • Gibbs
    • Comparing two groups
    • Linear Regression
    • Simulation of Exp-Abs-xy
  • Markov Chain
  • MC Approximation
  • MC Integration
    • Rao-Blackwellization
  • MC Optimization
  • MCMC
    • MCMC diagnostics
  • Metropolis-Hastings
    • Metropolis
    • Independent MH
    • Random Walk MH
    • ARMS MH
  • PBMCMC
  • RJMCMC
  • Diagnosing Convergence
  • SMCTC
  • Tempering
    • Parallel Tempering
  • Misc
    • R vs. Julia
  • References
Powered by GitBook
On this page
  • Classical Monte Carlo Integration
  • Toy Example
  • Polar Simulation

MC Integration

PreviousMC ApproximationNextRao-Blackwellization

Last updated 6 years ago

Classical Monte Carlo Integration

introduces the classical Monte Carlo integration:

Toy Example

根据大数律有

并且由中心极限定理有

另外,注意到

则可以直接利用上一节的结论。

## estimate pi
##
## I = \int H(x, y)dxdy
## where H(x, y) = 1 when x^2+y^2 <= 1;
## otherwise H(x, y) = 0

## volume
V = 4

n = 100000
x = runif(n, -1, 1)
y = runif(n, -1, 1)
H = x^2+y^2
H[H<=1] = 1
H[H>1] = 0
I = V* mean(H)
cat("I = ", I, "\n")

## n = 100, I = 2.96
## n = 1000, I = 3.22
## n = 10000, I = 3.1536
## n = 100000, I = 3.14504

Polar Simulation

The following Julia program can be used to do polar simulation.

using Statistics

function polarsim(x::Array)
    while true
        phi1 = rand()*2*pi
        phi2 = rand()*pi - pi/2
        u = rand()
        xdotxi = x[1]*cos(phi1) + x[2]*sin(phi1)*cos(phi2) + x[3]*sin(phi1)*sin(phi2)
        if u <= exp(xdotxi - sum(x.^2)/2)
            return(randn() + xdotxi)
        end
    end
end

# approximate E^pi
function Epi(m, x = [0.1, 1.2, -0.7])
    rhos = ones(m)
    for i = 1:m
        rhos[i] = 1/(2*polarsim(x)^2+3)
    end
    return(mean(rhos))
end

# example
Epi(10)

为了计算积分 我们通常采用MC模拟:

计算区域的volume:

近似:

其中

I=∫Dg(x)dx=∫DVg(x)1Vdx=Ef[Vg(x)],I = \int_Dg(\mathbf x)d\mathbf x = \int_D Vg(\mathbf x)\frac 1V d\mathbf x = {\mathbb E}_f[Vg(\mathbf x)],I=∫D​g(x)dx=∫D​Vg(x)V1​dx=Ef​[Vg(x)],

Note: Unless otherwise stated, the algorithms and the corresponding screenshots are adopted from .

Robert and Casella (2013)
Robert and Casella (2013)