大数定律与统计模拟（MCMC）

一、历史注记
    大数定律以数学形式表达了大量重复出现的随机现象的统计规律性，即频率的稳定性和平均结果的稳定性。该定律表述，在试验不变的条件下，重复试验多次，随机事件的频率近似于它的概率。偶然中包含着必然，必然的规律与特性在大量的样本中得以体现。

    最早的大数定律的表述可以追溯到公元1500年左右的意大利数学家Cardano；1713年，著名数学家James(Jacob) Bernouli正式提出并证明了最初的大数定律，不过当时现代概率论还没有建立起来，测度论、实分析的工具还没有出现，因此当时的大数定律是以``独立事件的概率''作为对象的；后来，历代数学家如Poisson（``大数定律''的名字来自于他）、Chebyshev、Markov、Khinchin（``强大数定律''的名字来自于他）、Borel、Cantelli等都对大数定律的发展做出了贡献；直到1930年，现代概率论奠基人、数学大师Kolgomorov才真正证明了最后的强大数定律。

    对于一般人来说，大数定律的非严格表述是这样的： \( X_{1},X_{2},\cdots,X_{n} \) 独立同分布，均值为 \( \mu \)，\( Y_{n}=X_{1}+X_{2}+\cdots+X_{n} \)，则 \( Y_{n}/n \) 收敛到 \( \mu \)。

    谈到``弱''与``强''的大数定律的区别，主要表现在结果上有差异。弱大数定律指 \( Y_{n}/n \) 依概率(in probability)收敛到 \( \mu \)，即对
\[\forall \varepsilon >0,\ \lim\limits_{n\rightarrow\infty}P\left\{|\frac{Y_{n}}{n}-\mu|\geq\varepsilon\right\}=0\]
强大数定律是指 \( Y_{n}/n \) 几乎处处收敛到 \( \mu \)，即 \( P\{\lim_{n\rightarrow\infty}Y_{n}/n=\mu\}=1 \) （注：几乎处处收敛意味着以概率收敛，《概率论教程》钟开莱 P71）。

假设 \( X_{1},X_{2},\cdots,X_{n} \) 是独立同分布随机变量序列，均值为 \( \mu \)。\( Y_{n}=\sum_{i=1}^n X_{i} \)，则弱大数定律成立的条件有以下几种：
(1) 方差 \( DX_{i} \) 一致有界；
证明方法：Chebyshev不等式，这个证明是Chebyshev给出的。（教学视频中的证明即是它，切比雪夫大数定律）。
(2) 均值 \( \mu \) 存在；
证明方法：用Taylor展开特征函数，证明其收敛到常数，得到依分布收敛，然后再用依分布收敛到常数等价于依概率收敛。

    强大数定律成立的条件有以下几种：
(3) 4阶矩存在；
(4) 2阶矩存在；
(5) 1阶矩存在。
从(1)到(5)，条件越来越弱，结果越来越强，证明也越来越难。

二、统计模拟------MCMC
    蒙特卡罗方法又称统计模拟法，也称为随机模拟法，它的理论基础就是大数定律。该方法的特点是：将所求解的问题与某概率模型相联系，用电子计算机实现统计抽样，按模拟方法获得问题的近似解。蒙特卡罗方法起源于于20世纪40年代，由美国在第二次世界大战中研制原子弹的``曼哈顿计划''的成员S.M.乌拉姆和J.冯诺依曼首先提出。
    数学家冯诺依曼用驰名世界的赌城------摩纳哥的Monte Carlo来命名这种方法，为它蒙上了一层神秘色彩。实际上，早在1777年蒙特卡罗方法就已经存在，当时法国数学家Buffon提出用投针实验的方法求圆周率 \( \pi \)。

1. MCMC求积分
例1  求曲线 \( \sqrt{x} \) 与直线 \( y=x \) 所围成的平面区域 \(A\) 的面积。
理论求解
\[ S_{A}=\int_0^1(\sqrt{x}-x)dx=\frac{1}{6} \]
MCMC 求解
    考虑两条直线 \( x=1,y=1 \) 和坐标轴围成边长为1的正方形；在正方形内随机投点，所投点的横、纵坐标 \( x \) 、\( y \) 均服从[0,1]上的均匀分布；所投点落到区域A内的概率 \( S_{A} \) 与正方形面积之比，即 \( S_{A}/1=S_{A} \)；当投点总数 \( n \) 足够大时，落到区域A的点的频率 \( m/n \) 近似表示概率，则有 \( S_{A}=m/n \)。

