第二话 统计计算之蒙特卡洛积分和方差缩减技术（未完待续）
统计计算
第一话 统计计算之随机变量生成方法
第二话 统计计算之蒙特卡洛积分和方差缩减技术
提示：写完文章后，目录可以自动生成，如何生成可参考右边的帮助文档
蒙特卡洛积分和方差缩减技术
- 统计计算
- 前言
- 1 蒙特卡洛积分
- 1.1 基本原理
- 1.2 [0，1]区间上的简单蒙特卡洛估计
- 1.3 [a,b]区间上的MC估计
- 1.3.1 变量转换
- 1.3.2 直接抽样
- 1.4 无界区间上的MC估计
- 1.5 hit-or-miss 方法
- 2 方差缩减
- 2.1 对偶变量法（寻找一对负相关的变量，而且要保证 g g g函数是单调的）
- 2.1.1 对偶变量思想
- 2.2 控制变量法
- 2.2.1 普通控制变量法
- 2.2.2 对偶变量作为控制变量
前言
蒙特卡罗积分是一种基于随机抽样的统计模拟的方法，本博客重点介绍几种常见的蒙特卡罗积分的方法以及几种常见的方差缩减的技术，在实际中有着广泛的应用。
1 蒙特卡洛积分
1.1 基本原理
假设 g ( x ) g(x) g(x)是一个函数，并且我们想要计算 ∫ a b g ( x ) d x \int_{a}^{b} g(x) d x ∫abg(x)dx。假设 X X X是一个随机变量，密度为 f ( x ) f(x) f(x)，则随机变量 Y = g ( X ) Y=g(X) Y=g(X)的数学期望为：
E [ g ( X ) ] = ∫ − ∞ ∞ g ( x ) f ( x ) d x E[g(X)]=\int_{-\infty}^{\infty} g(x) f(x) d x E[g(X)]=∫−∞∞g(x)f(x)dx
假设随机样本是从 X X X的分布中抽取的，那么 E [ g ( X ) ] E[g(X)] E[g(X)]的无偏估计就是样本均值。这是蒙特拉罗积分的理论基础。
1.2 [0，1]区间上的简单蒙特卡洛估计
假设现在考虑一个积分问题：
θ = ∫ 0 1 g ( x ) d x \theta=\int_{0}^{1} g(x) d x θ=∫01g(x)dx
由于积分区间是0到1，因此我们假设随机变量 X 1 , … , X m X_{1}, \ldots, X_{m} X1,…,Xm是来自均匀分布 U ( 0 , 1 ) U(0,1) U(0,1)的样本，根据1.1节的基本原理，则有：
θ ^ = g m ( X ) ‾ = 1 m ∑ i = 1 m g ( X i ) \hat{\theta}=\overline{g_{m}(X)}=\frac{1}{m} \sum_{i=1}^{m} g\left(X_{i}\right) θ^=gm(X)=m1i=1∑mg(Xi)
且 θ ^ \hat{\theta} θ^以概率1收敛于 E [ g ( X ) ] = θ E[g(X)]=\theta E[g(X)]=θ，因此 ∫ 0 1 g ( x ) d x \int_{0}^{1} g(x) d x ∫01g(x)dx的简单MC估计就是 g m ( X ) ‾ \overline{g_{m}(X)} gm(X)。
例1：
m <- 10000
x <- runif(m)
theta.hat <- mean(exp(-x))
print(theta.hat)
print(1 - exp(-1))
1.3 [a,b]区间上的MC估计
本节着重讨论 ∫ a b g ( t ) d t \int_{a}^{b} g(t) d t ∫abg(t)dt型的积分，针对这样的积分主要有两种思路：
1.3.1 变量转换
该方法通过变量转化的思想，将积分区间由[a,b]转换到[0,1]区间上，再根据1.2节[0,1]区间上的积分思路来求解。
线性变换如下： y = ( t − a ) / ( b − a ) y=(t-a) /(b-a) y=(t−a)/(b−a)，并且 d y = ( 1 / ( b − a ) ) d t d y=(1 /(b-a)) d t dy=(1/(b−a))dt，替换原来的 t t t，得到如下表达式：
∫ a b g ( t ) d t = ∫ 0 1 g ( y ( b − a ) + a ) ( b − a ) d y \int_{a}^{b} g(t) d t=\int_{0}^{1} g(y(b-a)+a)(b-a) d y ∫abg(t)dt=∫01g(y(b−a)+a)(b−a)dy
另 Y ∼ U ( 0 , 1 ) Y \sim U(0,1) Y∼U(0,1)，根据1.2节即可得到积分估计。
1.3.2 直接抽样
假设 t t t是来自均匀分布 U ( a , b ) U(a,b) U(a,b)，则构造均匀分布 U ( a , b ) U(a,b) U(a,b)的概率密度函数，则表达式如下：
∫ a b g ( t ) d t = ( b − a ) ∫ a b g ( t ) 1 b − a d t = ( b − a ) ∫ a b g ( t ) f ( t ) d t \int_{a}^{b} g(t) d t=(b-a) \int_{a}^{b} g(t) \frac{1}{b-a} d t =(b-a) \int_{a}^{b} g(t) f(t) d t ∫abg(t)dt=(b−a)∫abg(t)b−a1dt=(b−a)∫abg(t)f(t)dt
基于上式可知， ∫ a b g ( t ) d t \int_{a}^{b} g(t) d t ∫abg(t)dt其实就是 b − a b-a b−a倍的 E [ g ( Y ) ] E[g(Y)] E[g(Y)],而 Y ∼ U ( a , b ) Y \sim U(a,b) Y∼U(a,b).
总结简单的蒙特卡洛积分 ∫ a b g ( t ) d t \int_{a}^{b} g(t) d t ∫abg(t)dt估计的步骤如下：
例2：
m <- 10000
x <- runif(m, min=2, max=4)
theta.hat <- mean(exp(-x)) * 2
print(theta.hat)
print(exp(-2) - exp(-4))
1.4 无界区间上的MC估计
假如使用蒙特卡罗估计标准正态分布的累积密度函数：
Φ ( x ) = ∫ − ∞ x 1 2 π e − t 2 / 2 d t \Phi(x)=\int_{-\infty}^{x} \frac{1}{\sqrt{2 \pi}} e^{-t^{2} / 2} d t Φ(x)=∫−∞x2π 1e−t2/2dt
思路：由于积分限是无界的，因此我们不能直接使用上述1.3节的简单MC估计，由于标准正态分布是关于0对称的，因此，我们可以将这个问题划分为两个部分， x > = 0 和 x < 0 x>=0和x<0 x>=0和x<0，当 x > = 0 x>=0 x>=0时，问题其实就转化为有界积分了 θ = ∫ 0 x e − t 2 / 2 d t \theta=\int_{0}^{x} e^{-t^{2} / 2} d t θ=∫0xe−t2/2dt，可以根据1.3节的内容来对上述积分进行估计，我们可以生成 U ( 0 , x ) U(0,x) U(0,x)上的随机数。
x <- seq(.1, 2.5, length = 10)
m <- 10000
u <- runif(m)
cdf <- numeric(length(x))
for (i in 1:length(x)) {g <- x[i] * exp(-(u * x[i])^2 / 2)
cdf[i] <- mean(g) / sqrt(2 * pi) + 0.5
}
1.5 hit-or-miss 方法
命题：定义 I ( ⋅ ) I(·) I(⋅)是示性函数，并且 Z ∼ N ( 0 , 1 ) Z \sim N(0,1) Z∼N(0,1),对于任意的常数 x x x,有 E [ I ( Z ≤ x ) ] = P ( Z ≤ x ) = Φ ( x ) E[I(Z \leq x)]=P(Z \leq x)=\Phi(x) E[I(Z≤x)]=P(Z≤x)=Φ(x)，标准正态分布的cdf。
基于上述命题，假设抽取一组来自于标准正态分布的随机样本 z 1 , . . . , z m z_1,...,z_m z1,...,zm，那么则有：
Φ ( x ) ^ = 1 m ∑ i = 1 m I ( z i ≤ x ) \widehat{\Phi(x)}=\frac{1}{m} \sum_{i=1}^{m} I\left(z_{i} \leq x\right) Φ(x) =m1i=1∑mI(zi≤x)
收敛于 E [ I ( Z ≤ x ) ] = P ( Z ≤ x ) = Φ ( x ) E[I(Z \leq x)]=P(Z \leq x)=\Phi(x) E[I(Z≤x)]=P(Z≤x)=Φ(x)。
x <- seq(.1, 2.5, length = 10)
m <- 10000
z <- rnorm(m)
dim(x) <- length(x)
p <- apply(x, MARGIN = 1,
FUN = function(x, z) {mean(z < x)}, z = z)
hit-or-miss方法来估计 F ( x ) = ∫ − ∞ x f ( t ) d t F(x)=\int_{-\infty}^{x} f(t) d t F(x)=∫−∞xf(t)dt的步骤如下：
2 方差缩减
简单的MC积分是对函数求期望，本节主要考虑几种重要的方差缩减的方法来减少估计量的方差。
2.1 对偶变量法（寻找一对负相关的变量，而且要保证 g g g函数是单调的）
2.1.1 对偶变量思想
2.2 控制变量法
2.2.1 普通控制变量法
主要思想：寻找与 g ( X ) g(X) g(X)紧密相关的随机变量 f ( X ) f(X) f(X)作为控制变量，且 f ( X ) f(X) f(X)的均值是已知的。
2.2.2 对偶变量作为控制变量
第二话 统计计算之蒙特卡洛积分和方差缩减技术（未完待续）相关推荐
- 网络爬虫_第二章_提取_第四单元_BeautifulSoup库入门(未完待续)
1.Beautiful Soup库的基本元素 Beautiful Soup库的理解 <p>..</p>: 标签:Tag<p> p指的是标签的名字, Name成对出现 ...
- 热力学与统计物理全揽（未完待续）
三种分布及其热力学量 量子态密度计算 MB分布 BE分布 FD分布 微观状态数 ΩMB=N!∏lal!∏lglal\Omega_{MB}=\frac{N!}{\prod_l a_l!}\prod_lg ...
- 跟我一起学jQuery——第二集（未完待续..）
<锋利的JQuery>第二版阅读笔记-第二章 跟我一起学jQuery--第二集 代码风格 jQuery选择器 接下来,就要开始正式学习jQuery的各种使用了.但是没有规矩不成方圆,所以我 ...
- CentOS7.3部署OpenStack-Ocata版本手记（计算节点） - 未完待续
一.环境准备 1. 网络环境 节点名称 IP 域名 控制节点 192.168.159.34 node1.example.local 计算节点 192.168.159.35 node2.example. ...
- 2020牛客暑期多校训练营（第二场）未完待续......
F. Fake Maxpooling 题目: 题目大意: 输入n,m,k.矩阵的尺寸为nm,其中每一个元素为A[i][j] = lcm( i , j ).从中找出所有kk的子矩阵中元素最大的数之和. ...
- 3D数学系列之——再谈蒙特卡洛积分和重要性采样
目录 一.前篇文章回顾 二.积分的黎曼和形式 三.积分的概率形式(蒙特卡洛积分) 四.误差 五.蒙特卡洛积分计算与收敛速度 六.重要性采样 七.重要性采样方法和过程 八.重要性采样的优缺点 一.前篇文 ...
- 3D数学系列之——从“蒙的挺准”到“蒙的真准”解密蒙特卡洛积分！
目录 1.前言 2.积分概念简单回顾 3.积分在程序计算上的困难 4.蒙特卡洛积分 5.一些扩展应用 1.前言 在学习3D数学的过程中,或者说在学习游戏开发.引擎开发.渲染器开发.Shader开发 ...
- AI学习者必备 | 圣母大学公开统计计算课程讲义（视频+PPT+作业）
翻译 | AI科技大本营(rgznai100) 参与 | 刘畅 近日,圣母大学(University of Notre Dame)公开了一门统计学课程资源,包括:课程笔记和授课视频,课后作业(以及解决 ...
- 概率论和蒙特卡洛积分
引用自: GAMES101-现代计算机图形学入门-闫令琪_哔哩哔哩_bilibili 概率论的介绍 离散随机变量 连续随机变量与概率密度函数(PDF) 因为是连续的,所以它的期望是用积分计算的 从变量 ...
最新文章
- C技巧：结构体参数转成不定参数
- 来人呐，有人又要抢钱啦！
- ccleaner无法更新_CCleaner正在静默更新关闭自动更新的用户
- 【JavaScript】apply和call的区别在哪？
- html页面内分栏显示不全,怎么消除Word文档分栏后栏间不平衡现象
- 计算机专业是理科吗,计算机类和普通理科有什么区别？
- 记号（notation）的学习
- java xsi type_java – JAXB – 如何根据XML值设置XML元素的xsi：type？
- 计算机材料学常用计算软件,计算机在材料科学中的应用-用MaterialsStudio计算简单材料的能带.doc...
- 推荐两款好用的视频压缩工具（在保证画质的情况下最大限度地压制）
- 《数据库原理》——知识点总结（期末复习）
- 王者战力查询教程，每天可查，数据准确~
- 【小猿說】以小刀会“的成败论当今创业成败
- You called this URL via POST, but the URL doesn‘t end in a slash and you hav
- 如何通过API方式集成金蝶ERP
- 南卡电容笔好还是ideo好？高性价比的电容笔测评
- 高清摄像机的发展历程及市场现状分析
- 计算机网络知识全面讲解：使用Telnet命令发送电子邮件
- 情理之中 - Macs do Windows
- 微信小程序-访问豆瓣电影api400错误