It's in my mind a while ago..To get the same values

Copyimport numpy as np
import scipy.stats
ar = np.arange(20)
print(np.std(ar, ddof=1))
print(scipy.stats.tstd(ar))

output #

Copy5.916079783099616
5.916079783099616

My mentor use to say

-->ddof=1 if you're calculating np.std() for a sample taken from your complete dataset.

---> ddof=0 if you're calculating for the full population

The basic reason is that you're taking the standard deviation of two quite different things there. I think you're misunderstanding what scipy.stats.binom does. From the documentation:

The probability mass function for binom is:

binom.pmf(k) = choose(n,k) * p**k * (1-p)**(n-k)

for k in {0,1,...,n}.

binom takes n and p as shape parameters.

When you do binom(189, 100/189), you are creating a distribution that could take on any value from 0 to 189. This distribution unsurprisingly has a much larger variance than the other sample data you're using, which is restricted to values of either zero or one.

It looks like what you want would be scipy.stats.binom(1, 100/189).std(). However, you still can't expect the exact same value as what you're getting with your sample data, because the binom.std is computing the standard deviation of the overall distribution, whereas the other version (scipy.stats.tstd([1]*100 + [0]*89)) is computing the standard deviation only of a sample. If you increase the size of your sample (e.g., do scipy.stats.tstd([1]*1000 + [0]*890)), the sample standard deviation will approach the value you're getting from binom.std.

You can also get the population (not sample) std by using scipy.std or numpy.std instead of scipy.stats.tstd. scipy.stats.tstd doesn't have a ddof option to let you choose the degrees of freedom, and always computes a sample std.

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.9444564779275075
mu = 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