MCMC 实现的MATLAB程序：
function [S0,Sm]=quad1mont(n)
% [S0,Sm]=quad1mont(n),求曲线y=sqrt(x)与直线y=x所围成的区域A的面积
% 的理论值S0与蒙特卡洛模拟值Sm。输入参数n是随机投点的个数，可以是正整数
% 标量或向量
% S0=int('sqrt(x)-x',0,1);   %面积理论值（解析解）
S0=quad(@(x)sqrt(x)-x,0,1);   %面积理论值（数值解）
%计算区域A的面积的蒙特卡洛模拟值
for i=1:length(n)
    x=rand(n(i),1);    %点的横坐标
    y=rand(n(i),1);    %点的纵坐标
    m=sum(sqrt(x)>=y & y>=x);   %落到区域A内点的频数
    Sm(i)=m/n(i);
end

模拟结果：
>>[S0,Sm]=quad1mont ([100,1000,10000,100000,1000000])
S0 = 0.1667
Sm = 0.2200    0.1670    0.1712    0.1681    0.1666

例2 （二重积分）求球体 \( x^2+y^2+z^2=4 \) 被圆柱面 \( x^2+y^2=2x \) 所截得的（含在圆柱面内的部分）立体的体积。
                                               

理论求解:  记 \( D \) 为半圆周 \( y=\sqrt{2x-x^2} \) 及 \( x \) 轴所围成的闭区域（\(D\) 即是所截立体第一卦限部分在xoy面上的投影区域），则所求体积为
\[
V=4\iint\limits_{D}\sqrt{4-x^2-y^2}dxdy=4\iint\limits_{D}\sqrt{4-\rho^2}\rho d\rho d\theta=\frac{32}{3}(\frac{\pi}{2}-\frac{2}{3})\approx 9.644
\]
MCMC 求解：
    记 \( \Omega=\{ (x,y,z)|0\leq x\leq2,0\leq y\leq1,0\leq z\leq2\} \)，\( \Omega \) 是一个长方体区域；记所求立体在第一卦限部分为 \( T \)，\( T\subset\Omega \) 且 \( V=4V_{T} \)；在 \( \Omega \) 内随机投点，即所投点的坐标 \( X \) 、\( Y \) 、\( Z \) 分别服从[0,2]、[0,1]和[0,2]上的均匀分布；所投点落在 \( T \) 内的概率等于 \( T \) 的体积与 \( \Omega \) 的体积之比；当随机投点总数 \( n \) 足够大时，落在 \( T \) 内点的频率 \( m/n \) 近似于概率，即有 \( V_{T}/1=m/n \)，从而可得 \( V=16m/n \)。

MCMC 实现的MATLAB程序：
function [V0,Vm]=quad2mont(n)
% [V0,Vm]=quad2mont(n),求球面x^2+y^2+z^2=4被圆柱面x^2+y^2=2*x所
% 截得（含在圆柱面内的部分）立体的体积的理论值V0与蒙特卡洛模拟值Vm。
% 输入参数n是随机投点的个数，可以是正整数标量或向量
% V0=32*(pi/2-2/3)/3;   %面积理论值（解析解）
% 调用quad2d函数（MATLAB2009a中的新函数）求体积的理论值（数值解）
%V0=4*quad2d(@(x,y)sqrt(4-x.^2-y.^2),0,2,0,@(x)sqrt(1-(1-x).^2));
V0=4*quad(@(x)arrayfun(@(xx)quad(@(y)sqrt(4-xx.^2-y.^2),...
    0,sqrt(1-(1-xx).^2)),x),0,2);   %面积理论值（数值解）
%计算体积的蒙特卡洛模拟值
for i=1:length(n)
    x=2*rand(n(i),1);    %点的横坐标
    y=rand(n(i),1);      %点的纵坐标
    z=2*rand(n(i),1);    %点的竖坐标
