数据拟合分布吗?
我不是统计学家,更多的是一个研究型的网页开发者,但最近我听说很多关于 scipy 和 R 的事情。出于好奇,我想问这个问题(虽然对这里的专家来说可能听起来有点傻),因为我对这个领域的进展不太了解,想知道那些没有扎实统计背景的人是怎么处理这些问题的。
假设我们有一组从实验中观察到的真实数字,这些数字可能属于很多种分布中的一种(比如威布尔分布、厄朗分布、柯西分布、指数分布等等),有没有什么自动化的方法可以找到合适的分布和数据的分布参数?有没有好的教程可以带我一步步了解这个过程?
现实场景: 比如说,我进行了一项小调查,记录了每个人每天和多少人交谈的信息,假设我调查了300个人,得到了以下信息:
1 10
2 5
3 20
...
...
这里的X和Y告诉我,X这个人和Y个人在调查期间进行了交谈。现在,我想利用这300个人的信息来建立一个模型。问题就是,是否有自动化的方法可以找出这个数据的合适分布和分布参数?如果没有,是否有一个好的逐步流程可以实现这个目标?
6 个回答
这可能有点超出你的需求,但也许能给你一些启发。
从随机数据中估计概率密度函数的一种方法是使用边缘展开(Edgeworth)或巴特沃斯展开(Butterworth)。这些方法利用了一些叫做 累积量 的特性(它们的无偏估计是 k-统计量),并把密度函数表示为从高斯分布(正态分布)的小扰动。
不过,这些方法有一些严重的缺点,比如可能会产生发散的密度函数,甚至在某些区域出现负值的密度函数。不过,有些人发现它们在数据高度聚集时还是有用的,或者可以作为进一步估计的起点,或者用于分段估计的密度函数,或者作为启发式方法的一部分。
M. G. Kendall 和 A. Stuart 的《统计学高级理论》,第一卷,Charles Griffin,1963年 是我找到的最完整的参考资料,专门有一整页讲这个话题;而其他大多数书籍最多只提到一句,或者是用矩而不是累积量来列出展开式,这样就没什么用处了。希望你能找到一本,不过我当时还得让我的大学图书馆员去档案馆找... 不过那是很多年前的事了,也许现在在网上能找到更多帮助。
你问题的最一般形式属于一个叫做 非参数密度估计 的领域,给定:
- 来自未知分布的随机过程的数据,以及
- 对基础过程的约束
...你可以生成一个最有可能产生这些数据的密度函数。(更现实地说,你创建了一种方法,可以在任何给定点计算这个函数的近似值,这样你就可以用它进行进一步的工作,比如比较两组随机数据的密度函数,看看它们是否可能来自同一个过程。)
不过,个人来说,我在使用非参数密度估计方面并没有太多成功,但如果你有足够的耐心,可以考虑研究一下。
看看 fitdistrplus
这个包吧,链接在这里:http://cran.r-project.org/web/packages/fitdistrplus/index.html。
有几个简单的要点需要注意:
- 可以试试
descdist
这个函数,它会给你画出数据的偏斜度和峰度的图,还会显示一些常见的分布情况。 fitdist
这个函数让你可以拟合任何你能用密度和累积分布函数(cdf)定义的分布。- 然后你可以使用
gofstat
,这个函数会计算KS和AD统计量,用来衡量拟合结果和数据之间的差距。
这个问题比较复杂,没有完美的答案。我会尽量给你一个主要概念的概述,并推荐一些有用的阅读材料。
假设你有一组一维的数据,并且你认为这些数据可能是从一组有限的概率分布中生成的。你可以独立地考虑每个分布,并尝试找到一些合理的参数来描述你的数据。
根据数据设置概率分布函数的参数有两种方法:
根据我的经验,近年来最大似然法更受欢迎,尽管在某些领域可能并非如此。
下面是一个在R语言中估计参数的具体例子。考虑一组从均值为0、标准差为1的高斯分布生成的随机点:
x = rnorm( n = 100, mean = 0, sd = 1 )
假设你知道这些数据是通过高斯过程生成的,但你忘记了(或者从未知道过!)高斯的参数。你想用这些数据来合理估计均值和标准差。在R中,有一个标准库可以让这个过程变得非常简单:
library(MASS)
params = fitdistr( x, "normal" )
print( params )
这给我输出了以下结果:
mean sd
-0.17922360 1.01636446
( 0.10163645) ( 0.07186782)
这些结果相当接近正确答案,括号中的数字是参数的置信区间。记住,每次你生成一组新的点时,估计的结果都会不同。
从数学上讲,这使用最大似然法来估计高斯的均值和标准差。似然性在这里的意思是“在给定参数值的情况下,数据出现的概率”。最大似然性则是“使得生成我输入数据的概率最大的参数值”。最大似然估计是一个算法,用于寻找能够最大化生成输入数据概率的参数值,对于某些分布,这可能涉及到一些数值优化算法。在R中,大部分工作是由fitdistr完成的,在某些情况下会调用optim。
你可以像这样提取参数的对数似然性:
print( params$loglik )
[1] -139.5772
通常我们更常用对数似然性而不是似然性,以避免舍入误差。估计数据的联合概率涉及到相乘概率,而这些概率都小于1。即使是小数据集,联合概率也会很快接近0,而将数据的对数概率相加等同于相乘概率。当对数似然性接近0时,似然性达到最大,因此更负的数字对数据的拟合效果更差。
有了这样的计算工具,估计任何分布的参数变得很简单。考虑这个例子:
x = x[ x >= 0 ]
distributions = c("normal","exponential")
for ( dist in distributions ) {
print( paste( "fitting parameters for ", dist ) )
params = fitdistr( x, dist )
print( params )
print( summary( params ) )
print( params$loglik )
}
指数分布不会生成负数,所以我在第一行中去掉了负数。输出(这是随机的)看起来像这样:
[1] "fitting parameters for normal"
mean sd
0.72021836 0.54079027
(0.07647929) (0.05407903)
Length Class Mode
estimate 2 -none- numeric
sd 2 -none- numeric
n 1 -none- numeric
loglik 1 -none- numeric
[1] -40.21074
[1] "fitting parameters for exponential"
rate
1.388468
(0.196359)
Length Class Mode
estimate 1 -none- numeric
sd 1 -none- numeric
n 1 -none- numeric
loglik 1 -none- numeric
[1] -33.58996
实际上,指数分布生成这些数据的可能性比正态分布稍高,这可能是因为指数分布不需要给负数分配任何概率密度。
当你尝试将数据拟合到更多分布时,所有这些估计问题都会变得更糟。参数更多的分布更灵活,因此它们会比参数较少的分布更好地拟合你的数据。此外,一些分布是其他分布的特例(例如,指数分布是伽马分布的特例)。因此,使用先前的知识来限制你的模型选择在所有可能模型的子集内是非常常见的。
解决参数估计中一些问题的一个技巧是生成大量数据,并将一些数据留出用于交叉验证。为了交叉验证你的参数与数据的拟合,留出一些数据不参与估计过程,然后在留出的数据上测量每个模型的似然性。