R语言学习笔记-蒙特卡洛模拟方法

前言

 

二十世纪最伟大的10大算法之一，数学家冯·诺伊曼用驰名世界的赌城—摩纳哥的Monte Carlo—来命名这种方法，为它蒙上了一层神秘色彩。所谓蒙特卡洛方法，简单地说就是将问题转化成一个概率问题.并用计算机模拟产生一堆随机数据，之后就是对随机数据的统计工作了！ 

蒙特卡洛模拟方法=建立概率模型+计算机模拟+数理统计

什么是MCMC，什么时候使用它？

Markov Chain Monte Carlo，以下简称MCMC，MCMC只是蒙特卡洛方法的一种，尽管可以将许多其他常用方法看作是MCMC的简单特例。

实例分析

1.应用蒙特卡洛模拟方法计算Π值

1）基本原理

P=圆的面积/正方形的面积

P=Π/4

2）用计算机模拟，产生0—1之间的二维的随机数

在正方形的内部产生10000个点，统计好多个点（红色）落在圆中，算得一个比值（相当于古典概率）

3）对产生的随机数样本进行统计分析，算出pi值，并计算误差

这里我们看计算公式

P=Π/4，如果将比值乘上4，就可以得到pi值，即：

那么问题来了，怎么让计算机知道哪些点在圆中，哪些点在圆外？？？

计算该点到圆心的距离，然后与R（半径）比较大小即可。

2.代码实现

想对你说的话都在代码里面

# 生成10000个点
# 圆心
n = 10000
center = c(2.5,2.5)
# 半径
radius = 2.5
# distance From Center
dfc = function(a){
  sqrt(sum((center - a)^2))
  #　欧式距离计算公式
}
# 设生成n个点
n = 10000
# 生成2n个区间是[0,5]的均匀分布产生的随机数作为A矩阵的元素

# A 是n行2列
A = matrix(runif(n*2,min = 0,max = 5),nrow = n,ncol = 2,byrow = T)
# 对矩阵的每行的元素使用求距离函数，得到向量B
# B 是存放这点到圆心的距离，注意B是向量
B = apply(A, 1, dfc)
# 求向量b中的元素值（距离）小于半径的个数

num = mean(B < radius)
#　这是一个逻辑表达式的结果要么为０，要么为１
# 所以对这个向量求均值，就得到了比例
print(num)


# 绘图
# 设置背景色
par(bg = 'beige')
# 回值模型生成的所有点，让横纵轴不显示名称
# asp = 1的色湖之是为了不使后面的draw.circle函数绘制元的时候出现偏差
plot(A,col='azure3',asp=1,xlab = '',ylab = '',main = 'MC method application')

# 添加四条线,勾勒出正方形
abline(h=0,col='goldenrod4',lty = 'dotdash',lwd=3)
abline(h=5,col='goldenrod4',lty = 'dotdash',lwd=3)
abline(v=0,col='goldenrod4',lty = 'dotdash',lwd=3)
abline(v=5,col='goldenrod4',lty = 'dotdash',lwd=3)


# 绘制圆内和边缘上的点
points(A[B<radius,],col='aquamarine3')
#导入绘圆的包
library(plotrix)
draw.circle(2.5,2.5,2.5,border='coral2',lty='dashed',lwd=3)
# 添加圆心
points(x=2.5,y=2.5,col='brown1',pch=19,cex=1.5,lwd=1.5)
# 美化
points(x=2.5,y=2.5,col='olivedrab3',pch=11,cex=3,lwd=3)

# 模拟多次，用循环
# 创建向量 
piVec = c()
# 循环2000次
for(i in 1:2000){
  n = i
  # i个区间是在[0,5]的均匀分布产生的随机数构成n行2列的矩阵Ａ
  A = matrix(runif(n*2,min=0,max=5),nrow = n,ncol=2,byrow=T)
  # 求点到圆心的距离
  # B 是存放这点到圆心的距离，注意B是向量
  b = apply(A, 1, dfc)
  # 求向量b中的元素值（距离）小于半径的个数
  
  d = subset(b,b<radius)
  # 得到比例
  num = length(d)/length(b)
  piVec[i]=num*4
}
  #　这是一个逻辑表达式的结果要么为０，要么为１
  # 所以对这个向量求均值，就得到了比例
library(data.table)
  # 将piVec转为数据框
Pi = data.frame(piVec)
Pi = data.table(Pi)
  # ind列中储存的是模拟次数
Pi[,ind:=seq(0,1999)]
  # error列是误差
Pi[,error:=abs(pi-piVec)]
Pi = data.frame(Pi)
  # pi数据捐给的三列分别为序号、实际值、误差
names(Pi) = c('guess','ind','error')
  # 画出误差
library(ggplot2)
  # x轴是模拟次数，y是误差
ggplot(Pi,aes(x=ind,y=error))+
  geom_line(colour='#388E8E')+
  ggtitle('Error')+
  xlab('Sample Size')+
  ylab('Error')