Solución:
Tiene la moda y la desviación estándar de la distribución logarítmica normal. Usar el rvs() método de scipy lognorm, tienes que parametrizar la distribución en términos del parámetro de forma s, que es la desviación estándar sigma de la distribución normal subyacente, y la scale, cual es exp(mu), dónde mu es la media de la distribución subyacente.
Usted señaló que hacer esta reparametrización requiere resolver un polinomio cuártico. Para eso, podemos usar el numpy.poly1d clase. Las instancias de esa clase tienen un roots atributo.
Un poco de álgebra muestra que exp(sigma**2) es la raíz real positiva única del polinomio
x**4 - x**3 - (stddev/mode)**2 = 0
dónde stddev y mode son la desviación estándar dada y la moda de la distribución logarítmica normal, y para esa solución, la scale (es decir exp(mu)) es
scale = mode*x
Aquí hay una función que convierte el modo y la desviación estándar a la forma y la escala:
def lognorm_params(mode, stddev):
"""
Given the mode and std. dev. of the log-normal distribution, this function
returns the shape and scale parameters for scipy's parameterization of the
distribution.
"""
p = np.poly1d([1, -1, 0, 0, -(stddev/mode)**2])
r = p.roots
sol = r[(r.imag == 0) & (r.real > 0)].real
shape = np.sqrt(np.log(sol))
scale = mode * sol
return shape, scale
Por ejemplo,
In [155]: mode = 123
In [156]: stddev = 99
In [157]: sigma, scale = lognorm_params(mode, stddev)
Genere una muestra utilizando los parámetros calculados:
In [158]: from scipy.stats import lognorm
In [159]: sample = lognorm.rvs(sigma, 0, scale, size=1000000)
Aquí está la desviación estándar de la muestra:
In [160]: np.std(sample)
Out[160]: 99.12048952171304
Y aquí hay algo de código matplotlib para trazar un histograma de la muestra, con una línea vertical dibujada en el modo de distribución de la que se extrajo la muestra:
In [176]: tmp = plt.hist(sample, normed=True, bins=1000, alpha=0.6, color="c", ec="c")
In [177]: plt.xlim(0, 600)
Out[177]: (0, 600)
In [178]: plt.axvline(mode)
Out[178]: <matplotlib.lines.Line2D at 0x12c5a12e8>
El histograma:
Si desea generar la muestra usando numpy.random.lognormal() en lugar de scipy.stats.lognorm.rvs(), Puedes hacerlo:
In [200]: sigma, scale = lognorm_params(mode, stddev)
In [201]: mu = np.log(scale)
In [202]: sample = np.random.lognormal(mu, sigma, size=1000000)
In [203]: np.std(sample)
Out[203]: 99.078297384090902
No he investigado cuán robusto poly1d‘s roots algoritmo es, así que asegúrese de probar una amplia gama de posibles valores de entrada. Alternativamente, puede usar un solucionador de scipy para resolver el polinomio anterior para x. Puede unir la solución usando:
max(sqrt(stddev/mode), 1) <= x <= sqrt(stddev/mode) + 1