erzeugen N Zufallszahlen aus einer Skew-Normalverteilung mit numpy

Question

Ich brauche eine Funktion in Python, um N zufällige Zahlen aus einer Skew-Normalverteilung zurückzugeben. Der Skew muss als Parameter genommen werden.

z.B. meine aktuelle Nutzung ist

x = numpy.random.randn(1000)

und die ideale Funktion wäre z.B.

x = randn_skew(1000, skew=0.7)

Lösung muss mit entsprechen: Python-Version 2.7, numpy v.1.9

Eine ähnliche Antwort ist hier: skew normal distribution in scipy jedoch nicht die Zufallszahlen ein PDF erzeugt.

Answer 1

Ich beginne mit der PDF-Kurven als Referenz zu erzeugen:

NUM_SAMPLES = 100000 
SKEW_PARAMS = [-3, 0] 

def skew_norm_pdf(x,e=0,w=1,a=0): 
    # adapated from: 
    # http://stackoverflow.com/questions/5884768/skew-normal-distribution-in-scipy 
    t = (x-e)/w 
    return 2.0 * w * stats.norm.pdf(t) * stats.norm.cdf(a*t) 

# generate the skew normal PDF for reference: 
location = 0.0 
scale = 1.0 
x = np.linspace(-5,5,100) 

plt.subplots(figsize=(12,4)) 
for alpha_skew in SKEW_PARAMS: 
    p = skew_norm_pdf(x,location,scale,alpha_skew) 
    # n.b. note that alpha is a parameter that controls skew, but the 'skewness' 
    # as measured will be different. see the wikipedia page: 
    # https://en.wikipedia.org/wiki/Skew_normal_distribution 
    plt.plot(x,p) 

Next ich eine VB Implementierung des Abtastens Zufallszahlen aus der Skew Normalverteilung gefunden und konvertiert es zu Python:

# literal adaption from: 
# http://stackoverflow.com/questions/4643285/how-to-generate-random-numbers-that-follow-skew-normal-distribution-in-matlab 
# original at: 
# http://www.ozgrid.com/forum/showthread.php?t=108175 
def rand_skew_norm(fAlpha, fLocation, fScale): 
    sigma = fAlpha/np.sqrt(1.0 + fAlpha**2) 

    afRN = np.random.randn(2) 
    u0 = afRN[0] 
    v = afRN[1] 
    u1 = sigma*u0 + np.sqrt(1.0 -sigma**2) * v 

    if u0 >= 0: 
     return u1*fScale + fLocation 
    return (-u1)*fScale + fLocation 

def randn_skew(N, skew=0.0): 
    return [rand_skew_norm(skew, 0, 1) for x in range(N)] 

# lets check they at least visually match the PDF: 
plt.subplots(figsize=(12,4)) 
for alpha_skew in SKEW_PARAMS: 
    p = randn_skew(NUM_SAMPLES, alpha_skew) 
    sns.distplot(p) 

Und dann schrieb eine schnelle Version, die (ohne umfangreiche Tests) korrekt zu sein scheint:

def randn_skew_fast(N, alpha=0.0, loc=0.0, scale=1.0): 
    sigma = alpha/np.sqrt(1.0 + alpha**2) 
    u0 = np.random.randn(N) 
    v = np.random.randn(N) 
    u1 = (sigma*u0 + np.sqrt(1.0 - sigma**2)*v) * scale 
    u1[u0 < 0] *= -1 
    u1 = u1 + loc 
    return u1 

# lets check again 
plt.subplots(figsize=(12,4)) 
for alpha_skew in SKEW_PARAMS: 
    p = randn_skew_fast(NUM_SAMPLES, alpha_skew) 
    sns.distplot(p) 

Answer 2

Angepasst von rsnorm Funktion von fGarch R Paket

def random_snorm(n, mean = 0, sd = 1, xi = 1.5): 
    def random_snorm_aux(n, xi): 
     weight = xi/(xi + 1/xi) 
     z = numpy.random.uniform(-weight,1-weight,n) 
     xi_ = xi**numpy.sign(z) 
     random = -numpy.absolute(numpy.random.normal(0,1,n))/xi_ * numpy.sign(z) 
     m1 = 2/numpy.sqrt(2 * numpy.pi) 
     mu = m1 * (xi - 1/xi) 
     sigma = numpy.sqrt((1 - m1**2) * (xi**2 + 1/xi**2) + 2 * m1**2 - 1) 
     return (random - mu)/sigma 

    return random_snorm_aux(n, xi) * sd + mean