Sample from a distribution with known skew and kurtosis

sample_higher_moments.py

import statsmodels.sandbox.distributions.extras as extras
import scipy.interpolate as interpolate
import numpy as np

def generate_normal_four_moments(mu, sigma, skew, kurt, size=1000, sd_wide = 10):
    f = extras.pdf_mvsk([mu, sigma, skew, kurt])
    x = np.linspace(mu - sd_wide * sigma, mu + sd_wide * sigma, num=500)
    y = [f(i) for i in x]
    yy = np.cumsum(y) / np.sum(y)
    inv_cdf = interpolate.interp1d(yy, x, fill_value="extrapolate")
    rr = np.random.rand(size)
    return inv_cdf(rr)

just for anyone else who stumbles upon this. you're use of the term sigma is somewhat misleading. the function asks for v (variance). should probably be sigma ** 2.

this gist comes up on searches for use of pdf_mvsk so thought i'd point this out. @raddy, feel free to disagree.

I think you're correct looking at the pdf_mvsk implementation:

    N = len(mvsk)
    if N < 4:
        raise ValueError("Four moments must be given to "
                         "approximate the pdf.")

    mu, mc2, skew, kurt = mvsk

    totp = poly1d(1)
    sig = sqrt(mc2)
    if N > 2:
        Dvals = _hermnorm(N + 1)
        C3 = skew / 6.0
        C4 = kurt / 24.0
        # Note: Hermite polynomial for order 3 in _hermnorm is negative
        # instead of positive
        totp = totp - C3 * Dvals[3] + C4 * Dvals[4]

    def pdffunc(x):
        xn = (x - mu) / sig
        return totp(xn) * np.exp(-xn * xn / 2.0) / np.sqrt(2 * np.pi) / sig

    return pdffunc

Let me play with some numbers to sanity check!