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

我正在尝试将一些数据拟合到对数正态分布,并使用优化的参数生成随机对数正态分布.
经过一番搜索,我发现了一些解决方案,但没有人说服：

使用fit函数的solution1：

import  numpy as np
from scipy.stats      import lognorm

mydata = [1,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,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,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,'bx')

我不确定我是否理解使用发行版的方式,或者我是否遗漏了其他内容……

我虽然在这里找到解决方案：Does anyone have example code of using scipy.stats.distributions?
但我无法从我的数据中得到形状…我是否在使用fit函数时遗漏了一些东西？

谢谢

编辑：

这是一个例子,以便更好地理解我的问题：

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,scale = lognorm.fit(mydata)
for _ in xrange(1000):
    rnd_log = lognorm.rvs (shape,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:',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,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:',mean(means),sep=' '))
print (paste('observed std :',sdlog,mean(stdes),sep=' '))

我得到：

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

这更接近,所以我猜我在使用scipy时做错了…