使用观测数据的形状生成随机对数正态分布

Question

我正在尝试把一些数据拟合到对数正态分布上，然后用优化后的参数生成随机的对数正态分布。经过一番搜索，我找到了一些解决方案，但都不是很令人信服：

解决方案1：使用拟合函数：

import  numpy as np
from scipy.stats      import lognorm

mydata = [1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,3,3,3,3,3,3,3,3,3,3,3,4,4,4,4,4,4,4,4,5,5,5,5,5,6,6,6,6,6,7,7,7,8,8,8,8,8,9,9,9,10,10,11,12,13,14,14,15,19,19,21,23,25,27,28,30,31,36,41,45,48,52,55,60,68,75,86,118,159,207,354]

shape, loc, scale = lognorm.fit(mydata)
rnd_log = lognorm.rvs (shape, loc=loc, scale=scale, size=100)

或者解决方案2：使用原始数据中的mu和sigma：

import  numpy as np
from scipy.stats      import lognorm

mydata = [1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,3,3,3,3,3,3,3,3,3,3,3,4,4,4,4,4,4,4,4,5,5,5,5,5,6,6,6,6,6,7,7,7,8,8,8,8,8,9,9,9,10,10,11,12,13,14,14,15,19,19,21,23,25,27,28,30,31,36,41,45,48,52,55,60,68,75,86,118,159,207,354]

mu    = np.mean([np.log(i) for i in mydata])
sigma = np.std([np.log(i) for i in mydata])

distr   = lognorm(mu, sigma)
rnd_log = distr.rvs (size=100)

这些解决方案都没有很好地拟合：

import pylab
pylab.plot(sorted(mydata, reverse=True), 'ro')
pylab.plot(sorted(rnd_log, reverse=True), 'bx')

我不确定自己是否理解了如何使用分布，或者我是不是漏掉了什么...

我以为可以在这里找到解决方案：有没有人有使用scipy.stats.distributions的示例代码？，但我无法从我的数据中得到形状...我在使用拟合函数时是不是漏掉了什么？

谢谢

编辑：

这是一个例子，帮助我更好地理解我的问题：

print 'solution 1:'
means = []
stdes = []
distr   = lognorm(mu, sigma)
for _ in xrange(1000):
    rnd_log = distr.rvs (size=100)
    means.append (np.mean([np.log(i) for i in rnd_log]))
    stdes.append (np.std ([np.log(i) for i in rnd_log]))
print 'observed mean:',mu   , 'mean simulated mean:', np.mean (means)
print 'observed std :',sigma, 'mean simulated std :', np.mean (stdes)

print '\nsolution 2:'
means = []
stdes = []
shape, loc, scale = lognorm.fit(mydata)
for _ in xrange(1000):
    rnd_log = lognorm.rvs (shape, loc=loc, scale=scale, size=100)
    means.append (np.mean([np.log(i) for i in rnd_log]))
    stdes.append (np.std ([np.log(i) for i in rnd_log]))
print 'observed mean:',mu   , 'mean simulated mean:', np.mean (means)
print 'observed std :',sigma, 'mean simulated std :', np.mean (stdes)

结果是：

solution 1:
observed mean: 1.82562655734 mean simulated mean: 1.18929982267
observed std : 1.39003773799 mean simulated std : 0.88985924363

solution 2:
observed mean: 1.82562655734 mean simulated mean: 4.50608084668
observed std : 1.39003773799 mean simulated std : 5.44206119499

而如果我在R中做同样的事情：

mydata <- c(1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,3,3,3,3,3,3,3,3,3,3,3,4,4,4,4,4,4,4,4,5,5,5,5,5,6,6,6,6,6,7,7,7,8,8,8,8,8,9,9,9,10,10,11,12,13,14,14,15,19,19,21,23,25,27,28,30,31,36,41,45,48,52,55,60,68,75,86,118,159,207,354)
meanlog <- mean(log(mydata))
sdlog <- sd(log(mydata))
means <- c()
stdes <- c()
for (i in 1:1000){
  rnd.log <- rlnorm(length(mydata), meanlog, sdlog)
  means <- c(means, mean(log(rnd.log)))
  stdes <- c(stdes, sd(log(rnd.log)))
}

print (paste('observed mean:',meanlog,'mean simulated mean:',mean(means),sep=' '))
print (paste('observed std :',sdlog  ,'mean simulated std :',mean(stdes),sep=' '))

我得到：

[1] "observed mean: 1.82562655733507 mean simulated mean: 1.82307191072317"
[1] "observed std : 1.39704049131865 mean simulated std : 1.39736545866904"

这要接近得多，所以我想我在使用scipy时做错了什么...

Answer 1

在scipy这个库里，lognormal分布的参数设置方式跟我们通常见到的有点不同。你可以查看一下scipy.stats.lognorm的文档，特别是里面的“注意事项”部分。

下面是如何得到你想要的结果的方法（注意在拟合的时候我们把位置参数设为0）：

In [315]: from scipy import stats

In [316]: x = np.array([1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,3,3,3,3,3,3,3,3,3,3,3,4,4,4,4,4,4,4,4,5,5,5,5,5,6,6,6,6,6,7,7,7,8,8,8,8,8,9,9,9,10,10,11,12,13,14,14,15,19,19,21,23,25,27,28,30,31,36,41,45,48,52,55,60,68,75,86,118,159,207,354])

In [317]: mu, sigma = stats.norm.fit(np.log(x))

In [318]: mu, sigma
Out[318]: (1.8256265573350701, 1.3900377379913127)

In [319]: shape, loc, scale = stats.lognorm.fit(x, floc=0)

In [320]: np.log(scale), shape
Out[320]: (1.8256267737298788, 1.3900309739954713)

现在你可以生成样本来确认你的预期结果了：

In [321]: dist = stats.lognorm(shape, loc, scale)

In [322]: means, sds = [], []

In [323]: for i in xrange(1000):
   .....:     sample = dist.rvs(size=100)
   .....:     logsample = np.log(sample)
   .....:     means.append(logsample.mean())
   .....:     sds.append(logsample.std())
   .....:

In [324]: np.mean(means), np.mean(sds)
Out[324]: (1.8231068508345041, 1.3816361818739145)