R语言 Simulation

Question

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

目錄 Show

例1 蒙特卡洛模拟求积分
例2 蒙特卡洛模拟在项目管理中的应用
获得代码

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

例1 蒙特卡洛模拟求积分

I = ∫ 0 1 ln ⁡ ( 1 + x ) 1 + x 2 d x I = \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为落在红色区域内的次数。随着投点的次数增多，我们就可以根据大数定律认为 k n = S r e d S 正方形 \frac kn = \frac {S_{red}}{S_{正方形}} nk=S正方形Sred，那么这个 k n \frac kn nk就是我们所需要的统计量。又因为这个正方形面积为1，所以 k n \frac kn nk就是我们要求的定积分。

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次来跟正确答案 π 8 ln ⁡ 2 \frac \pi8 \ln 2 8π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()