统计力学中的概率论基础(二)

技术背景

上一篇文章,我们继续记录统计力学中的一些基础的概率论知识。这一篇文章主要介绍的是一些常用的概率密度函数的对应参数计算,如期望值、方差等。

伯努利分布

在离散分布中,最简单的分布为伯努利(Bernoulli)分布,也叫0-1分布。伯努利分布的随机变量就跟抛硬币一样只有两种:0(失败)和1(成功),对应的概率可以表示为:

p=P(X=1),q=P(X=0)=1p

因为伯努利分布是一个单次试验,因此根据期望值定义可以计算得:

μ=E(X)=i=0nXiP(Xi)=X0P(X=0)+X1P(X=1)=0q+1p=p

方差为:

E[(Xμ)2]=E(X22μX+μ2)=E(X2)2μE(X)+μ2=E(X2)p2E(X2)=i=0nXi2P(Xi)=X02P(X=0)+X12P(X=1)=pσ2=E[(Xμ)2]=pp2=pq

如果p=12,q=1p=12,那么对应的伯努利分布的期望值为μ=p=12,方差为:σ2=pq=14

二项分布

二项分布n次伯努利试验中成功次数的分布,记为B(n,p),事件Xi表示在n次的伯努利试验中成功了i次。关于二项分布的期望值推导,可以使用期望值的定义:

E(X)=i=0nXiP(Xi)=i=0nCniXipiqni

写到此处,我们应该要想起两数求和的n次方形式:

(p+q)n=i=0npiqni=1

那么回头再计算期望值就可以直接使用这个结果:

(1)μ=E(X)=i=0nCniXipiqni(2)=i=1nCniXipiqni+0(3)=i=1nn!i!(ni)!ipiq(ni)(4)=i=1nn(n1)!(i1)![(n1)(i1)]!ppi1q[(n1)(i1)](5)=npk=0n1(n1)!k!(n1k)!pkq(n1k)(6)=np(p+q)n1(7)=np

类似的,二项分布的方差为:

(8)σ2=E[(Xμ)2]=E(X2)μ2(9)=i=0nCniXi2piqniμ2(10)=i=1nn!i!(ni)!i2piqni+0μ2(11)=μk=0n1(n1)!k!(n1k)!(k+1)pkq(n1k)μ2(12)=μk=0mm!(k)!(mk)!kpkq(mk)+μμ2(13)=μ(n1)p+μμ2(14)=npnp2(15)=npq

关于这里面的概率系数Cnipiqni,我们可以把它画出来看一下:

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;λ)=eλλkk!

同样的我们也可以把这个分布画出来:

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()

得到的结果如下:

这是一个λ越小,整体采样分布越靠近0的概率函数。那么我们可以使用期望值的定义来计算泊松分布的期望值:

μ=E(X)=k=0NXkeλλkk!=λk=1Neλλ(k1)(k1)!=λn=0Neλλnn!=λn=0NPr(Xn)=λ

类似地,泊松分布的方差为:

(16)σ2=E(X2)μ2=k=0Nk2eλλkk!μ2(17)=k=1Nkeλλk(k1)!μ2(18)=eλλ+λk=2Nkeλλ(k1)(k1)!μ2(19)=eλλ+λm=1N(m+1)eλλmm!μ2(20)=eλλ+λm=1Nmeλλmm!+λm=1Neλλmm!μ2(21)=eλλ+λ2+λ(1eλ)λ2(22)=λ

正态分布

正态分布是一个极其重要的连续分布,其概率密度函数为:

f(x)=12πσe12(xμσ)2

看到这个形式的概率密度函数,我们应该率先想到高斯积分:

I=+ex2dx

对该积分取一个平方,形式变成了一个二元积分:

I2=++e(x2+y2)dxdy

要求解这个积分,常规的方法是做一个极坐标变换,在知乎上有一个问题专门讨论了直角坐标系与极坐标系两个不同的矢量空间的积分问题。我们知道上述积分表达式中的dxdy表示的是一个面积元dS,那么能直接定义dS=dxdy吗?其实我们画两个图就看明白了:

在这个图中可以看出,其实只有在x,y为正交的坐标系下(即矢量dxdy=0时)才成立。关于面积元更加准确的定义应该是这样的:

dS=dx×dy=dxdysinθ

因此对于一个坐标变换x=f(m,n),y=g(m,n)来说,对应的面积元的变换形式为:

(23)dS=dx×dy=[df(m,n)dmdm+df(m,n)dndn]×[dg(m,n)dmdm+dg(m,n)dndn](24)=df(m,n)dmdg(m,n)dndm×dn+df(m,n)dndg(m,n)dmdn×dm(25)=[df(m,n)dmdg(m,n)dndf(m,n)dndg(m,n)dm]dm×dn

此时代入一个极坐标变换x=f(r,θ)=rcosθ,y=g(r,θ)=rsinθ有:

