如何在 R 中生成给定、均值、SD、偏斜和峰度的分布？

Posted 2023-02-16

【中文标题】如何在 R 中生成给定、均值、SD、偏斜和峰度的分布？

【英文标题】：How to generate distributions given, mean, SD, skew and kurtosis in R?

【发布时间】：2011-06-15 23:55:38

【问题描述】：

是否可以在 R 中生成均值、SD、偏斜和峰度已知的分布？到目前为止，最好的方法似乎是创建随机数并相应地转换它们。 如果有专门用于生成可以调整的特定发行版的软件包，我还没有找到。 谢谢

【问题讨论】：

如前所述，这些并不能唯一地描述一个分布。即使你定义了所有的时刻，你也不能保证唯一地定义一个分布。我认为你需要解释你到底想要做什么。你为什么要这样做？您能否设置更多限制以定义分布？

啊，是的，我们想要一个单一维度的单峰连续分布。由此产生的分布最终将在数值上进行转换，以通过模拟测试生态位理论的变化。

在交叉验证 (stats.SE) 上，以下内容有些相关，读者可能对此感兴趣：How to simulate data that satisfy specific constraints such as having specific mean and standard deviation?

【参考方案1】：

熵方法是个好主意，但是如果您有数据样本，则与仅使用矩相比，您使用的信息更多！因此，力矩拟合通常不太稳定。如果您没有关于分布的更多信息，那么熵是一个很好的概念，但是如果您有更多信息，例如关于支持，然后使用它！如果您的数据偏斜且为正，那么使用对数正态模型是一个好主意。如果您也知道上尾是有限的，那么不要使用对数正态分布，而可能使用 4 参数 Beta 分布。如果对支持或尾部特征一无所知，那么缩放和移位的对数正态模型可能很好。如果您需要更多关于峰度的灵活性，那么例如具有缩放 + 移位的 logT 通常很好。如果您知道拟合应该接近正态，它也会有所帮助，如果是这种情况，则使用包含正态分布的模型（无论如何通常都是如此），否则您可能会例如使用广义割线双曲线分布。如果你想做这一切，那么在某些时候模型会有一些不同的情况，你应该确保没有间隙或不良的过渡效果。

【参考方案2】：

您的一个解决方案可能是 PearsonDS 库。它允许您使用前四个矩与峰度 > 偏度^2 + 1 的限制的组合。

要从该分布中生成 10 个随机值，请尝试：

library("PearsonDS")
moments <- c(mean = 0,variance = 1,skewness = 1.5, kurtosis = 4)
rpearson(10, moments = moments)

【参考方案3】：

这个问题是 3 年前提出的，所以我希望我的回答不会太晚。

有一种方法可以在知道某些时刻时唯一地识别分布。这就是最大熵的方法。这种方法产生的分布是最大化你对分布结构的无知的分布，考虑到你所知道的。任何其他也具有您指定的矩但不是 MaxEnt 分布的分布都隐含地假设比您输入的结构更多。最大化的泛函是香农的信息熵，$S[p(x)] = - \int p(x)log p(x) dx$。知道均值、sd、偏度和峰度，分别转化为对分布的一阶、二阶、三阶和四阶矩的约束。

然后问题是最大化 S 受约束： 1) $\int x p(x) dx = "第一时刻"$, 2) $\int x^2 p(x) dx = "第二时刻"$, 3) ... 等等

我推荐“Harte, J., Maximum Entropy and Ecology: A Theory of Abundance, Distribution, and Energetics (Oxford University Press, New York, 2011)”一书。

这是一个尝试在 R 中实现此功能的链接： https://stats.stackexchange.com/questions/21173/max-entropy-solver-in-r

【参考方案4】：

我同意您需要密度估计来复制任何分布。但是，如果您有数百个变量，这在 Monte Carlo 模拟中很典型，您就需要做出妥协。

一种建议的方法如下：

