泊松分布与指数分布以及一道贝叶斯推断例题

泊松分布与指数分布

泊松分布是一个离散型分布，其概率质量函数写作：
P(Z=k)=λke−λk!,k=0,1,2,...P(Z=k)=\frac{\lambda^k e^{-\lambda}}{k!}, k=0,1,2,...P(Z=k)=k!λke−λ​,k=0,1,2,...
其中，λ>0\lambda>0λ>0为参数，kkk为非负整数

如果Z∼Poi(λ)Z \sim Poi(\lambda)Z∼Poi(λ)，则有E[Z∣λ]=λ,Var[Z∣λ]=λE[Z|\lambda] = \lambda, Var[Z|\lambda] = \lambdaE[Z∣λ]=λ,Var[Z∣λ]=λ

指数分布是一个连续型分布，其概率密度函数写作：
fZ(z∣λ)=λe−λz,z≥0f_Z(z|\lambda) = \lambda e^{-\lambda z}, z\geq 0fZ​(z∣λ)=λe−λz,z≥0
其中，λ>0\lambda>0λ>0为参数
其累积分布函数写作：
FZ(z∣λ)=∫0∞λe−λzdz=1−e−λzF_Z(z|\lambda) = \int_0^\infty \lambda e^{-\lambda z}dz = 1-e^{-\lambda z}FZ​(z∣λ)=∫0∞​λe−λzdz=1−e−λz

如果Z∼Exp(λ)Z \sim Exp(\lambda)Z∼Exp(λ)，则有E[Z∣λ]=1/λ,Var[Z∣λ]=1/λ2E[Z|\lambda] = 1/\lambda, Var[Z|\lambda] = 1/\lambda^2E[Z∣λ]=1/λ,Var[Z∣λ]=1/λ2

λ\lambdaλ是否固定可以区分频率派与贝叶斯派，后者认为分布的参数是随机变量。

泊松分布与指数分布的关系

泊松分布描述单位时间内事件发生的次数，指数分布描述事件发生的时间间隔

令单位时间内事件发生的次数为λ\lambdaλ，时间ttt内共发生λt\lambda tλt次
时间ttt内的泊松分布如下：
P(Z=k)=(λt)ke−λtk!P(Z=k)=\frac{(\lambda t)^k e^{-\lambda t}}{k!}P(Z=k)=k!(λt)ke−λt​
P(Z=0)=e−λtP(Z=0)= e^{-\lambda t}P(Z=0)=e−λt

时间ttt内无事发生 等价于 第一次事件发生的时间T>tT>tT>t，因此有：
P(T>t)=P(x=0)=e−λtP(T>t) = P(x=0) = e^{-\lambda t}P(T>t)=P(x=0)=e−λt
时间ttt内至少有一件事发生的概率可写作：
P(T≤t)=1−P(T>t)=1−e−λtP(T\leq t) = 1- P(T>t) = 1 - e^{-\lambda t}P(T≤t)=1−P(T>t)=1−e−λt
这就是指数分布的累计函数FT(t)F_T(t)FT​(t)

The Poisson and Exponential Distribution, Jhon.C.B.Cooper. https://neurophysics.ucsd.edu/courses/physics_171/exponential.pdf

例题

例题来源于贝叶斯方法 概率编程与贝叶斯推断[M]
数据与代码可以在链接中找到https://github.com/CamDavidsonPilon/Probabilistic-Programming-and-Bayesian-Methods-for-Hackers/blob/master/Chapter1_Introduction/Ch1_Introduction_PyMC_current.ipynb

现在有某用户70天的短信接收量数据，对其建模并根据数据估计参数的分布？

注：我只是直观感受一下贝叶斯方法如何解决问题，先不作考虑为什么λ\lambdaλ服从指数分布以及为什么设计两个λ\lambdaλ但又令α\alphaα等于短信数的均值分之一。

建模：
1、每天接收的短信数CiC_iCi​服从泊松分布，即Ci∼Poi(λ)C_i \sim Poi(\lambda)Ci​∼Poi(λ)
2、λ\lambdaλ是一个分段函数，λ=λ1,t<τ;λ2,t≥τ\lambda = \lambda_1, t<\tau; \lambda_2, t\geq \tauλ=λ1​,t<τ;λ2​,t≥τ
3、λ∼Exp(α)\lambda \sim Exp(\alpha)λ∼Exp(α)，令α=1/ΣiCi70\alpha = 1/\Sigma_i \frac{C_i}{70}α=1/Σi​70Ci​​
4、τ∼DiscreteUniform(1,70)\tau \sim DiscreteUniform(1,70)τ∼DiscreteUniform(1,70)

importnumpyasnpimportmatplotlib.pyplotaspltimportpymcaspmif__name__=='__main__':# 加载数据count_data=np.loadtxt("data.csv")# 先验：定义似然中未知参数的分布alpha=1/count_data.mean()withpm.Model()asmodel:# lambda服从指数分布，tau服从离散均匀分布lambda_1=pm.Exponential("lambda_1",alpha)lambda_2=pm.Exponential("lambda_2",alpha)tau=pm.DiscreteUniform("tau",lower=0,upper=len(count_data)-1)withmodel:#model没有被销毁# lambda是分段函数idx=np.arange(len(count_data))lambda_=pm.math.switch(tau>idx,lambda_1,lambda_2)# 似然：给定参数得到数据的概率withmodel:observation=pm.Poisson("obs",lambda_,observed=count_data)#每日短信数服从泊松分布# 后验：结合数据计算似然中未知参数服从的分布withmodel:step=pm.Metropolis()# 选择Metropolis采样算法trace=pm.sample(10000,tune=5000,step=step,return_inferencedata=False)# 获取先验withmodel:prior=pm.sample_prior_predictive(samples=10000)lambda_1_prior_samples=prior.prior["lambda_1"].values.flatten()lambda_2_prior_samples=prior.prior["lambda_2"].values.flatten()tau_prior_samples=prior.prior["tau"].values.flatten()# 绘图显示plt.subplot(3,1,1)plt.hist(lambda_1_prior_samples,bins=30)plt.subplot(3,1,2)plt.hist(lambda_2_prior_samples,bins=30)plt.subplot(3,1,3)plt.hist(tau_prior_samples,bins=30)plt.show()# 获取后验lambda_1_samples=trace['lambda_1']lambda_2_samples=trace['lambda_2']tau_samples=trace['tau']# 绘图显示plt.subplot(3,1,1)plt.hist(lambda_1_samples,bins=30)plt.subplot(3,1,2)plt.hist(lambda_2_samples,bins=30)plt.subplot(3,1,3)plt.hist(tau_samples,bins=30)plt.show()