蒙特卡洛方法简介 _数据科学

介绍数据分析过程中出现的两大类数字问题是优化和积分问题。并非总是能够分析性地计算与给定模型关联的估计量, 因此经常导致我们考虑数值解。避免该问题的一种方法是使用仿真。蒙特卡洛估计是指从概率分布中模拟假设抽取, 以计算该分布的大量数量。
蒙特卡洛(Monte Carlo)的基本思想是将积分写成某个概率分布的期望值, 然后使用矩估计器($ E [g(X)] \ approx \ overline {g(X)} = \ dfrac {1} {n} \ sum g(X_ {i})$)。
如果我们有一个连续函数$ g(\ theta)$, 并且想要对区间(a, b)进行积分, 则可以将积分重写为均匀分布$ U \ sim U [a, b]的期望值$, 即：

文章图片
使用矩估计器的方法, 我们的积分近似为：

文章图片
哪里

文章图片
是来自均匀分布的模拟值。
示例1：指数积分逼近
1)给定一个函数$ f(x)= e ^ {x} $, 区间[3, 5]中的积分为：

文章图片
积分的蒙特卡洛近似为：

#Declaring the desired function f = function(x){return(exp(x))} #Declaring the absolute error function error = function(x, y){return(abs(x-y))} #The actual integral answer ans = exp(5)-exp(3)set.seed(6971) #number of iterations n = 10^2 #simulated uniform data x= runif(n, 3, 5) #MonteCarlo approximation MCa= (5-3)*mean(f(x)) #Approximation error e = error(ans, MCa)rest =data.frame(n = n, MCapprox = MCa, error = e) set.seed(6971) for(k in 3:6){ n = 10^k x = runif(n, 3, 5) mca = (5-3)*mean(f(x)) rest = rbind(rest, c(n, mca, error(ans, mca) ) ) }kable(rest, digits = 5, align = 'c', caption = "Integral Monte Carlo approximation results", col.names =c("Number of simulations", "Monte Carlo approximation", "Error approximation"))

广义蒙特卡洛逼近
在一般情况下, 给定分布f的积分近似为：

文章图片
可以通过以下步骤描述构建$ \ widehat {I} $的算法：
1)产生

文章图片
从f分布
2)计算：

文章图片
3)获得样本均值：$$ \ overline {I} = \ dfrac {1} {n} \ sum_ {k = 1} ^ {n} \ dfrac {g(\ theta_k)} {f(\ theta_k)} $$
在下一个块中, 将提供简单的蒙特卡洛逼近函数以说明算法的工作原理, 其中a和b是均匀密度参数, n是所需模拟的数量, f是我们要集成的函数。

# The simple Monte Carlo function MCaf = function(n, a, b, f){ x = runif(n, a, b) MCa = (b-a)*mean(f(x)) return(MCa) }

贝叶斯数据分析中的蒙特卡洛方法
贝叶斯数据分析的主要思想是使用贝叶斯推理方法拟合模型(例如回归模型或时间序列模型)。我们假设感兴趣的参数具有理论分布, 使用观察数据的分布(似然性)更新该分布(后验), 并使用贝叶斯定理更新有关我们参数的先前或外部信息(优先分布)。
$$ p(\ theta / X)\ text {} \ alpha \ text {} p(X / \ theta)p(\ theta)$$其中：
$ p(\ theta / X)$是参数后验分布
【蒙特卡洛方法简介】$ p(X / \ theta)$是观测数据的采样分布(似然)。
$ p(\ theta)$是参数先验分布。
贝叶斯方法的主要问题是估计后验分布。马尔可夫链蒙特卡洛方法(mcmc)生成后验分布的样本, 并使用蒙特卡洛方法近似估计期望值, 概率或分位数。
在接下来的两节中, 我们提供两个示例以近似理论分布的概率和分位数。上面介绍的过程是贝叶斯方法中常用的方法。
示例2：伽马分布的概率近似
让我们假设我们要计算一个随机变量$ \ theta $在零和5之间的概率$ P(0 < \ theta < 5)$, 其中$ \ theta $具有伽玛分布, 参数a = 2和b = 1 / 3($ \ theta \ sim Gamma(a = 2, b = 1/3)$), 因此概率为：

文章图片
如果$ \ theta $属于区间[0, 5], 则$ I _ {[0, 5]}(\ theta)= 1 $。蒙特卡罗近似法的思想是计算属于该区间的观测数。间隔[0, 5], 然后将其除以模拟数据的总和。

set.seed(6972) #number of iterations n = 10^2 #simulated uniform data x= rgamma(n, shape = 2, 1/3) #MonteCarlo approximation MCa= mean(x < = 5) #Approximation error e = error(pgamma(5, 2, 1/3), MCa)rest =data.frame(n = n, MCapprox = MCa, error = e)for(k in 3:6){ n = 10^k x= rgamma(n, shape = 2, 1/3) mca= mean(x < = 5) rest = rbind(rest, c(n, mca, error(pgamma(5, 2, 1/3), mca) ) ) }kable(rest, digits = 5, align = 'c', caption = "Probability Monte Carlo approximation results", col.names =c("Number of simulations", "Monte Carlo approximation", "Error approximation"))

示例3：正态分布的分位数逼近
假设我们要计算具有正态分布且参数$ \ mu = 20 $和$ \ sigma = 3 $的随机变量$ \ theta $的0.95分位数($ \ theta \ sim normal(\ mu = 20, \ sigma ^ {2} = 9)$), 因此0.95分位数为：

文章图片
主要思想是找到最大概率等于或小于0.95的样本值, 使用模拟数据的Quantile()函数估算蒙特卡洛分位数近似值。

set.seed(6973) #number of iterations n = 10^2 #simulated uniform data x= rnorm(n, 20, 3) #MonteCarlo approximation MCa= quantile(x, 0.95) #Approximation error e = error(qnorm(0.95, 20, 3), MCa)rest =data.frame(n = n, MCapprox = MCa, error = e)for(k in 3:6){ n = 10^k x= rnorm(n, 20, 3) mca= quantile(x, 0.95) rest = rbind(rest, c(n, mca, error(qnorm(0.95, 20, 3), mca) ) ) }kable(rest, digits = 5, align = 'c', caption = "Quantile Monte Carlo approximation results", col.names =c("Number of simulations", "Monte Carlo approximation", "Error approximation"), row.names = FALSE)

讨论与结论
蒙特卡洛逼近法为积分逼近提供了另一种工具, 并且是贝叶斯推断方法中的重要工具, 尤其是当我们使用复杂模型时。在我们所有的三个示例中, 蒙特卡洛方法都提供了极好的逼近度, 但要使逼近误差接近零, 就需要大量的仿真。
参考文献
1)在R, Springer 2004, Christian P.Robert和George Casella中介绍了蒙特卡洛方法。
2)《马尔可夫链手册》蒙特卡罗, 查普曼和霍尔, 史蒂夫·布鲁克斯, 安德鲁·盖尔曼, 加林·L·琼斯和小李·孟的手册。
3)数理统计概论, 皮尔逊, 罗伯特·霍格, 约瑟夫·麦基恩和艾伦·克雷格。
4)统计推论, 一种综合方法, 查普曼和霍尔, Helio S. Migon, Dani Gamerman, Francisco Louzada。