Distribution Mean And Standard Deviation Using Scipy.stats
Solution 1:
The scipy.stats.lognorm
lognormal distribution is parameterised in a slightly unusual way, in order to be consistent with the other continuous distributions. The first argument is the shape parameter, which is your sigma
. That's followed by the loc
and scale
arguments, which allow shifting and scaling of the distribution. Here you want loc=0.0
and scale=exp(mu)
. So to compute the mean, you want to do something like:
>>>import numpy as np>>>from scipy.stats import lognorm>>>mu = 0.4104857306>>>sigma = 3.4070874277012617>>>lognorm.mean(sigma, 0.0, np.exp(mu))
500.0000010889041
Or more clearly: pass the scale
parameter by name, and leave the loc
parameter at its default of 0.0
:
>>>lognorm.mean(sigma, scale=np.exp(mu))
500.0000010889041
As @coldspeed says in his comment, your expected value for the standard deviation doesn't look right. I get:
>>> lognorm.std(sigma, scale=np.exp(mu))
165831.2402402415
and I get the same value calculating by hand.
To double check that these parameter choices are indeed giving the expected lognormal, I created a sample of a million deviates and looked at the mean and standard deviation of the log of that sample. As expected, those give me back values that look roughly like your original mu
and sigma
:
>>>samples = lognorm.rvs(sigma, scale=np.exp(mu), size=10**6)>>>np.log(samples).mean() # should be close to mu
0.4134644116056518
>>>np.log(samples).std(ddof=1) # should be close to sigma
3.4050012251732285
In response to the edit: you've got the formula for the variance of a lognormal slightly wrong - you need to replace the exp(sigma**2 - 1)
term with (exp(sigma**2) - 1)
. If you do that, and rerun the fsolve
computation, you get:
sigma = 0.9444564779275075mu = 5.768609079062494
And with those values, you should get the expected mean and standard deviation:
>>>from scipy.stats import lognorm>>>import numpy as np>>>sigma = 0.9444564779275075>>>mu = 5.768609079062494>>>lognorm.mean(sigma, scale=np.exp(mu))
499.9999999949592
>>>lognorm.std(sigma, scale=np.exp(mu))
599.9999996859631
Rather than using fsolve
, you could also solve analytically for sigma
and mu
, given the desired mean and standard deviation. This gives you more accurate results, more quickly:
>>>mean = 500.0>>>var = 600.0**2>>>sigma = np.sqrt(np.log1p(var/mean**2))>>>mu = np.log(mean) - 0.5*sigma*sigma>>>mu, sigma
(5.768609078769636, 0.9444564782482624)
>>>lognorm.mean(sigma, scale=np.exp(mu))
499.99999999999966
>>>lognorm.std(sigma, scale=np.exp(mu))
599.9999999999995
Post a Comment for "Distribution Mean And Standard Deviation Using Scipy.stats"