hello，大家好，今天带给大家的是一种计算机模拟的方法——蒙特卡洛模拟（Monte Carlo）。这是一种基于概率统计模型所衍生的一种计算机模拟的方法，而它的原理就是概率论中所涉及的“大数定律”，也就是在实验次数非常多时，频率会依概率收敛，即频率非常接近概率。

蒙特卡洛模拟的基本思路可以归纳为三步：构造问题的概率模型->从已知概率分布中抽样->建立所需的统计量。后面我会用两个例子来做详细的解释。因为本人后续的学习内容要以R语言为主，因此这两个例子我都是基于R语言来实现的。

例1 蒙特卡洛模拟求积分

I=∫01ln⁡(1+x)1+x2dxI = \int_0^1 \frac {\ln(1+x)}{1+x^2} \mathrm{d} x I=∫011+x2ln(1+x)dx
首先我们可以先画下图，看下这个所要求的积分函数是什么样子的。

f <- function(x){log(1+x)/(1+x^2)
}
x <- seq(0,1,length=50)
y <- f(x)
df <- data.frame(x,y)
ggplot(df, mapping = aes(x=x,y=y))+geom_line()

画完后发现，整个函数在[0,1]的区间内都是大于0的，那么我们就可以画一个边长为1含有这个函数的正方形，往该区域内投点，看看落在该函数内的点有多少个（即图中红色区域），这样我们就完成了第一步，构造问题的概率模型。

ggplot(df, mapping = aes(x=x,y=y))+geom_line()+geom_ribbon(aes(ymin=0, ymax=y, x = x), fill="red", alpha=0.2)+geom_hline(yintercept = c(0,1))+geom_vline(xintercept = c(0,1))

接下来，第二步，从已知概率分布中抽样，往该区域随机投点属于均匀分布，因此我们用runif()函数即可生成n个随机点。每有一个点落在红色区域，便记为一次，记k为落在红色区域内的次数。随着投点的次数增多，我们就可以根据大数定律认为kn=SredS正方形\frac kn = \frac {S_{red}}{S_{正方形}}nk=S正方形Sred，那么这个kn\frac knnk就是我们所需要的统计量。又因为这个正方形面积为1，所以kn\frac knnk就是我们要求的定积分。

MC1 <- function(n){k <- 0; x <- runif(n, 0, 1); y <- runif(n, 0, 1) #从已知概率分布中抽样for (i in 1:n){if (y[i] < f(x[i]))k <- k+1}k/n #建立所需的统计量
}

这样我们蒙特卡洛模拟就完成了，一般来说，模拟次数越多，最后的精度就越高，因此我们就模拟100000次来跟正确答案π8ln⁡2\frac \pi8 \ln 28πln2作比较。

MC1(100000)
## [1] 0.27057
pi/8*log(2)
## [1] 0.2721983

从结果中我们可以看出准确度还是非常不错的。

经过这个简单的例子，大家应该就能对蒙特卡洛模拟有个基本的印象，接下来再举一个在知乎上看到的例子（https://www.zhihu.com/question/263316961）。

例2 蒙特卡洛模拟在项目管理中的应用

现在有个项目，该项目共有三个WBS要素分别是设计、建造和测试，为了简单起见我们假设这三个WBS要素的预估的工期概率分布都呈标准正态分布，各自的平均工期、标准差以及最悲观、最可能和最乐观的估计工期如下表所示

	最乐观	最可能	最悲观	平均工期	标准差	分布
设计	8	14	20	14	2	正态分布
建造	14	23	32	23	3	正态分布
测试	10	22	34	22	4	正态分布

现在我们要通过这些信息推断整个工期的分布参数。

首先，我们可以先画出这三个的概率密度函数，有一个直观的感受。

x <- seq(7,35,length = 100)
y1 <- dnorm(x, mean = 14, sd = 2)
y2 <- dnorm(x, mean = 23, sd = 3)
y3 <- dnorm(x, mean = 22, sd = 4)
data <- data.frame(x,y1,y2,y3)
colnames(data) <- c("x","y1","y2","y3")
ggplot(data)+geom_line(aes(x=x,y=y1), color = "red")+geom_line(aes(x=x,y=y2), color = "blue")+geom_line(aes(x=x,y=y3), color = "green")+theme_classic()

紧接着构建蒙特卡洛模拟

MC2 <- function(n){y1 <- rnorm(n , mean = 14, sd = 2) #从已知概率分布中抽样y2 <- rnorm(n , mean = 23, sd = 3)y3 <- rnorm(n , mean = 22, sd = 4)y <- y1 + y2 + y3 #构造问题的概率模型result <- c(mean(y),var(y)) #建立所需的统计量，即样本均值和样本方差return(result)
}

