Python – Use the mean and standard deviation of the distribution of scipy.stats

I

am trying to get the mean and standard deviation of the lognormal distribution, where mu=0.4104857306 and sigma=3.4070874277012617, I expect mean=500 and std=600. I’m not sure what I’m doing wrong. Here is the code:

import scipy.stats as stats
import numpy as np
a = 3.4070874277012617
b = 0.4104857306
c = stats.lognorm.mean(a,b)
d = stats.lognorm.var(a,b)
e = np.sqrt(d)
print("Mean =",c)
print("std =",e)

The output is here:

Mean = 332.07447304207903
sd = 110000.50047821256

Thank you in advance.

Edit:

Thanks for your help. I checked it and found that there were some calculation errors. I can get mean=500 now but still can’t get std=600. Here is the code I used:

import numpy as np
import math
from scipy import exp
from scipy.optimize import fsolve

def f(z):
    mean = 500
    std = 600
    sigma = z[0]
    mu = z[1]
    f = np.zeros(2)
    f[0] = exp(mu + (sigma**2) / 2) - mean
    f[1] = exp(2*mu + sigma**2) * exp(sigma**2 - 1) - std**2
    return f
z = fsolve (f,[1.1681794012855686,5.5322865416282365])
print("sigma =",z[0])
print("mu =",z[1])
print(f(z))

sigma = 1.1681794012855686
mu = 5.5322865416282365

I

tried checking the results with my calculator and I could get std=600 as I needed and I still got 853.5698320847896 and lognorm.std(sigma, scale=np.exp (mu)).

scipy.stats.lognorm The lognormal distribution is parameterized in a slightly unusual way in order to be consistent with other continuous distributions. The first parameter is the shape parameter, which is your sigma. Next are the loc and scale parameters, which allow the distribution to be moved and scaled. Here you need loc=0.0 and scale=exp(mu). So, to calculate the mean, you need to do the following:

>>> 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 to be clearer: pass the scale parameter by name and leave the loc parameter at the default value of 0.0:

>>> lognorm.mean(sigma, scale=np.exp(mu))
500.0000010889041

As @coldspeed said in his comment, your expected value for standard deviation does not look correct. I get:

>>> lognorm.std(sigma, scale=np.exp(mu))
165831.2402402415

I calculated the same value by hand.

To double-check whether these parameter choices actually give the expected lognormal distribution, I created a sample containing a million biases and looked at the mean and standard deviation of the logarithms of that sample. As expected, these returned values look roughly similar to 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 editor: you get a slightly wrong formula for the lognormal distribution variance – you need to replace the exp(sigma**2 – 1) term with (exp(sigma**2) - 1). If you do this, and rerun the fsolve calculation, you get:

sigma = 0.9444564779275075
mu = 5.768609079062494

Using these 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

In addition to using fsolve, you can analyze and solve for sigma and mu, given the desired mean and standard deviation. This gives you more accurate results faster:

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