• 设为首页
  • 点击收藏
  • 手机版
    手机扫一扫访问
    迪恩网络手机版
  • 关注官方公众号
    微信扫一扫关注
    公众号

R语言蒙特卡罗估计π

原作者: [db:作者] 来自: [db:来源] 收藏 邀请

 初学蒙特卡洛

 

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

  


鲜花

握手

雷人

路过

鸡蛋
该文章已有0人参与评论

请发表评论

全部评论

专题导读
热门推荐
热门话题
阅读排行榜

扫描微信二维码

查看手机版网站

随时了解更新最新资讯

139-2527-9053

在线客服(服务时间 9:00~18:00)

在线QQ客服
地址:深圳市南山区西丽大学城创智工业园
电邮:jeky_zhao#qq.com
移动电话:139-2527-9053

Powered by 互联科技 X3.4© 2001-2213 极客世界.|Sitemap