создание случайных логнормальных распределений по форме наблюдаемых данных

Я пытаюсь подогнать некоторые данные под логнормальное распределение и на основе этого сгенерировать случайное логнормальное распределение с оптимизированными параметрами. После некоторого поиска я нашел несколько решений, но ни одно не убедительно:

решение 1 с использованием функции fit:

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? , но я не могу получить форму из моих данных... я что-то упустил в использовании функции fit?

спасибо

EDIT:

это пример, чтобы лучше понять мою проблему:

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...