正文
技术背景
接上一篇文章，我们继续记录统计力学中的一些基础的概率论知识。这一篇文章主要介绍的是一些常用的概率密度函数的对应参数计算，如期望值、方差等。
伯努利分布
在离散分布中，最简单的分布为伯努利（Bernoulli）分布，也叫0-1分布。伯努利分布的随机变量就跟抛硬币一样只有两种：0（失败）和1（成功），对应的概率可以表示为：
\[p=P(X=1),q=P(X=0)=1-p
\]
因为伯努利分布是一个单次试验，因此根据期望值定义可以计算得：
\[\mu=E(X)=\sum_{i=0}^nX_iP(X_i)=X_0P(X=0)+X_1P(X=1)=0*q+1*p=p
\]
方差为：
\[E[(X-\mu)^2]=E(X^2-2\mu X+\mu^2)=E(X^2)-2\mu E(X)+\mu^2=E(X^2)-p^2\\
E(X^2)=\sum_{i=0}^nX_i^2P(X_i)=X_0^2P(X=0)+X_1^2P(X=1)=p\\
\sigma^2=E[(X-\mu)^2]=p-p^2=pq
\]
如果\(p=\frac{1}{2},q=1-p=\frac{1}{2}\)，那么对应的伯努利分布的期望值为\(\mu=p=\frac{1}{2}\)，方差为：\(\sigma^2=pq=\frac{1}{4}\)。
二项分布
二项分布是\(n\)次伯努利试验中成功次数的分布，记为\(B(n,p)\)，事件\(X_i\)表示在\(n\)次的伯努利试验中成功了\(i\)次。关于二项分布的期望值推导，可以使用期望值的定义：
\[E(X)=\sum_{i=0}^nX_iP(X_i)=\sum_{i=0}^nC_{n}^{i}X_ip^iq^{n-i}
\]
写到此处，我们应该要想起两数求和的\(n\)次方形式：
\[(p+q)^n=\sum_{i=0}^np^iq^{n-i}=1
\]
那么回头再计算期望值就可以直接使用这个结果：
\[\begin{align}
\mu=E(X)&=\sum_{i=0}^nC_{n}^{i}X_ip^iq^{n-i}\\
&=\sum_{i=1}^nC_n^iX_ip^iq^{n-i}+0\\
&=\sum_{i=1}^n\frac{n!}{i!(n-i)!}ip^iq^{(n-i)}\\
&=\sum_{i=1}^n\frac{n(n-1)!}{(i-1)![(n-1)-(i-1)]!}pp^{i-1}q^{[(n-1)-(i-1)]}\\
&=np\sum_{k=0}^{n-1}\frac{(n-1)!}{k!(n-1-k)!}p^kq^{(n-1-k)}\\
&=np(p+q)^{n-1}\\
&=np
\end{align}
\]
类似的，二项分布的方差为：
\[\begin{align}
\sigma^2=E[(X-\mu)^2]&=E(X^2)-\mu^2\\
&=\sum_{i=0}^nC_{n}^{i}X_i^2p^iq^{n-i}-\mu^2\\
&=\sum_{i=1}^n\frac{n!}{i!(n-i)!}i^2p^iq^{n-i}+0-\mu^2\\
&=\mu\sum_{k=0}^{n-1}\frac{(n-1)!}{k!(n-1-k)!}(k+1)p^kq^{(n-1-k)}-\mu^2\\
&=\mu\sum_{k=0}^m\frac{m!}{(k)!(m-k)!}kp^kq^{(m-k)}+\mu-\mu^2\\
&=\mu(n-1)p+\mu-\mu^2\\
&=np-np^2\\
&=npq
\end{align}
\]
关于这里面的概率系数\(C_{n}^{i}p^iq^{n-i}\)，我们可以把它画出来看一下：
import numpy as np
import matplotlib.pyplot as plt
def C(n,k,p):
return p**n if k==0 else (1-p)**n if k==n else p**k*(1-p)**(n-k)*np.prod(np.arange(1,n+1))/np.prod(np.arange(1,k+1))/np.prod(np.arange(1,n-k+1))
vC = lambda n,kk, p: [C(n,k,p) for k in kk]
N=20
p=0.5
kk = np.arange(N+1)
res = vC(N,kk,p)
plt.figure()
plt.plot(kk,res,color='black')
plt.plot(kk,res,'o',color='red')
plt.show()
泊松分布
泊松分布可以看作是二项分布的变种，事件都是一系列的伯努利分布，但是每个事件发生的概率由如下形式给出：
\[Pr(X=k)=p(k;\lambda)=\frac{e^{-\lambda}\lambda^k}{k!}
\]
同样的我们也可以把这个分布画出来：
import numpy as np
import matplotlib.pyplot as plt
N = 20
k = np.arange(N+1)
ps = lambda lbd, kk: [np.exp(-lbd) if k==0 else np.exp(-lbd)*(lbd**k)/np.prod(np.arange(1, k+1)) for k in kk]
plt.figure()
plt.plot(k, ps(0.1, k), label=r'$\lambda=0.1$')
plt.plot(k, ps(0.5, k), label=r'$\lambda=0.5$')
plt.plot(k, ps(1, k), label=r'$\lambda=1$')
plt.plot(k, ps(2, k), label=r'$\lambda=2$')
plt.plot(k, ps(3, k), label=r'$\lambda=3$')
plt.plot(k, ps(4, k), label=r'$\lambda=4$')
plt.plot(k, ps(5, k), label=r'$\lambda=5$')
plt.plot(k, ps(6, k), label=r'$\lambda=6$')
plt.legend()
plt.show()
得到的结果如下：
这是一个\(\lambda\)越小，整体采样分布越靠近\(0\)的概率函数。那么我们可以使用期望值的定义来计算泊松分布的期望值：
\[\mu=E(X)=\sum_{k=0}^NX_k\frac{e^{-\lambda}\lambda^k}{k!}=\lambda\sum_{k=1}^N \frac{e^{-\lambda}\lambda^{(k-1)}}{(k-1)!}=\lambda\sum_{n=0}^N \frac{e^{-\lambda}\lambda^n}{n!}=\lambda\sum_{n=0}^N Pr(X_n)=\lambda
\]
类似地，泊松分布的方差为：
\[\begin{align}
\sigma^2=E(X^2)-\mu^2&=\sum_{k=0}^Nk^2\frac{e^{-\lambda}\lambda^k}{k!}-\mu^2\\
&=\sum_{k=1}^Nk\frac{e^{-\lambda}\lambda^k}{(k-1)!}-\mu^2\\
&=e^{-\lambda}\lambda+\lambda\sum_{k=2}^Nk\frac{e^{-\lambda}\lambda^(k-1)}{(k-1)!}-\mu^2\\
&=e^{-\lambda}\lambda+\lambda\sum_{m=1}^N(m+1)\frac{e^{-\lambda}\lambda^m}{m!}-\mu^2\\
&=e^{-\lambda}\lambda+\lambda\sum_{m=1}^Nm\frac{e^{-\lambda}\lambda^m}{m!}+\lambda\sum_{m=1}^N\frac{e^{-\lambda}\lambda^m}{m!}-\mu^2\\
&=e^{-\lambda}\lambda+\lambda^2+\lambda(1-e^{-\lambda})-\lambda^2\\
&=\lambda
\end{align}
\]
正态分布
正态分布是一个极其重要的连续分布，其概率密度函数为：
\[f(x)=\frac{1}{\sqrt{2\pi}\sigma}e^{-\frac{1}{2}(\frac{x-\mu}{\sigma})^2}
\]
看到这个形式的概率密度函数，我们应该率先想到高斯积分：
\[I=\int_{-\infty}^{+\infty}e^{-x^2}dx
\]
对该积分取一个平方，形式变成了一个二元积分：
\[I^2=\int_{-\infty}^{+\infty}\int_{-\infty}^{+\infty}e^{-(x^2+y^2)}dxdy
\]
要求解这个积分，常规的方法是做一个极坐标变换，在知乎上有一个问题专门讨论了直角坐标系与极坐标系两个不同的矢量空间的积分问题。我们知道上述积分表达式中的\(dxdy\)表示的是一个面积元\(dS\)，那么能直接定义\(dS=dxdy\)吗？其实我们画两个图就看明白了：
在这个图中可以看出，其实只有在\(x,y\)为正交的坐标系下（即矢量\(d\textbf{x}\cdot d\textbf{y}\)=0时）才成立。关于面积元更加准确的定义应该是这样的：
\[d\textbf{S}=d\textbf{x}\times d\textbf{y}=dxdy\sin\theta
\]
因此对于一个坐标变换\(x=f(m,n),y=g(m,n)\)来说，对应的面积元的变换形式为：
\[\begin{align}
d\textbf{S}=d\textbf{x}\times d\textbf{y}&=\left[\frac{df(m,n)}{dm}d\textbf{m}+\frac{df(m,n)}{dn}d\textbf{n}\right]\times\left[\frac{dg(m,n)}{dm}d\textbf{m}+\frac{dg(m,n)}{dn}d\textbf{n}\right]\\
&=\frac{df(m,n)}{dm}\frac{dg(m,n)}{dn}d\textbf{m}\times d\textbf{n}+\frac{df(m,n)}{dn}\frac{dg(m,n)}{dm}d\textbf{n}\times d\textbf{m}\\
&=\left[\frac{df(m,n)}{dm}\frac{dg(m,n)}{dn}-\frac{df(m,n)}{dn}\frac{dg(m,n)}{dm}\right]d\textbf{m}\times d\textbf{n}
\end{align}
\]
此时代入一个极坐标变换\(x=f(r,\theta)=r\cos\theta,y=g(r,\theta)=r\sin\theta\)有：
\[\begin{align}
d\textbf{x}\times d\textbf{y}&=(\cos\theta\cdot r\cos\theta+r\sin\theta\cdot\sin\theta)d\vec{r}\times d\vec{\theta}\\
&=rd\vec{r}d\vec{\theta}
\end{align}
\]
所以常规的写法就是\(dxdy=rdrd\theta\)，那么此时再回到高斯积分：
\[\begin{align}
I^2&=\int_{-\pi}^{\pi}\int_{0}^{+\infty}e^{-r^2}rdrd\theta\\
&=-\frac{1}{2}\int_{-\pi}^{\pi}\left(\left.e^{-r^2}\right|_{0}^{+\infty}\right)d\theta\\
&=\frac{1}{2}\int_{-\pi}^{\pi}d\theta\\
&=\pi
\end{align}
\]
即\(I=\int_{-\infty}^{+\infty}e^{-x^2}dx=\sqrt{\pi}\)。此时回头看正态分布有：
\[\begin{align}
\int_{-\infty}^{+\infty}e^{-\frac{1}{2}(\frac{x-\mu}{\sigma})^2}dx&=\int_{-\infty}^{+\infty}e^{-(\frac{x}{\sqrt{2}\sigma})^2}dx\\
&=\sqrt{2}\sigma\int_{-\infty}^{+\infty}e^{-y^2}dy\\
&=\sqrt{2\pi}\sigma
\end{align}
\]
也就是说，\(\int_{-\infty}^{+\infty}\frac{1}{\sqrt{2\pi}\sigma}e^{-\frac{1}{2}(\frac{x-\mu}{\sigma})^2}dx\)，这表示正态分布给出的概率函数，在\(x\)轴上的积分为1，这符合我们对概率密度函数的预期。很显然的是，正态分布函数\(f(x)\)是关于\(x=\mu\)对称的一个函数，即\(P(x<\mu)=P(x>\mu)\)，按照连续分布的期望值定义：
\[\begin{align}
E(x)=\int xf(x)dx&=\int_{x<\mu} xf(x)dx+\int_{x=\mu}xf(x)dx+\int_{x>\mu}xf(x)dx\\
&=\int_{x<\mu} xf(x)dx+0+\int_{x<\mu} (2\mu-x)f(x)dx\\
&=2\mu\int_{x<\mu}f(x)dx\\
&=\mu
\end{align}
\]
至于方差，根据定义有：
\[\begin{align}
E[(x-\mu)^2]=\int_{-\infty}^{+\infty}(x-\mu)^2f(x)dx&=\int_{-\infty}^{+\infty}x^2f(x)dx-\mu^2\\
&=\int_{-\infty}^{+\infty}x^2\frac{1}{\sqrt{2\pi}\sigma}e^{-\frac{1}{2}(\frac{x-\mu}{\sigma})^2}dx-\mu^2
\end{align}
\]
考虑到积分项\(x\)带了一个平方，而概率幅在指数项上也带了一个\(x^2\)，因此我们可以考虑构造一个分部积分\((uv)'=u'v+uv'\rightarrow uv'=(uv)'-u'v\)来进行求解。首先我们定义一个变量\(y=\frac{x-\mu}{\sqrt{2}\sigma}\)便于求解，那么有：\(x=\sqrt{2}\sigma y+\mu, dx=\sqrt{2}\sigma dy\)，代入方差：
\[\begin{align}
E(x^2)&=\frac{1}{\sqrt{2\pi}\sigma}\int_{-\infty}^{+\infty}(\sqrt{2}\sigma y+\mu)^2e^{-y^2}\sqrt{2}\sigma dy\\
&=\frac{1}{\sqrt{\pi}}\int_{-\infty}^{+\infty}2\sigma^2y^2e^{-y^2}dy+\frac{1}{\sqrt{\pi}}\int_{-\infty}^{+\infty}2\sqrt{2}\mu\sigma ye^{-y^2}dy+\frac{1}{\sqrt{\pi}}\int_{-\infty}^{+\infty}\mu^2e^{-y^2}dy\\
&=\frac{2\sigma^2}{\sqrt{\pi}}\int_{-\infty}^{+\infty}y^2e^{-y^2}dy+\frac{2\sqrt{2}\mu\sigma}{\sqrt{pi}}\int_{-infty}^{+\infty}ye^{-y^2}dy+\frac{\mu^2}{\sqrt{\pi}}\int_{-\infty}^{+\infty}e^{-y^2}dy\\
&=\frac{2\sigma^2}{\sqrt{\pi}}\int_{-\infty}^{+\infty}y^2e^{-y^2}dy+\mu^2\\
&=\frac{\sigma^2}{\sqrt{\pi}}\int_{-\infty}^{+\infty}e^{-y^2}dy-\frac{\sigma^2}{\sqrt{\pi}}\left.\left(ye^{-y^2}\right)\right|_{-\infty}^{+\infty}+\mu^2\\
&=\sigma^2+\mu^2
\end{align}
\]
因此最终的方差为：
\[E[(x-\mu)^2]=E(x^2)-\mu^2=\sigma^2
\]
这里面用到了一个洛必达法则：\(\lim_{x\to -\infty}xe^{-x^2}=\lim_{x\to +\infty}xe^{-x^2}=\lim_{x\to \infty}\frac{1}{2xe^{x^2}}=0\)。因此，正态分布函数中的\(\mu\)和\(\sigma\)也是它的期望值和标准差。特别地，当\(\mu=0,\sigma=1\)时，该分布称为标准正态分布。
均匀分布
均匀分布，其实可以认为是所有分布的基础。我们可以通过均匀随机数结合一个累积分布函数，去构造其概率密度采样。均匀随机数的概率密度函数为：
\[f(x)=\left\{
\begin{matrix}
\frac{1}{b-a}, && a\leq x\leq b\\
0, && others
\end{matrix}
\right.
, a<b
\]
那么这个函数在整个轴上的积分为：\(\int_{-\infty}^{+\infty}\frac{1}{b-a}dx=\frac{1}{b-a}\left.(x)\right|_{a}^b=1\)。类似的，也可以简单的推导均匀分布的期望值：
\[\mu=\frac{1}{b-a}\int_{a}^bxdx=\frac{1}{2(b-a)}\left.(x^2)\right|_{a}^b=\frac{a+b}{2}
\]
还有均匀分布的方差：
\[\sigma^2=\frac{1}{b-a}\int_{a}^bx^2dx=\frac{1}{3(b-a)}\left.(x^3)\right|_{a}^b=\frac{(b-a)^2}{12}
\]
其实关于均匀分布，更多的时候选用的是\(x~U(0,1)\)，也就是0~1区间内的均匀随机采样，这也是各种编程框架下默认使用的均匀分布随机数。其实别看这个均匀随机数形式上简简单单，在蒙特卡洛采样中有非常多的应用，比如可以用均匀采样估计圆周率\(\pi\)的值：
import numpy as np
res = lambda n: np.sum(np.linalg.norm(np.random.random(size=(n,2)),axis=-1)<1)*4./n
for i in [10,100,1000,10000]: print ('The {} samples evaluation value of π is: {}'.format(i,res(i)))
# The 10 samples evaluation value of π is: 3.2
# The 100 samples evaluation value of π is: 3.44
# The 1000 samples evaluation value of π is: 3.284
# The 10000 samples evaluation value of π is: 3.1808
这个算法就是在一个边长为1的正方形内均匀采样，然后取\(\frac{1}{4}\)圆内的点做统计，最终落在圆形内部的点会占总样本数的\(\frac{\pi}{4}\)，以此来估计圆周率\(\pi\)。
指数分布
指数分布是一个“单调”分布，其概率密度函数形式为：
\[f(x)=\lambda e^{-\lambda x}, x\geq 0
\]
首先还是对这个概率密度函数做一个归一化的校验：\(\int_{0}^{+\infty}\lambda e^{-\lambda x}dx=-\left.e^{-\lambda x}\right|_{0}^{+\infty}=1\)，然后照例计算其期望值：
\[\mu=\lambda\int_{0}^{+\infty}xe^{-\lambda x}dx=\int_{0}^{+\infty}e^{-\lambda x}dx-\left.(xe^{-\lambda x})\right|_{0}^{+\infty}=\frac{1}{\lambda}
\]
指数分布的期望值为：
\[\sigma^2=\lambda\int_{0}^{+\infty}x^2e^{-\lambda x}dx-\mu^2=\int_0^{+\infty}2xe^{-\lambda x}dx-\left.(x^2e^{-\lambda x})\right|_0^{+\infty}-\mu^2=\frac{2}{\lambda^2}-\frac{1}{\lambda^2}=\frac{1}{\lambda^2}
\]
类似于正态分布中的计算，这里也用到了分部积分和洛必达法则。
大数定理
简单理解大数定理就是，在执行的采样数量\(n\)足够大的时候，样本概率\(\frac{m}{n}\)会趋近于真实概率\(p\)，也叫依概率收敛：
\[\lim_{n\rightarrow \infty}P(\left|\frac{m}{n}-p\right|<\epsilon)=1
\]
这里\(\epsilon\rightarrow0\)是一个小数。除了事件概率收敛外，还有期望值大数定理和方差大数定理，都是依概率收敛：
\[\lim_{n\rightarrow \infty}P(\left|\frac{1}{n}\sum_{i=1}^nX_i-\mu\right|<\epsilon)=1,\lim_{n\rightarrow \infty}P(\left|\frac{1}{n}\sum_{i=1}^n\left(X_i-\mu\right)^2-\sigma^2\right|<\epsilon)=1
\]
总结概要
可以理解的是，概率密度函数，一般情况下都是连续的。但是对于采样或者随机试验来说，其实都是离散采样。大数定理通过取一个极限，将概率密度函数跟试验联系了起来。这篇文章主要介绍的是常用的几个概率密度函数的期望值和方差的计算，以及大数定理的基本概念。
版权声明
本文首发链接为：https://www.cnblogs.com/dechinphy/p/prob-2.html
作者ID：DechinPhy
更多原著文章：https://www.cnblogs.com/dechinphy/
请博主喝咖啡：https://www.cnblogs.com/dechinphy/gallery/image/379634.html
参考资料
- 《统计力学导引》--郑伟谋
