Monte Carlo估计pi值
圆和外接正方形的面积比,是$ \frac{\pi r^2}{(2r)^2} = \frac{\pi}{4}$.
通过这一比值,可以使用蒙特-卡罗方法来估计Pi,这是Monte Carlo方法的最经典的一个例子。
getPI <- function(N) {
x <- runif(N)
y <- runif(N)
hits <- sum(sqrt(x^2+y^2) < 1)
pi <- 4*hits/N
return(pi)
}
options(digits=15)
set.seed(12345)
n <- c(seq(from=100, to=1000, by=100),
seq(from=1000, to=10000, by=1000)
)
require(ggplot2)
p <- ggplot() +
aes(x=n, y=sapply(n, getPI)) +
geom_point() +
geom_hline(aes(yintercept=pi, colour="red")) +
xlab("")+ ylab("") +
theme(legend.position="none")
print(p)