search

GaryBAYLOR / R-code

A collection of algorithms written by myself for solving statistical problems
0 stars 0 forks source link

Use Monte Carlo to calculate the normal CDF #3

Open GaryBAYLOR opened 9 years ago

GaryBAYLOR commented 9 years ago

If we want to calculate int_0^1 1/sqrt(2pi) (-x^2/2) dx, there is no easy integral table for use to get the exact value. While, we can use Monte Carlo simulation, based on the law of large numbers. The following is R code.

MCnormal <- function(n) {
     a <- runif(n)
     b <- dnorm(a)
     c <- mean(b)
     real <- pnorm(1) - pnorm(0)
     diff <- abs(real - c)/ real
     return(list(simulated = c, real = real, diff = diff))
}

The large the n, the more precise the estimation.

> MCnormal(100)
$simulated
[1] 0.3388052

$real
[1] 0.3413447

$diff
[1] 0.007439889

> MCnormal(1000)
$simulated
[1] 0.3418248

$real
[1] 0.3413447

$diff
[1] 0.00140629

> MCnormal(1000000)
$simulated
[1] 0.3413338

$real
[1] 0.3413447

$diff
[1] 3.214687e-05

> MCnormal(1000000)
$simulated
[1] 0.3413239

$real
[1] 0.3413447

$diff
[1] 6.112997e-05

When n comes to 1e6, the difference is only .006%, which is pretty negligible.