将数据拟合到分布中？

Question

我不是一名统计学家（更倾向于研究型的网页开发人员），但最近我经常听到关于 scipy 和 R 的话题。因此出于好奇心，我想问一个问题（尽管这对于这里的专家来说可能听起来很傻），因为我不确定在这个领域的进展，并想知道没有扎实统计学背景的人如何解决这些问题。

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

现实场景： 例如，假设我进行了一项小型调查并记录了大约300人每天与多少人交谈的信息，我有以下信息：

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

在这里，X Y 表示在调查期间，人 X 与 Y 个人交谈过。现在使用来自这 300 个人的信息，我想将其纳入模型中。问题归结为是否有任何自动化的方法来找到适合此数据的正确分布和分布参数，如果没有，是否有一个好的逐步程序来实现相同的目标？

- Legend

---

你没有描述问题中最重要的部分——你想用这个模型做什么？请说明。 - hadley

这个问题最好在 stats.se 上提问。 - csgillespie

嗯，我愿意允许他不必处理他将如何处理参数模型的问题。甚至仅使用从足够参数模型生成的合成数据就足以提出这样的问题。自助法很棒，但您必须保留或发送数据。 - Iterator

Answer 1

这是一个复杂的问题，没有完美的答案。我将尝试概述主要概念，并为您指引一些有用的阅读材料。假设您有一维数据集，并且有一组有限的概率分布函数，您认为这些数据可能是从中生成的。您可以独立考虑每个分布，并尝试找到合理的参数，以便给定您的数据。有两种方法可以根据数据设置概率分布函数的参数：

1. 最小二乘法
2. 最大似然法

根据我的经验，最大似然法在最近几年中更受欢迎，尽管在每个领域可能并非如此。以下是一个具体的示例，介绍如何在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。为了交叉验证参数对数据的拟合，可以在估计过程中留出一些数据，然后测量每个模型在留出的数据上的可能性。

- James Thompson

---

+1 @James，非常好。你有没有除正态分布之外的二维分布链接或建议？ - denis

当您假设您的概率分布是高斯分布时，最大似然目标函数会简化为（加权）最小二乘函数（即最小二乘法是最大似然的一种特殊情况）。 - Andre Holzner

在交叉验证中，您如何比较不参与训练的数据上不同模型的似然值？即是否有一种统计测试可以接受似然值并说明哪个模型更适合？ - Jacob

@Jacob - 有一些方法，但没有一个完美的。我们很快会遇到过拟合问题，即具有更多参数的模型在拟合数据方面表现更好，但泛化能力不佳。阅读关于AIC和BIC的内容，如果您有更多问题，请在交叉验证Stack Exchange上发布。 - James Thompson

括号中的数字是估计值的标准偏差，而非置信区间。 - James

Answer 2

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

有几个要点需要注意：

• 尝试使用函数 descdist，它提供了数据偏度和峰度的图形，并展示了一些常见的分布。
• fitdist 允许你拟合任何你可以用概率密度函数和累积分布函数定义的分布。
• 然后可以使用 gofstat 计算KS值和AD值，这些值可以测量拟合结果与数据之间的距离。

- Ramnath

---

自从我提出这个问题已经三年了，现在我意识到您的四行代码包含了很多有用的信息。我仔细看了一下，只看了descdist。您能否详细说明如何使用fitdist和gofstat以及descdist来正式呈现分析结果？如果您至少能指引我参考一些现有的教程，我将非常感激。谢谢您的时间！ - Legend

@Ramnath，这里有人需要你的帮助。 - aloha

我也在研究这些函数，你找到了现有的教程吗？ - Uis234

Answer 3

这可能比你需要的更加通用，但可能会给你提供一些启示。从随机数据中估计概率密度函数的一种方法是使用Edgeworth或Butterworth扩展。这些近似使用称为累积量（其无偏估计量为k-统计量）的密度函数属性，并将密度函数表示为与高斯分布的扰动。它们都有一些相当严重的缺点，例如产生发散的密度函数，甚至在某些区域上是负的密度函数。然而，一些人发现它们对于高度聚集的数据，或作为进一步估计的起点，或用于分段估计的密度函数，或作为启发式的一部分非常有用。

"M. G. Kendall和A. Stuart，The advanced theory of statistics，vol. 1，Charles Griffin，1963"是我找到的最全面的参考资料，其中有一个完整的页面专门讨论了这个主题；大多数其他文本最多只有一句话或列出以矩为基础的展开，而这有点无用。不过祝你好运能找到副本，我曾经不得不让我的大学图书管理员前往档案馆...但那是几年前的事了，也许今天互联网会更有帮助。

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

• 来自具有未知分布的随机过程的数据，以及
• 对底层过程的约束

...你产生一个密度函数，它最有可能产生这些数据。（更现实地说，你创建一种方法，在任何给定点上计算该函数的近似值，你可以用它进行进一步的工作，例如比较两组随机数据的密度函数，以查看它们是否可能来自同一个过程）。

"

然而，就我个人而言，使用非参数密度估计并没有什么实际用处，但如果你有足够的理智，可以研究一下。

- detly

Answer 4

我不是科学家，但如果你用铅笔和纸做的话，显而易见的方法就是绘制图表，然后将其与已知标准分布之一进行比较。

进一步来说，“比较”就是要看标准分布曲线和你的曲线是否相似。

三角学、正切等可能是我的最后想法。

我不是专家，只是一名谦卑的Web开发者 =）

- TomWilde

---

我是一名科学家，你构建数据图并将其与众所周知的分布进行比较的想法真的很好 - 它是最大似然和最小二乘拟合的基础。两者之间的区别在于它们如何评估您的数据和分布之间的适合度，但两种方法都基于你直观有吸引力的想法。 :) - James Thompson

Answer 5

您想要将您的真实世界数据与一组理论分布进行比较。在基本R中有qqnorm()函数，它可以对正态分布进行这样的操作，但我更喜欢e1071中的probplot函数，它允许您测试其他分布。下面是一个代码片段，它将把您的真实数据与我们粘贴到列表中的每个理论分布进行比较。我们使用plyr来遍历该列表，但也有其他几种方法可以遍历该列表。

library("plyr") 
library("e1071")

realData <- rnorm(1000) #Real data is normally distributed

distToTest <- list(qnorm = "qnorm", lognormal = "qlnorm", qexp =  "qexp")

#function to test real data against list of distributions above. Output is a jpeg for each distribution.
testDist <- function(x, data){
    jpeg(paste(x, ".jpeg", sep = ""))
    probplot(data, qdist = x)
    dev.off()
    }

l_ply(distToTest, function(x) testDist(x, realData))

- Chase

---

请问是否可以将“负二项式”分布添加到测试列表中？我一直在尝试，但不确定如何添加有空格的东西。例如，R网站上的参考资料表明我需要输入Negative Binomial，但不确定如何将其添加到列表中。 - Legend

Answer 6

就算只是一点点价值，看起来你可能想要研究泊松分布。

- Ghalib Suleiman