我们根据正态分布的可加性，可以知道总工期的分布应该是均值为59，方差为29的正态分布。然后进行100000次模拟后，看看模拟的结果与真实结果相差多少。

result <- MC2(100000)
print(result)
## [1] 59.00243 29.11788
x <- seq(7,80,length = 1000)
data <- data.frame(x,y1 <- dnorm(x, mean = 14, sd = 2),dnorm(x, mean = 23, sd = 3),dnorm(x, mean = 22, sd = 4),dnorm(x, mean = result[1], sd = result[2]^0.5))
colnames(data) <- c("x","y1","y2","y3","y")
ggplot(data)+geom_line(aes(x=x,y=y1), color = "red")+geom_line(aes(x=x,y=y2), color = "blue")+geom_line(aes(x=x,y=y3), color = "green")+geom_line(aes(x=x,y=y))+theme_classic()

从结果中可以看出，精度还是非常高的。

结语

从以上两个例子不难看出，蒙特卡洛模拟在参数估计方面还是十分有效的，除了无偏估计以外，极大似然估计、渐近估计等也都可以用蒙特卡洛模拟进行计算。

获得代码

以下是我的个人公众号，本文完整代码已上传，关注公众号回复“蒙特卡洛模拟”，即可获得，谢谢大家支持。

R语言：蒙特卡洛模拟相关推荐

Logit Beta分布及其R语言随机模拟算法
Logit Beta分布及其R语言随机模拟算法 Logit Beta分布 Logit Beta分布的采样算法 Logit Beta分布是一个在广义线性模型中时常遇到的分布,通常是作为模型算法的一个中间 ...
R plot图片背景设置为透明_数据科学06 | R语言程序设计模拟和R分析器
模拟simulation ➢概率函数概率函数通常用来生成特征已知的模拟数据,以及在统计函数中计算概率值. 对于任意分布有四种基本函数: 前缀作用 d 产生随机数 r 估计概率分布的密度 p 估计累 ...
r语言计算模型的rmse_直播丨R语言与作物模型高级应用实战技术应用
随着基于过程的作物生长模型(Process-based Crop Growth Simulation Model)的发展,R语言在作物生长模型和数据分析.挖掘和可视化中发挥着越来越重要的作用.想要成为 ...
R语言解读自回归模型
R的极客理想系列文章,涵盖了R的思想,使用,工具,创新等的一系列要点,以我个人的学习和体验去诠释R的强大. R语言作为统计学一门语言,一直在小众领域闪耀着光芒.直到大数据的爆发,R语言变成了一门炙手可 ...
蒙特卡洛模拟计算风险价值VAR之R语言实现
一.解析VAR 当在分析方法中计算风险价值(VAR)时,我们需要假设金融工具的返回遵循一定的概率分布.最常用的是正态分布,这也是为什么我们通常称它为delta normal方法.要计算VAR,我们需要 ...
蒙特卡洛 pi C语言,python R 实现蒙特卡洛算法计算pi值
算法说明下面是摘抄内容,说明的非常详细,如果不懂就不要往下看了: 蒙特卡洛算法是通过概率来计算pi的值的.对于一个单位为1的正方形,以其某一个顶点为圆心,边为半径在正方形内画扇形(一个1/4的圆形的 ...
用蒙特卡洛方法计算派－python和R语言
用蒙特卡洛方法算pi-基于python和R语言最近follow了MOOC上一门python课,开始学Python.同时,买来了概率论与数理统计,准备自学一下统计.(因为被鄙视过不是统计专业却想搞数据 ...
R语言统计分布及模拟
#R语言中统计分布和模拟 #R中的各种概率统计分布 #汉文名称英文名称 R对应的名字附加参数 #β分布 beta beta shape1, shape2, ncp(偏态指数(non-central ...
R语言模拟疫情传播-gganimate包
本文用gganimate包展示模拟疫情数据本文篇幅较长,分为以下几个部分: 前言效果展示小结附录:代码前言前文<R语言模拟疫情传播-RVirusBroadcast>已经介绍了一 ...
蒙特卡洛python求解派_用蒙特卡洛方法计算派－python和R语言
标签: 用蒙特卡洛方法算pi-基于python和R语言最近follow了MOOC上一门python课,开始学Python.同时,买来了概率论与数理统计,准备自学一下统计.(因为被鄙视过不是统计专业却 ...

R语言：蒙特卡洛模拟

例1 蒙特卡洛模拟求积分

例2 蒙特卡洛模拟在项目管理中的应用

结语

获得代码

R语言：蒙特卡洛模拟相关推荐

最新文章

热门文章