将数据拟合到分布？

Question

我不是一个统计学家（更像是一个研究性的web开发人员），但这些天我听到了很多关于scipy和R的消息。所以出于好奇，我想问这个问题（尽管在这里的专家看来可能听起来很傻），因为我不确定这方面的进展，想知道没有良好统计背景的人是如何处理这些问题的。

给定一组从实验中观察到的实数，假设它们属于许多分布中的一个（如Weibull、Erlang、Cauchy、Exponential等），有没有自动找到数据的正确分布和分布参数的方法？有什么好的教程可以指导我完成这个过程吗？

现实世界场景： 例如，假设我发起了一项小型调查，记录了一个人每天与300人交谈的人数，我有以下信息：

1 10
2 5
3 20
...
...

其中X Y告诉我，在调查期间X和Y人交谈过。现在利用300人提供的信息，我想把它做成一个模型。问题归结为是否有任何自动方法可以找到正确的数据分布和分布参数？如果没有，是否有一个很好的逐步过程来实现相同的目标？

Answer 1

看看fitdistrplus（http://cran.r-project.org/web/packages/fitdistrplus/index.html）。

有几点需要注意：

• 试试函数descdist，它提供了数据的倾斜与峰度的关系图，还显示了一些常见的分布。
• fitdist允许您拟合任何可以根据密度和cdf定义的分布。
• 然后，您可以使用gofstat来计算KS和AD统计，这些统计可以测量拟合与数据之间的距离。

Answer 2

这可能比你需要的更一般，但可能会给你一些东西继续下去。

从随机数据中估计概率密度函数的一种方法是使用Edgeworth或Butterworth展开。这些近似使用称为cumulants（其无偏估计为k-statistics）的密度函数特性，并将密度函数表示为高斯分布的扰动。

它们都有一些相当可怕的弱点，比如产生发散的密度函数，甚至在某些区域产生负的密度函数。然而，有些人发现它们对于高度聚集的数据，或者作为进一步估计的起点，或者对于分段估计的密度函数，或者作为启发式的一部分，是有用的。

M.G.肯德尔和A.斯图尔特，《统计学高级理论》，第一卷， 查尔斯格里芬，1963年，是我找到的最完整的参考文献，有一整页的篇幅都是关于这个主题的；大多数其他文本最多有一个句子，或是列出了以时间为单位的扩展，而不是累加，这有点没用。不过，祝你好运，我不得不派我的大学图书馆员去档案馆找一份。。。但这是几年前的事了，所以也许今天互联网会更有帮助。

您问题的最一般形式是一个称为非参数密度估计的领域的主题，其中：

• 来自具有未知分布的随机进程的数据，以及
• 对基础流程的约束

…生成一个密度函数，它最有可能生成数据。（更实际地说，您创建了一种方法，用于计算任意给定点处此函数的近似值，您可以使用该方法进行进一步的工作，例如比较两组随机数据的密度函数，以查看它们是否可能来自同一过程）。

不过，就我个人而言，我在使用非参数密度估计做任何有用的事情上运气都不好，但如果你有一个稳定的理智供应，你应该研究一下。

Answer 3

这是一个复杂的问题，没有完美的答案。我将试着给你一个主要概念的概述，并为你指出一些有用的阅读方向。

假设你有一组一维的数据，并且你有一组有限的概率分布函数，你认为这些数据可能是从中产生的。您可以独立地考虑每个分布，并尝试根据您的数据找到合理的参数。 对于给定数据的概率分布函数，有两种设置参数的方法：

1. Least Squares
2. Maximum Likelihood

根据我的经验，近年来，最大似然法被优先考虑，尽管在每个领域可能都不是这样。

下面是一个如何估计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)

这些非常接近正确答案，括号中的数字是参数周围的置信区间。记住，每次生成一组新的点时，都会得到估计值的新答案。

在数学上，这是使用最大似然来估计高斯分布的均值和标准差。似然是指（在这种情况下）“数据给定参数值的概率”。最大似然是指“使生成输入数据的概率最大化的参数值”。最大似然估计是用于查找使生成输入数据的概率最大化的参数值的算法，对于某些分布，它可能涉及numerical optimization算法。在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

指数分布实际上比正态分布产生这些数据的可能性稍大，这可能是因为指数分布不必给负数分配任何概率密度。

当您试图将数据拟合到更多分布时，所有这些估计问题都会变得更糟。具有更多参数的分布更灵活，因此它们比具有更少参数的分布更适合您的数据。另外，有些分布是其他分布的特例（例如，Exponential是Gamma的特例）。因此，使用先验知识将选择模型约束到所有可能模型的子集是非常常见的。

解决参数估计中一些问题的一个技巧是生成大量数据，并将一些数据留给cross-validation。交叉验证你的pa适合性参数到数据，在估计过程中保留一些数据，然后在遗漏的数据上测量每个模型的可能性。