R语言蒙特卡罗估计π

Posted 2021-03-19 super-yb

　初学蒙特卡洛

　

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‘)