R語言蒙特卡羅估計π

　初學蒙特卡洛

　

center <- c(2.5,2.5)            #圓心
radius <- 2.5   #半徑 
distancefromcenter <- function(a){    # 點到直線的距離
  sqrt(sum((a-center)^2))
}
n <- 1000000
A <- matrix(runif(2*n,0,5),nrow = n,byrow = T)  #一行為一個點，n行
b <- apply(A, 1, distancefromcenter)   #對矩陣A按行計算點到直線的距離
num <- mean(b<radius)  #括號里為1或0，求均數相當於計算了1占n的比例 table(b<radius)
pai <- num*4  #圓面積與正方形面積之比 為π/4
pai
#畫圖
par(bg='beige') #背景色
plot(A,col='azure3',xlab = '',ylab = '',asp = 1) #asp讓x和y軸的刻度量度一樣
abline(h=0,col='goldenrod4',lty='dotdash',lwd=3) #畫上下左右的直線
abline(h=5,col='goldenrod4',lty='dotdash',lwd=3)
abline(v=5,col='goldenrod4',lty='dotdash',lwd=3)
abline(v=0,col='goldenrod4',lty='dotdash',lwd=3)
points(A[b<radius,],col='aquamarine3')
install.packages('plotrix')
library(plotrix)
draw.circle(2.5,2.5,2.5,border='coral2',lty='dashed',lwd=3) #畫圓
points(2.5,2.5,col='brown1',pch=19,cex=1.5,lwd=1.5)  #圓心

#模擬
paivector <- c()
for(i in 1:10000){
  n <- i
  A <- matrix(runif(2*n,0,5),nrow = n,byrow = T)  #一行為一個點，n行
  b <- apply(A, 1, distancefromcenter) 
  d <- subset(b,b<radius)
  num <- length(d)/length(b)
  paivector[i] <- num*4
}

install.packages('data.table')
library(data.table)
class(paivector)
pa <- data.frame(paivector)
pa1 <- data.table(pa)   #enhanced data frame
pai <- pa1[,id:=seq(0,9999)] #可以直接加添加一列，但是pa不行
pai <- pa1[,error := abs(pi-paivector)]
names(pai) <- c('guess','id','error')

library(ggplot2)
ggplot(pai,aes(x=id,y=error))+geom_line(color='#388E8E')+ggtitle('error')+xlab('sample size')+ylab('error')