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

泊松分布与指数分布泊松分布是一个离散型分布其概率质量函数写作P(Zk)λke−λk!,k0,1,2,...P(Zk)\frac{\lambda^k e^{-\lambda}}{k!}, k0,1,2,...P(Zk)k!λke−λ​,k0,1,2,...其中λ0\lambda0λ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\lambda0λ0为参数其累积分布函数写作FZ(z∣λ)∫0∞λe−λzdz1−e−λzF_Z(z|\lambda) \int_0^\infty \lambda e^{-\lambda z}dz 1-e^{-\lambda z}FZ​(z∣λ)∫0∞​λe−λzdz1−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(Zk)(λt)ke−λtk!P(Zk)\frac{(\lambda t)^k e^{-\lambda t}}{k!}P(Zk)k!(λt)ke−λt​P(Z0)e−λtP(Z0) e^{-\lambda t}P(Z0)e−λt时间ttt内无事发生 等价于 第一次事件发生的时间TtTtTt因此有P(Tt)P(x0)e−λtP(Tt) P(x0) e^{-\lambda t}P(Tt)P(x0)e−λt时间ttt内至少有一件事发生的概率可写作P(T≤t)1−P(Tt)1−e−λtP(T\leq t) 1- P(Tt) 1 - e^{-\lambda t}P(T≤t)1−P(Tt)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_datanp.loadtxt(data.csv)# 先验定义似然中未知参数的分布alpha1/count_data.mean()withpm.Model()asmodel:# lambda服从指数分布tau服从离散均匀分布lambda_1pm.Exponential(lambda_1,alpha)lambda_2pm.Exponential(lambda_2,alpha)taupm.DiscreteUniform(tau,lower0,upperlen(count_data)-1)withmodel:#model没有被销毁# lambda是分段函数idxnp.arange(len(count_data))lambda_pm.math.switch(tauidx,lambda_1,lambda_2)# 似然给定参数得到数据的概率withmodel:observationpm.Poisson(obs,lambda_,observedcount_data)#每日短信数服从泊松分布# 后验结合数据计算似然中未知参数服从的分布withmodel:steppm.Metropolis()# 选择Metropolis采样算法tracepm.sample(10000,tune5000,stepstep,return_inferencedataFalse)# 获取先验withmodel:priorpm.sample_prior_predictive(samples10000)lambda_1_prior_samplesprior.prior[lambda_1].values.flatten()lambda_2_prior_samplesprior.prior[lambda_2].values.flatten()tau_prior_samplesprior.prior[tau].values.flatten()# 绘图显示plt.subplot(3,1,1)plt.hist(lambda_1_prior_samples,bins30)plt.subplot(3,1,2)plt.hist(lambda_2_prior_samples,bins30)plt.subplot(3,1,3)plt.hist(tau_prior_samples,bins30)plt.show()# 获取后验lambda_1_samplestrace[lambda_1]lambda_2_samplestrace[lambda_2]tau_samplestrace[tau]# 绘图显示plt.subplot(3,1,1)plt.hist(lambda_1_samples,bins30)plt.subplot(3,1,2)plt.hist(lambda_2_samples,bins30)plt.subplot(3,1,3)plt.hist(tau_samples,bins30)plt.show()