从 R 中的对数正态分布生成随机数的正确方法

Question

我有一组家庭收入数据。我想对数据拟合对数正态分布，然后根据该分布生成随机数。

下面使用 MASS 包的方法似乎给出了不正确的结果。当我使用来自

fitdistr

的参数从对数正态分布创建大样本时，平均值和标准差都远高于原始数据。

library(MASS)

mean(hhi)    # = 124,027.8
sd(hhi)      # = 127,492.1

params <- fitdistr(hhi,'lognormal')  # mu = 11.31322, sigma = 1.016148

random_hhi <- rlnorm(1000000, params$estimate['meanlog'], params$estimate['sdlog'])
mean(random_hhi)    # returns 137,131.5
sd(random_hhi)      # returns 182,139

这篇文章，在对类似问题的评论中建议，提供了一种基于对数正态分布的均值和方差定义的不同方法。这种方法在我的例子中似乎效果更好。请注意，平均值和 SD 更接近数据中的值。

m <- mean(hhi)
s <- sd(hhi)
location <- log(m^2 / sqrt(s^2 + m^2))  # returns 11.3677234406609
shape <- sqrt(log(1 + (s^2 / m^2)))     # returns 0.8491615
random_hhi <- rlnorm(n=1000000, location, shape)
mean(random_hhi)                        # returns 124,024.1
sd(random_hhi)                          # returns 127,291.8

所以，我的问题是：为什么第一种方法失败而第二种方法成功？我希望

fitdistr

返回最适合数据的参数

mu

和

sigma

，以便从具有

rlnorm

的这些参数生成的数字将遵循相同的分布。我错过了什么？

Answer 1

获取这些参数的两种方式并不相同。

MASS::fitdistr

执行最大似然估计，而在第二种情况下，您已采用理论恒等式，就好像它们反映了样本的确切事实一样。

需要明确的是，

rlnorm

（或任何其他

stats

随机采样函数）会准确返回您所要求的分布。通过反转恒等式，您可以简单地验证平均值 = 11.313 和 SD = 1.016 的对数正态分布将大致产生您在原始尺度中获得的参数（不幸的是，对数正态分布相当严重，这使得原始尺度的 SD 很差）总结措施）。问题在于如何获得这些分布参数。

一个小例子：

set.seed(123)
## Doesn't really follow any univariate distribution
x <- c(rnorm(15, 15, 5), rnorm(5, 150, 30))

## What is the lognormal distribution that is most likely to have produced this sample?
MASS::fitdistr(x, "lognormal")
#> logmean=3.292, logsd=1.020

## What is the mean & SD of the lognormal distribution following the observed mean & SD?
c(logmean=log(mean(x)^2 / sqrt(sd(x)^2 + mean(x)^2)),
  logsd=sqrt(log(1 + (sd(x)^2 / mean(x)^2))))
#> logmean=3.429, logsd=0.985

再次强调，这些都是“正确的”，但回答了不同的问题。由此应该清楚为什么第二种方法更接近您观察到的参数。使用正态分布也会发生同样的情况：

MASS::fitdistr(x, "normal")
#> mean=50.142, sd=62.584

c(mean=mean(x), sd=sd(x))
#> mean=50.141, sd=64.210

从这些参数集进行模拟将清楚地为您提供不同的 SD（这就是您作为输入提供的），第二个将重现观察到的 SD，而第一个则不会。

我无法告诉您这两个中的哪一个是您想要的，因为您没有提供更多详细信息，但您似乎不太可能只想准确地重现观察到的摘要统计数据（您已经拥有这些）。如果有的话，观察到的摘要与 MLE 之间的不匹配表明对数正态不完全适合您的数据。

作为附录，这里的代码可以让您重现

MASS::fitdistr

的作用，即如何计算最大似然参数估计值：

mle.normal <- function(x) {
   mu <- mean(x)
   sig <- sqrt(mean((x-mu)**2))
   c(mu=mu, sig=sig)
}

mle.lognormal <- function(x) {
   stopifnot(all(x > 0))
   lx <- log(x)
   lmu <- mean(lx)
   lsig <- sqrt(mean((lx-lmu)**2))
   c(logmu=lmu, logsig=lsig)
}