Geek Answers Handbook

Generate N random numbers from an oblique normal distribution using numpy

I need a function in python to return N random numbers from an oblique normal distribution. As a parameter, it is necessary to accept a bias.

eg. my current use

x = numpy.random.randn(1000)

and the ideal function would be for example

x = randn_skew(1000, skew=0.7)

The solution should match: python version 2.7, numpy v.1.9

A similar answer can be found here: change the normal distribution in scipy. However, this creates a PDF not random numbers.

+6
source share
3 answers

I'll start by creating PDF curves for reference:

 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 = (xe) / 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) # nb 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)

Normal Distributed PDF Files

Next, I found a VB implementation of random numbers of a sample from the normal distribution of oblique and converted it to 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)

histograms from oblique normal distributions as generated

And then he wrote a quick version that (without extensive testing) looks right:

 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)

histograms from oblique normal distributions generated by the fast method

+13
source

Adapted from rsnorm function from fGarch R package

 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
+1
source 
 from scipy.stats import skewnorm a=10 data= skewnorm.rvs(a, size=1000)

Here's a parameter you can refer to: https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.skewnorm.html

0
source

More articles:


All Articles