(26)dx×dy=(cosθrcosθ+rsinθsinθ)dr×dθ(27)=rdrdθ

所以常规的写法就是dxdy=rdrdθ,那么此时再回到高斯积分:

(28)I2=ππ0+er2rdrdθ(29)=12ππ(er2|0+)dθ(30)=12ππdθ(31)=π

I=+ex2dx=π。此时回头看正态分布有:

(32)+e12(xμσ)2dx=+e(x2σ)2dx(33)=2σ+ey2dy(34)=2πσ

也就是说,+12πσe12(xμσ)2dx,这表示正态分布给出的概率函数,在x轴上的积分为1,这符合我们对概率密度函数的预期。很显然的是,正态分布函数f(x)是关于x=μ对称的一个函数,即P(x<μ)=P(x>μ),按照连续分布的期望值定义:

(35)E(x)=xf(x)dx=x<μxf(x)dx+x=μxf(x)dx+x>μxf(x)dx(36)=x<μxf(x)dx+0+x<μ(2μx)f(x)dx(37)=2μx<μf(x)dx(38)=μ

至于方差,根据定义有:

(39)E[(xμ)2]=+(xμ)2f(x)dx=+x2f(x)dxμ2(40)=+x212πσe12(xμσ)2dxμ2

考虑到积分项x带了一个平方,而概率幅在指数项上也带了一个x2,因此我们可以考虑构造一个分部积分(uv)=uv+uvuv=(uv)uv来进行求解。首先我们定义一个变量y=xμ2σ便于求解,那么有:x=2σy+μ,dx=2σdy,代入方差:

(41)E(x2)=12πσ+(2σy+μ)2ey22σdy(42)=1π+2σ2y2ey2dy+1π+22μσyey2dy+1π+μ2ey2dy(43)=2σ2π+y2ey2dy+22μσpiinfty+yey2dy+μ2π+ey2dy(44)=2σ2π+y2ey2dy+μ2(45)=σ2π+ey2dyσ2π(yey2)|++μ2(46)=σ2+μ2

因此最终的方差为:

E[(xμ)2]=E(x2)μ2=σ2

这里面用到了一个洛必达法则:limxxex2=limx+xex2=limx12xex2=0。因此,正态分布函数中的μσ也是它的期望值和标准差。特别地,当μ=0,σ=1时,该分布称为标准正态分布。

均匀分布

均匀分布,其实可以认为是所有分布的基础。我们可以通过均匀随机数结合一个累积分布函数,去构造其概率密度采样。均匀随机数的概率密度函数为:

f(x)={1ba,axb0,others,a<b

那么这个函数在整个轴上的积分为:+1badx=1ba(x)|ab=1。类似的,也可以简单的推导均匀分布的期望值:

μ=1baabxdx=12(ba)(x2)|ab=a+b2

还有均匀分布的方差:

σ2=1baabx2dx=13(ba)(x3)|ab=(ba)212

其实关于均匀分布,更多的时候选用的是x U(0,1),也就是0~1区间内的均匀随机采样,这也是各种编程框架下默认使用的均匀分布随机数。其实别看这个均匀随机数形式上简简单单,在蒙特卡洛采样中有非常多的应用,比如可以用均匀采样估计圆周率π的值:

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的正方形内均匀采样,然后取14圆内的点做统计,最终落在圆形内部的点会占总样本数的π4,以此来估计圆周率π

指数分布

指数分布是一个“单调”分布,其概率密度函数形式为:

f(x)=λeλx,x0

首先还是对这个概率密度函数做一个归一化的校验:0+λeλxdx=eλx|0+=1,然后照例计算其期望值:

μ=λ0+xeλxdx=0+eλxdx(xeλx)|0+=1λ

指数分布的期望值为:

σ2=λ0+x2eλxdxμ2=0+2xeλxdx(x2eλx)|0+μ2=2λ21λ2=1λ2

类似于正态分布中的计算,这里也用到了分部积分和洛必达法则。

大数定理

简单理解大数定理就是,在执行的采样数量n足够大的时候,样本概率mn会趋近于真实概率p,也叫依概率收敛:

limnP(|mnp|<ϵ)=1

这里ϵ0是一个小数。除了事件概率收敛外,还有期望值大数定理和方差大数定理,都是依概率收敛:

limnP(|1ni=1nXiμ|<ϵ)=1,limnP(|1ni=1n(Xiμ)2σ2|<ϵ)=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

参考资料

  1. 《统计力学导引》--郑伟谋

本文作者:Dechin的博客

本文链接:https://www.cnblogs.com/dechinphy/p/18191203/prob-2

版权声明:本作品采用CC BY-NC-SA 4.0许可协议进行许可。

posted @   DECHIN  阅读(238)  评论(0编辑  收藏  举报
点击右上角即可分享
微信分享提示
评论
收藏
关注
推荐
深色
回顶
收起

喜欢请打赏

扫描二维码打赏

了解更多