We used the log-normal distribution configured to get the desired skewness and variance, then shifted it by its mean to center it at 0.
import numpy as np
import matplotlib.pyplot as plt
import scipy.stats
import sympy
from sympy import exp, ln, pi, sqrt
x, mu, sigma = sympy.symbols('x mu sigma')
# Equations for the mean, variance, and skewness of the log-normal distribution.
lognormal_mean_eqn = exp(mu + sigma**2 / 2)
lognormal_variance_eqn = (exp(sigma**2) - 1) * exp(2*mu + sigma**2)
lognormal_skew_eqn = (exp(sigma**2) + 2) * sqrt(exp(sigma**2) - 1)
# We set the skewness equation for the log-normal distribution to 2.7, then solved
# for sigma to get this value.
lognormal_sigma = 0.672337
# We then plugged sigma into the variance equation for the log-normal distribution
# and set it equal to 4.2 and solved for mu.
lognormal_mu = 0.77126
lognormal_params = {mu: lognormal_mu, sigma: lognormal_sigma}
# Plug the parameters we found into the variance/skewness equations to confirm
# they produce the desired values.
print 'variance', lognormal_variance_eqn.subs(lognormal_params)
print 'skewness', lognormal_skew_eqn.subs(lognormal_params)
variance 4.19999788269716 skewness 2.70000063930158
# The variance and skewness look good, but the mean is not 0.
mean = lognormal_mean_eqn.subs(lognormal_params)
print mean
2.71089413824342
# The log-normal PDF can be modified to shift it by an amount equal to the mean
# we calculated.
def get_shifted_lognormal_pdf(mu, sigma, shift=0.0):
"""Returns an expression for the shifted PDF of a log-normal distribution.
This is the standard equation for the PDF of a log-normal distribution, but
shifted by 'shift' to the right.
"""
normalizer = 1 / ((x - shift) * sqrt(2*pi) * sigma)
exp_part = exp(-((ln(x - shift) - mu)**2)/(2 * sigma**2))
pdf = normalizer * exp_part
return pdf
# Get the PDF with the desired mean/variance/skewness.
mean_shift = -mean
pdf = get_shifted_lognormal_pdf(lognormal_mu, lognormal_sigma, shift=mean_shift)
print pdf
0.743674675051351*2**(1/2)*exp(-1.10610404462546*(log(x + 2.71089413824342) - 0.77126)**2)/(pi**(1/2)*(x + 2.71089413824342))
# Plot the PDF
xs = np.linspace(mean_shift + 1e-9, 10, num=1000)
vectorized_pdf = np.vectorize(sympy.lambdify(x, pdf))
ys = vectorized_pdf(xs)
plt.plot(xs, ys)
plt.axvline(x=0, color='red', linestyle='--')
<matplotlib.lines.Line2D at 0xa75a88c>