R语言:蒙特卡洛模拟

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

蒙特卡洛模拟的基本思路可以归纳为三步:构造问题的概率模型->从已知概率分布中抽样->建立所需的统计量。后面我会用两个例子来做详细的解释。因为本人后续的学习内容要以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=SSred,那么这个 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要素的预估的工期概率分布都呈标准正态分布,各自的平均工期、标准差以及最悲观、最可能和最乐观的估计工期如下表所示

最乐观最可能最悲观平均工期标准差分布
设计81420142正态分布
建造142332233正态分布
测试102234224正态分布

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

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

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

在这里插入图片描述
从结果中可以看出,精度还是非常高的。

结语

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

获得代码

以下是我的个人公众号,本文完整代码已上传,关注公众号回复“蒙特卡洛模拟”,即可获得,谢谢大家支持。
在这里插入图片描述

  • 23
    点赞
  • 197
    收藏
    觉得还不错? 一键收藏
  • 1
    评论
蒙特卡洛模拟在R语言中的实现可以通过以下步骤进行: 1. 构造问题的概率模型:首先,需要确定要解决的问题,并将其转化为概率模型。例如,可以使用概率密度函数来表示要求解的积分或概率分布。 2. 从已知概率分布中抽样:使用R语言中的随机数生成函数,如runif()或rnorm(),从已知的概率分布中生成随机样本。这些样本将用于模拟问题的不确定性。 3. 建立所需的统计量:根据问题的要求,使用抽样数据计算所需的统计量。例如,可以使用抽样数据计算积分的近似值或估计概率分布的参数。 下面以两个例子来说明蒙特卡洛模拟在R语言中的实现: 例1:蒙特卡洛模拟求积分 在这个例子中,我们要使用蒙特卡洛模拟来估计函数 f(x) = x^2 在区间 [0, 1] 上的积分值。首先,我们可以使用runif()函数生成在 [0, 1] 区间上均匀分布的随机样本。然后,根据这些样本计算函数值,并求取平均值,最后乘以区间长度得到积分的近似值。 例2:蒙特卡洛模拟在项目管理中的应用 在这个例子中,我们要使用蒙特卡洛模拟来估计一个项目的完成时间。首先,我们需要确定项目中各个任务的完成时间的概率分布。然后,使用R语言中的随机数生成函数,如rnorm(),从这些概率分布中生成随机样本。接下来,根据这些样本计算项目的完成时间,并重复模拟多次以得到一个分布。最后,可以使用这个分布来估计项目的完成时间的概率分布。 总结起来,蒙特卡洛模拟在R语言中的实现包括构造问题的概率模型、从已知概率分布中抽样和建立所需的统计量。通过这些步骤,可以使用蒙特卡洛模拟来解决各种问题。

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论 1
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值