使用 Fleishman 变换获取给定偏斜和峰度的系数。 Fleishman 获取偏斜和峰度并为您提供系数 生成 N 个正态变量（均值 = 0，标准值 = 1） 用 Fleishman 系数转换 (2) 中的数据，将正常数据转换为给定的偏斜和峰度 在此步骤中，使用来自第 (3) 步的数据，并使用 new_data = 所需均值 +（第 3 步中的数据）* 所需标准差将其转换为所需的均值和标准差 (std)

第 4 步得到的数据将具有所需的均值、标准差、偏度和峰度。

注意事项：

Fleishman 不适用于偏度和峰度的所有组合 以上步骤假设不相关变量。如果要生成关联数据，则需要在 Fleishman 变换之前进行一步

【讨论】：

这个有R实现吗？

【参考方案5】：

SuppDists 包中有一个 Johnson 发行版。 Johnson 会给你一个匹配矩或分位数的分布。其他 cmets 是正确的，即 4 个矩不是一个分布。但约翰逊肯定会尝试。

这是一个将 Johnson 拟合到一些样本数据的示例：

require(SuppDists)

## make a weird dist with Kurtosis and Skew
a <- rnorm( 5000, 0, 2 )
b <- rnorm( 1000, -2, 4 )
c <- rnorm( 3000,  4, 4 )
babyGotKurtosis <- c( a, b, c )
hist( babyGotKurtosis , freq=FALSE)

## Fit a Johnson distribution to the data
## TODO: Insert Johnson joke here
parms<-JohnsonFit(babyGotKurtosis, moment="find")

## Print out the parameters 
sJohnson(parms)

## add the Johnson function to the histogram
plot(function(x)dJohnson(x,parms), -20, 20, add=TRUE, col="red")

最终的情节是这样的：

您可以看到其他人指出的一些问题，即 4 个矩不能完全捕获分布。

祝你好运！

编辑 正如哈德利在 cmets 中指出的那样，约翰逊的合身看起来不太好。我做了一个快速测试并使用moment="quant" 拟合约翰逊分布，它使用 5 个分位数而不是 4 个矩来拟合约翰逊分布。结果看起来好多了：

parms<-JohnsonFit(babyGotKurtosis, moment="quant")
plot(function(x)dJohnson(x,parms), -20, 20, add=TRUE, col="red")

这会产生以下内容：

任何人都知道为什么约翰逊在使用时刻适合时似乎有偏见？

【讨论】：

那条曲线看起来有些不对劲 - 一个简单的位置偏移将使拟合效果显着提高

我同意它看起来不错。当我有一点时间时，我可能会深入研究它。

【参考方案6】：

这是一个有趣的问题，实际上并没有很好的解决方案。我认为即使你不知道其他时刻，你也知道分布应该是什么样子。例如，它是单峰的。

有几种不同的方法可以解决这个问题：

假设一个基础分布和匹配时刻。有许多标准的 R 包可以做到这一点。一个缺点是多元泛化可能不清楚。

鞍点近似。在本文中：

Gillespie, C.S. 和 Renshaw, E. An improved saddlepoint approximation. 数学生物科学，2007。

我们只考虑在最初的几分钟内恢复 pdf/pmf。我们发现这种方法在偏度不太大的情况下有效。

拉盖尔扩展：

Mustapha, H. 和 Dimitrakopoulosa, R. Generalized Laguerre expansions of multivariate probability densities with moments。 计算机和数学及其应用，2010 年。

本文中的结果似乎更有希望，但我没有将它们编码。

【参考方案7】：

正如@David 和@Carl 上面所写，有几个包专门用于生成不同的发行版，参见例如the Probability distributions Task View on CRAN.

如果您对理论感兴趣（如何使用给定参数绘制适合特定分布的数字样本），那么只需寻找适当的公式，例如请参阅gamma distribution on Wiki，并使用提供的参数组成一个简单的质量系统来计算比例和形状。

查看一个具体示例 here，其中我根据均值和标准差计算了所需 beta 分布的 alpha 和 beta 参数。

【参考方案8】：

这些参数实际上并未完全定义分布。为此，您需要一个密度或等效的分布函数。