蒙特卡洛法求积分
问题一:我们如何用蒙特卡洛方法求积分?问题二:如何近似求一个随机变量的数学期望?问题三:估计的误差是多少?问题四:如何从理论上对蒙特卡洛估计做分析?结论
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
sns.set_style('whitegrid')
问题一:我们如何用蒙特卡洛方法求积分?
你眼中的蒙特卡洛方法求积分,可能是这样子的:
Image Name 最最经典的例子就是求 的近似值了,生成若干个均匀的点,然后统计在圆内的点的个数的比例,这个比例就是 的近似了! 但是这种方法的计算量非常大,而且随着维数的增长需要的计算量也会成倍上升,收敛速度也并不快。
但是我要讲的蒙特卡洛求积跟这个些许不一样。
假如我想求 ,我们来用概率的语言表达一下它。 设随机变量 ,即 上的均匀分布, 具有密度函数 。 那么就有: ,这个公式是下面推导中非常重要的一环。 事实上,借助这个公式,我们将求积分转化为求某个随机变量的数学期望!
接下来我们就可以用统计学的手段来处理了。
问题二:如何近似求一个随机变量的数学期望?
通常情况下, 都是一个比较复杂的函数。我们想要近似求期望,只能用统计学的手段。
设随机变量 ,一个常用的办法是,如果我们找到 个随机变量 的样本 那么 就是 一个好的近似!
容易知道,上式中的 服从 上的均匀分布。
所以我们的做法可以总结如下:
- 生成 个 上均匀分布的随机数 。
- 计算函数值
- 积分 具有近似值
接下来看一个例子来验证我们的结果:
, 用蒙特卡洛方法近似计算 ,我们知道真实值是 ,介于2.6到2.7之间 取
N = 1000
x = 2 * np.random.random(size=N) # [0, 2] 上的均匀分布
ans = 2 / N * np.sum(x ** 2)
print(ans)
2.6431046604313226
目前看来效果还算不错,但是问题就这样结束了么?
问题三:估计的误差是多少?
凡估计必有误差
每一次采样都可以得到一个估计值,我们多次采样,得到多个估计值,画出多个估计值的分布图,从图上就可以近似看出估计的误差了。
record = [] # 记录多次采样的估计值
N = 1000 # 每次采样取1000个点
for _ in range(1000):
x = 2 * np.random.random(size=N)
ans = 2 / N * np.sum(x ** 2)
record.append(ans)
plt.title('N=1000')
sns.distplot(record); plt.show()
惊讶的观察到:
- 虽然我们真实值是2.67左右,但是估计出2.5、2.9这种偏离程度大的值还是有可能的。
- 多次采样的结果分布像是正态分布。(这是巧合吗?)
把单次采样点增加到2000个,直觉告诉我们,误差会减小!
record = [] # 记录多次采样的估计值
N = 2000 # 每次采样取2000个点
for _ in range(1000):
x = 2 * np.random.random(size=N)
ans = 2 / N * np.sum(x ** 2)
record.append(ans)
plt.title('N=2000')
sns.distplot(record); plt.show()
明显看到分布更加紧凑了一些。
问题四:如何从理论上对蒙特卡洛估计做分析?
在统计学上,我们评价一个估计好不好的标准有哪些呢? 统计量有三大基本性质:
- 无偏性、有效性、相合性(一致性)
无偏性表示这个估计有没有 bias;有效性指这个估计的方差够不够小;相合性或者说一致性,说的是当样本容量非常大的时候,估计值是否趋近于真实值。
接下来我们从理论上来讨论这三点: 设 是 的一个估计量,则:
- 是无偏的,这点很好理解,对于样本 , 。
- 是相合估计量,根据大数定律,估计量 几乎绝对收敛与 。
- 有效性没法说明,但是我们可以知道 的方差为 (样本之间独立) 通过中心极限定理,还可以知道,当样本容量很大的时候, 渐进服从以 为均值, 为方差的正态分布。
最后,我想展示一下,本文所述的转化为估计随机变量期望的蒙特卡洛方法 与 传统的往正方形内投点计算落在圆内的点个数来估计 值的方法的不同。
Image Name
plt.figure(figsize=(12, 4))
record = []
for _ in range(1000):
x = np.random.random(size=(2, 200))
ans = 4 * np.mean(np.linalg.norm(x, axis=0) < 1) # 是否落在圆内
record.append(ans)
plt.subplot(121)
sns.distplot(record, kde=False)
record = []
for _ in range(1000):
x = np.random.random(size=200)
ans = 4 * np.mean(np.sqrt(1 - x ** 2)) # y = sqrt(1 - x ^ 2)
record.append(ans)
plt.subplot(122)
sns.distplot(record, kde=False); plt.show()
同样是取了2000个点(做200次计算),统计1000次结果。左图为传统方法,右图为本文所述转化为求期望的方法。
明显右边的效果更好!
结论
本文简要介绍了蒙特卡洛方法求积分的思路,以及相应的理论推导。蒙特卡洛求积分的本质是利用随机模拟估计一个随机变量的期望。理解好蒙特卡洛求积的思想有助于进一步学习MCMC方法。
进一步还可以思考:
- 如何用蒙特卡洛估计重积分?这种方法会随着维数的增大而出现计算困难吗?(即:是多项式级别的更困难还是指数级别的更困难)
- 说到用随机模拟估计一个随机变量的期望,如果我们能够获取到一定数量的某随机变量的样本,但是这个随机变量的分布非常畸形(考虑中国10亿劳动人口有6亿月收入不到1千,福布斯中国排行榜上前100的富豪能对中国人均收入产生巨大的影响),如何改进我们模拟的方法?
- 用贝叶斯判别分析方法预测股票涨跌
- 开发 | 在 Mac OS X 装不上 TensorFlow?看了这篇就会装
- 【答疑解惑】Java中的默认构造器和equals方法
- 原生JS | 当兔子遇到鸡
- 【Android基础】Activity的生命周期函数
- 七种常用回归技术,如何正确选择回归模型?
- 爬取拉勾网大数据相关岗位薪资信息存到excel,并作数据分析
- 【Windows编程】系列第五篇:GDI图形绘制
- 抓取链家官网北京房产信息并用python进行数据挖掘
- 用R语言做钻石价格预测
- 2分钟完成30*15页拉勾网职位需求关键词的抓取
- 【专业技术】Linux设备驱动第七篇:高级字符驱动操作之阻塞IO
- Python抓取上海各地区房价平均值
- R语言分析 老九门 到底谁是主角
- JavaScript 教程
- JavaScript 编辑工具
- JavaScript 与HTML
- JavaScript 与Java
- JavaScript 数据结构
- JavaScript 基本数据类型
- JavaScript 特殊数据类型
- JavaScript 运算符
- JavaScript typeof 运算符
- JavaScript 表达式
- JavaScript 类型转换
- JavaScript 基本语法
- JavaScript 注释
- Javascript 基本处理流程
- Javascript 选择结构
- Javascript if 语句
- Javascript if 语句的嵌套
- Javascript switch 语句
- Javascript 循环结构
- Javascript 循环结构实例
- Javascript 跳转语句
- Javascript 控制语句总结
- Javascript 函数介绍
- Javascript 函数的定义
- Javascript 函数调用
- Javascript 几种特殊的函数
- JavaScript 内置函数简介
- Javascript eval() 函数
- Javascript isFinite() 函数
- Javascript isNaN() 函数
- parseInt() 与 parseFloat()
- escape() 与 unescape()
- Javascript 字符串介绍
- Javascript length属性
- javascript 字符串函数
- Javascript 日期对象简介
- Javascript 日期对象用途
- Date 对象属性和方法
- Javascript 数组是什么
- Javascript 创建数组
- Javascript 数组赋值与取值
- Javascript 数组属性和方法
- 面向对象学习
- 常用模块
- 代理错误[WinError 10061]
- Linux系统JDK+Tomcat环境安装布署过程
- Python version 3.6 required, which was not found in the registry错误解决
- LNMP架构应用实战——Nginx服务介绍与安装
- 使用tidylib解决不规则网页问题
- LNMP架构应用实战——Nginx服务配置文件介绍
- Mac Sublime Text3快捷键
- Linux系统shell脚本编程——生产实战案例
- 学习python第一天总纲
- 学习python第二天数据库day1
- LNMP架构应用实战——Nginx配置虚拟主机
- 学习python第三天单行函数
- MySQL数据库入门——常用基础命令