m=sum((x.^2+y.^2+z.^2<=4) & ((x-1).^2+y.^2<=1));%落到区域T内点的频数
    Vm(i)=16*m/n(i);
end

模拟结果：
>> [v0,vm]=quad2mont([100,1000,10000,100000,1000000])
v0 = 9.6440
vm = 9.9200    9.7120    9.7312    9.6917    9.6444

2. 两个有趣问题的 MCMC

例1  三门问题（Monty Hall Problem）}
    假设你正在参加一个游戏节目，你被要求在三扇门中选择一扇：其中一扇后面有一辆车；其余两扇后面则是山羊。你选择了一道门，假设是一号门，然后确知门后面内容的主持人，开启了另一扇后面有山羊的门，假设是三号门。他然后问：你想选择二号门吗？

理论求解：引入两个随机事件： \( A= \) ``一开始选中汽车''，\( B= \) ``更换选择后选中汽车''，由全概率公式可求得
\[
P(B)=P(A)P(B|A)+P(\bar{A})P(B|\bar{A})=\frac{1}{3}\times 0+\frac{2}{3}\times1=\frac{2}{3}
\]
MCMC 求解：
    思路：设两只羊的编号分别为``1''和``2''，汽车编号为``3''；现从数字1、2、3中随机选取一个数字，若一开始选中1或2，则更换选择后选中``3''，即赢得汽车；若一开始选中``3''，则更换选择后选中1或2，即得不到汽车；将试验重复 \( n \) 次，记录一开始选中1或2的次数 \( m \)，从而可以确定更换选择后赢得汽车的频率为 \( m/n \)， \( n \) 足够大时，这个频率趋于更换选择后赢得汽车的概率。
    （请同学们自己编程试验模拟吧！）

例2  街头骗局揭秘
    街头常见一类"摸球游戏''，规则：一袋中装有16个大小、形状相同的玻璃球，其中8个红色、8个白色。游戏者从中一次摸出8个球，8个球中，当颜色出现以下比数时，摸球者可得到相应的``奖励''或惩罚，如下表
\[
\begin{array}{|c|c|c|c|c|c|}
\hline
\mbox{可能}&{A}&{B}&{C}&{D}&{E} \\
 \mbox{结果} &{8:0}&{7:1}&{6:2}&{5:3}&{4:4} \\\hline
\mbox{奖(罚)金}&{10}&{1}&{0.5}&{0.2}&{-3} \\
\hline
\end{array}
\]此游戏从表面上看非常有吸引力，5种可能出现的结果有4种可得到奖金，且最高奖达到10元，而只有1种情况受罚，罚金只有3元。试分析此游戏中谁是真正的赢家？

理论求解： 假设摸球者在一次游戏中得到的奖金（罚金）为 \( X \)，则
\[
P\{X=10\}=P(A)=2/C_{16}^{8}=0.0001554
\]\[
P\{X=1\}=P(B)=2C_{8}^{7}C_{8}^{1}/C_{16}^{8}=0.009946
\]\[
P\{X=0.5\}=P(C)=2C_{8}^{6}C_{8}^{2}/C_{16}^{8}=0.1218
\]\[
P\{X=0.2\}=P(D)=2C_{8}^{3}C_{8}^{5}/C_{16}^{8}=0.4873
\]\[
P\{X=-3\}=P(E)=C_{8}^{4}C_{8}^{4}/C_{16}^{8}=0.3807
\]则可得到 \( EX=-0.9723 \) 元，摸球者在一次游戏中平均要赔掉0.9723元。

MCMC 求解：
  思路：给16个球分别编号 \( 1\sim 16 \)，设8个红球编号为 \( 1\sim
8 \)。进行 \( n \) 次模拟，每次模拟生成8个随机整数（取值范围 \( 1\sim 16 \) ）作为一次抽取的8个球。统计 \( n \) 次模拟中各种可能的结果出现的频数 \( n_{A},n_{B},n_{C},n_{D},n_{E} \)，
从而可得摸球者在一次游戏中获得奖金的期望模拟值为：
\[
E_{m}=\frac{10n_{A}+n_{B}+0.5n_{C}+0.2n_{D}-3n_{E}}{n}
\]
    （请同学们自己编程试验模拟）