Drawing from a Laplace distribution using Scipy yields negatively skewed density-CodePudding

When I make draws from a Laplace distribution with mean zero and scale drawn from any distribution that maps into the positive orthant, the resulting empirical distribution is negatively skewed, regardless of the number of draws, distribution for the scale and seed. Regarding the large sample size symmetry is expected, however. See the following two examples that can be reproduced

Example 0:

import numpy as np
from scipy.stats import halfcauchy
from scipy.stats import laplace

lam_0 = halfcauchy.rvs(loc=0, scale=1, size=2000000, random_state=77)
lap_0 = laplace.rvs(loc=0, scale=1 / lam_0, random_state=77)
np.quantile(lap_0, 0.05)
-22.130260524443447
np.quantile(lap_0, 0.95)
0.38451887570738214
np.mean(lap_0) 
-67.06943091954444

Example 1:

import numpy as np
from scipy.stats import expon
from scipy.stats import laplace

lam_1 = expon.rvs(loc=0, scale=2, size=1000000, random_state=42)
lap_1 = laplace.rvs(loc=0, scale=1 / lam_1, random_state=42)
np.quantile(lap_1, 0.05)
-29.27074349002619
np.quantile(lap_1, 0.95)
0.2953765780255653
np.mean(lap_1) 
-71.64564905737133

CodePudding user response：

In both examples, you are giving the same integer random_state argument to the distribution. So, for example, the generation of lam_0 and lap_0 are based on the same sequence of samples from the uniform distribution generated by the underlying random number generator. This results in a correlation between the array that you give as the scale argument to laplace.rvs() and the uniform samples that laplace.rvs() uses to generate its samples--and that results in the negative values being large and the positive values being small.

If you use different random_state parameters, you get the distribution that you expect. If you want reproducibility, you can create one generator instance, and pass it to both functions (so the underlying generator is generating new uniform samples in the call of laplace.rvs()). For example,

In [252]: rng = np.random.default_rng(121263137472525314065)

In [253]: lam_0 = halfcauchy.rvs(loc=0, scale=1, size=2000000, random_state=rng)

In [254]: lap_0 = laplace.rvs(loc=0, scale=1/lam_0, random_state=rng)

In [255]: np.quantile(lap_0, 0.05)
Out[255]: -6.064955243909569

In [256]: np.quantile(lap_0, 0.95)
Out[256]: 6.111148574695757