如何用 python 对生存分析数据样本进行 Weibull 分布参数的最大似然估计 python 参数估计
在这一节,我们将讨论贝叶斯方法是如何思考数据的,我们怎样通过 MCMC 技术估计模型参数。
from IPython.display import Image
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import pymc3 as pm
import scipy
import scipy.stats as stats
import scipy.optimize as opt
import statsmodels.api as sm
%matplotlib inline
plt.style.use('bmh')
colors = ['#348ABD', '#A60628', '#7A68A6', '#467821', '#D55E00',
'#CC79A7', '#56B4E9', '#009E73', '#F0E442', '#0072B2']
messages = pd.read_csv('data/hangout_chat_data.csv')
贝叶斯方法如何思考数据?
当我开始学习如何运用贝叶斯方法的时候,我发现理解贝叶斯方法如何思考数据是很有用的。设想下述的场景:
一个好奇的男孩每天观察从他家门前经过的汽车的数量。他每天努力地记录汽车的总数。一个星期过去后,他的笔记本上记录着下面的数字:12,33,20,29,20,30,18
从贝叶斯方法的角度看,这个数据是由随机过程产生的。但是,既然数据被观测,它便固定了并且不会改变。这个随机过程有些模型参数被固定了。然而,贝叶斯方法用概率分布来表示这些模型参数的不确定性。
由于这个男孩调查的是计数(非负整数),一个通常的做法是用泊松分布对数据(如随机过程)建模。泊松分布只有一个参数 μ,它既是数据的平均数,也是方差。下面是三个不同 μ 值的泊松分布。
代码:
fig = plt.figure(figsize=(11,3))
ax = fig.add_subplot(111)
x_lim = 60
mu = [5, 20, 40]
for i in np.arange(x_lim):
plt.bar(i, stats.poisson.pmf(mu[0], i), color=colors[3])
plt.bar(i, stats.poisson.pmf(mu[1], i), color=colors[4])
plt.bar(i, stats.poisson.pmf(mu[2], i), color=colors[5])
_ = ax.set_xlim(0, x_lim)
_ = ax.set_ylim(0, 0.2)
_ = ax.set_ylabel('Probability mass')
_ = ax.set_title('Poisson distribution')
_ = plt.legend(['$mu$ = %s' % mu[0], '$mu$ = %s' % mu[1], '$mu$ = %s' % mu[2]])
在上一节中,我们引入我的 hangout 聊天数据集。特别地,我对我的回复时间(response_time
)感兴趣。鉴于 response_time
是计数数据,我们可以用泊松分布对其建模,并估计参数 μ 。我们将用频率论方法和贝叶斯方法两种方法来估计。
fig = plt.figure(figsize=(11,3))
_ = plt.title('Frequency of messages by response time')
_ = plt.xlabel('Response time (seconds)')
_ = plt.ylabel('Number of messages')
_ = plt.hist(messages['time_delay_seconds'].values,
range=[0, 60], bins=60, histtype='stepfilled')
频率论方法估计μ
在进入贝叶斯方法之前,让我们先看一下频率论方法估计泊松分布参数。我们将使用优化技术使似然函数最大。
下面的函数poisson_logprob()返回在给定泊松分布模型和参数值的条件下,观测数据总体的可能性。用方法opt.minimize_scalar找到在观测数据基础上参数值 μ 的最可信值(最大似然)。该方法的机理是,这个优化技术会自动迭代可能的mu值直到找到可能性最大的值。
y_obs = messages['time_delay_seconds'].values
def poisson_logprob(mu, sign=-1):
return np.sum(sign*stats.poisson.logpmf(y_obs, mu=mu))
freq_results = opt.minimize_scalar(poisson_logprob)
%time print("The estimated value of mu is: %s" % freq_results['x'])
所以,μ 的估计值是18.0413533867。优化技术没有对不确定度进行评估,它只返回一个点,效率很高。
下图描述的是我们优化的函数。对于每个μ值,图线显示给定数据和模型在μ处的似然度。优化器以登山模式工作——从曲线上随机一点开始,不停向上攀登直到达到最高点。
x = np.linspace(1, 60)
y_min = np.min([poisson_logprob(i, sign=1) for i in x])
y_max = np.max([poisson_logprob(i, sign=1) for i in x])
fig = plt.figure(figsize=(6,4))
_ = plt.plot(x, [poisson_logprob(i, sign=1) for i in x])
_ = plt.fill_between(x, [poisson_logprob(i, sign=1) for i in x],
y_min, color=colors[0], alpha=0.3)
_ = plt.title('Optimization of $mu$')
_ = plt.xlabel('$mu$')
_ = plt.ylabel('Log probability of $mu$ given data')
_ = plt.vlines(freq_results['x'], y_max, y_min, colors='red', linestyles='dashed')
_ = plt.scatter(freq_results['x'], y_max, s=110, c='red', zorder=3)
_ = plt.ylim(ymin=y_min, ymax=0)
_ = plt.xlim(xmin=1, xmax=60)
上述优化方法估计泊松分布的参数值为 18 .我们知道,对任意一个泊松分布,参数 μ 代表的是均值和方差。下图描述了这个分布。
fig = plt.figure(figsize=(11,3))
ax = fig.add_subplot(111)
x_lim = 60
mu = np.int(freq_results['x'])
for i in np.arange(x_lim):
plt.bar(i, stats.poisson.pmf(mu, i), color=colors[3])
_ = ax.set_xlim(0, x_lim)
_ = ax.set_ylim(0, 0.1)
_ = ax.set_xlabel('Response time in seconds')
_ = ax.set_ylabel('Probability mass')
_ = ax.set_title('Estimated Poisson distribution for Hangout chat response time')
_ = plt.legend(['$lambda$ = %s' % mu])
上述泊松分布模型和 μ 的估计表明,观测小于10 或大于 30 的可能性很小,绝大多数的概率分布在 10 和 30 之间。但是,这不能反映我们观测到的介于 0 到 60 之间的数据。
贝叶斯方法估计 μ
如果你之前遇到过贝叶斯理论,下面的公式会看起来很熟悉。我从没认同过这个框架直到我读了John K. Kruschke的书“贝叶斯数据分析”,从他优美简洁的贝叶斯图表模型视角认识这个公式。
代码:
Image('graphics/Poisson-dag.png', width=320)
上述模式可以解读如下(从下到上):
- 对每一个对话(i)观测计数数据(y)
- 数据由一个随机过程产生,我们认为可以用泊松分布模拟
- 泊松分布只有一个参数,介于 0 到 60 之间
由于我们没有任何依据对μ在这个范围内进行预估,就用均匀分布对 μ 建模
神奇的方法MCMC
下面的动画很好的演示了马尔科夫链蒙特卡罗方法(MCMC)的过程。MCMC采样器从先验分布中选取参数值,计算从这些参数值决定的某个分布中得到观测数据的可能性。
这个计算可以作为MCMC采样器的导引。从参数的先验分布中选取值,然后计算给定数据条件下这些参数值可能性——导引采样器到更高概率的区域。
与上面讨论的频率论优化技术在概念上有相似之处,MCMC采样器向可能性最高的区域运动。但是,贝叶斯方法不关心寻找绝对最大值,而是遍历和收集概率最高区域附近的样本。所有收集到的样本都可认为是一个可信的参数。
Image(url='graphics/mcmc-animate.gif')
with pm.Model() as model:
mu = pm.Uniform('mu', lower=0, upper=60)
likelihood = pm.Poisson('likelihood', mu=mu, observed=messages['time_delay_seconds'].values)
start = pm.find_MAP()
step = pm.Metropolis()
trace = pm.sample(200000, step, start=start, progressbar=True)
Applied interval-transform to mu and added transformed mu_interval to model.
[-----------------100%-----------------] 200000 of 200000 complete in 48.3 sec
上面的代码通过遍历 μ 的后验分布的高概率区域,收集了 200,000 个 μ 的可信样本。下图(左)显示这些值的分布,平均值与频率论的估值(红线)几乎一样。但是,我们同时也得到了不确定度,μ 值介于 17 到 19 之间都是可信的。我们稍后会看到这个不确定度很有价值。
_ = pm.traceplot(trace, varnames=['mu'], lines={'mu': freq_results['x']})
丢弃早期样本(坏点)
你可能疑惑上面 MCMC 代码中pm.find.MAP()
的作用。MAP 代表最大后验估计,它帮助 MCMC 采样器寻找合适的采样起始点。理想情况下,模型从高概率区域开始,但有时不是这样。因此,样本集中的早期样本(坏点)经常被丢弃。
模型的收敛性
1)样本轨迹
由于上面模型对 μ 的估计不能表明估计的效果很好,可以采用一些检验手段。首先,观察样本输出,样本的轨迹应该轻微上下浮动,就像一条毛虫。如果你看到轨迹上下跳跃或滞留在某点,这就是模型不收敛的标志,通过 MCMC 采样器得到的估计没有意义。
2)自相关图
也可采另一种测试,自相关测试(见下图),这是评估MCMC采样链中样本前后相关关系的一种方法。相关性弱的样本比相关性强的样本能够给参数估计提供更多的信息。
直观上看,自相关图应该快速衰减到0附近,然后在 0 附近上下波动。如果你的自相关图没有衰减,通常这是弱融合的标志,你需要重新审查你的模型选择(如似然度)和抽样方法(如Metropolis)
_ = pm.autocorrplot(trace[:2000], varnames=['mu'])
—End—
上一篇: 布尔(布林)代数及其运算
下一篇: 求最小化布尔函数 EN 的算法
推荐阅读
-
如何用 python 对生存分析数据样本进行 Weibull 分布参数的最大似然估计 python 参数估计
-
包婷婷 (201550484)作业一 统计软件简介与数据操作-SPSS(Statistical Product and Service Solutions),"统计产品与服务解决方案"软件。最初软件全称为"(SolutionsStatistical Package for the Social Sciences),但是随着SPSS产品服务领域的扩大和服务深度的增加,SPSS公司已于2000年正式将英文全称更改为"统计产品与服务解决方案",标志着SPSS的战略方向正在做出重大调整。为IBM公司推出的一系列用于统计学分析运算、数据挖掘、预测分析和决策支持任务的软件产品及相关服务的总称SPSS,有Windows和Mac OS X等版本。 1984年SPSS总部首先推出了世界上第一个统计分析软件微机版本SPSS/PC+,开创了SPSS微机系列产品的开发方向,极大地扩充了它的应用范围,并使其能很快地应用于自然科学、技术科学、社会科学的各个领域。世界上许多有影响的报刊杂志纷纷就SPSS的自动统计绘图、数据的深入分析、使用方便、功能齐全等方面给予了高度的评价。 R统计软件介绍 R是一套完整的数据处理、计算和制图软件系统。其功能包括:数据存储和处理系统;数组运算工具(其向量、矩阵运算方面功能尤其强大);完整连贯的统计分析工具;优秀的统计制图功能;简便而强大的编程语言:可操纵数据的输入和输出,可实现分支、循环,用户可自定义功能。 与其说R是一种统计软件,还不如说R是一种数学计算的环境,因为R并不是仅仅提供若干统计程序、使用者只需指定数据库和若干参数便可进行一个统计分析。R的思想是:它可以提供一些集成的统计工具,但更大量的是它提供各种数学计算、统计计算的函数,从而使使用者能灵活机动的进行数据分析,甚至创造出符合需要的新的统计计算方法。 该语言的语法表面上类似 C,但在语义上是函数设计语言(functional programming language)的变种并且和Lisp 以及 APL有很强的兼容性。特别的是,它允许在"语言上计算"(computing on the language)。这使得它可以把表达式作为函数的输入参数,而这种做法对统计模拟和绘图非常有用。 R是一个免费的*软件,它有UNIX、LINUX、MacOS和WINDOWS版本,都是可以免费下载和使用的。在R主页那儿可以下载到R的安装程序、各种外挂程序和文档。在R的安装程序中只包含了8个基础模块,其他外在模块可以通过CRAN获得。 二、R语言 R是用于统计分析、绘图的语言和操作环境。R是属于GNU系统的一个*、免费、源代码开放的软件,它是一个用于统计计算和统计制图的优秀工具。 R作为一种统计分析软件,是集统计分析与图形显示于一体的。它可以运行于UNIX,Windows和Macintosh的操作系统上,而且嵌入了一个非常方便实用的帮助系统,相比于其他统计分析软件,R还有以下特点: 1.R是*软件。这意味着它是完全免费,开放源代码的。可以在它的网站及其镜像中下载任何有关的安装程序、源代码、程序包及其源代码、文档资料。标准的安装文件身自身就带有许多模块和内嵌统计函数,安装好后可以直接实现许多常用的统计功能。[2] 2.R是一种可编程的语言。作为一个开放的统计编程环境,语法通俗易懂,很容易学会和掌握语言的语法。而且学会之后,我们可以编制自己的函数来扩展现有的语言。这也就是为什么它的更新速度比一般统计软件,如,SPSS,SAS等快得多。大多数最新的统计方法和技术都可以在R中直接得到。[2] 3. 所有R的函数和数据集是保存在程序包里面的。只有当一个包被载入时,它的内容才可以被访问。一些常用、基本的程序包已经被收入了标准安装文件中,随着新的统计分析方法的出现,标准安装文件中所包含的程序包也随着版本的更新而不断变化。在另外版安装文件中,已经包含的程序包有:base一R的基础模块、mle一极大似然估计模块、ts一时间序列分析模块、mva一多元统计分析模块、survival一生存分析模块等等.[2] 4.R具有很强的互动性。除了图形输出是在另外的窗口处,它的输入输出窗口都是在同一个窗口进行的,输入语法中如果出现错误会马上在窗口口中得到提示,对以前输入过的命令有记忆功能,可以随时再现、编辑修改以满足用户的需要。输出的图形可以直接保存为JPG,BMP,PNG等图片格式,还可以直接保存为PDF文件。另外,和其他编程语言和数据库之间有很好的接口。[2] 5.如果加入R的帮助邮件列表一,每天都可能会收到几十份关于R的邮件资讯。可以和全球一流的统计计算方面的专家讨论各种问题,可以说是全世界最大、最前沿的统计学家思维的聚集地.[2] R是基于S语言的一个GNU项目,所以也可以当作S语言的一种实现,通常用S语言编写的代码都可以不作修改的在R环境下运行。 R的语法是来自Scheme。R的使用与S-PLUS有很多类似之处,这两种语言有一定的兼容性。S-PLUS的使用手册,只要稍加修改就可作为R的使用手册。所以有人说:R,是S-PLUS的一个“克隆”。 但是请不要忘了:R是免费的(R is free)。R语言源代码托管在github,具体地址可以看参考资料。[3] 。 R语言的下载可以通过CRAN的镜像来查找。 R语言有域名为.cn的下载地址,有六个,其中两个由Datagurn,由 中国科学技术大学提供的。R语言Windows版,其中由两个下载地点是Datagurn和 USTC提供的。 三、stata Stata 是一套提供其使用者数据分析、数据管理以及绘制专业图表的完整及整合性统计软件。它提供许许多多功能,包含线性混合模型、均衡重复反复及多项式普罗比模式。用Stata绘制的统计图形相当精美。 新版本的STATA采用最具亲和力的窗口接口,使用者自行建立程序时,软件能提供具有直接命令式的语法。Stata提供完整的使用手册,包含统计样本建立、解释、模型与语法、文献等超过一万余页的出版品。 除此之外,Stata软件可以透过网络实时更新每天的最新功能,更可以得知世界各地的使用者对于STATA公司提出的问题与解决之道。使用者也可以透过Stata. Journal获得许许多多的相关讯息以及书籍介绍等。另外一个获取庞大资源的管道就是Statalist,它是一个独立的listserver,每月交替提供使用者超过1000个讯息以及50个程序。 四、PYTHON