CODEWITHSUNDEEP
pythonscipy
Dori Levanoni
Dori Levanoni

Reputation: 15

scipy skewnorm mean not matching theory?

I'm using scipy skewnorm to create a skewed distribution with a loc and scale.

I adjust the loc and scale passed to scipy.stats.skewnorm based on Adelchi Azzalini's page(Here is link), using the section at the bottom of that page on "mean value" and "delta".

The code I'm using is:

import math
import scipy.stats

skew = -2
mean = 0.05
stdev = 0.05

delta = skew / math.sqrt(1. + math.pow(skew, 2.))
adjMean = mean - stdev * math.sqrt(2. / math.pi) * delta
adjStdev = math.sqrt(math.pow(stdev, 2.) / (1. - 2. * math.pow(delta, 2.) / math.pi))

print 'target mean={:.4f} actual mean={:.4f}'.format(mean, float(scipy.stats.skewnorm.stats(skew, loc=adjMean, scale=adjStdev, moments='mvsk')[0]))
print 'target stdev={:.4f} actual stdev={:.4f}'.format(stdev, math.sqrt(float(scipy.stats.skewnorm.stats(skew, loc=adjMean, scale=adjStdev, moments='mvsk')[1])))

When I run it, though, I'm not getting the mean I expect, while the stdev is what I expect:

target mean=0.0500 actual mean=0.0347
target stdev=0.0500 actual stdev=0.0500

I feel like I'm missing something either about skewnorm or in scipy.stats.skewnorm...

I have numerically integrated the distribution and the mean matches the "actual mean" above.

Upvotes: 0

Views: 686

Answers (1)

Warren Weckesser
Warren Weckesser

Reputation: 114811

You have an algebra mistake. You have

adjMean = mean - stdev * math.sqrt(2. / math.pi) * delta

but on the right side, stdev should be adjStdev.

Here's a modified version of your code:

import math
import scipy.stats

skew = 2.0
mean = 1.5
stdev = 3.0

delta = skew / math.sqrt(1. + math.pow(skew, 2.))
adjStdev = math.sqrt(math.pow(stdev, 2.) / (1. - 2. * math.pow(delta, 2.) / math.pi))
adjMean = mean - adjStdev * math.sqrt(2. / math.pi) * delta

print('target mean={:.4f} actual mean={:.4f}'.format(mean, float(scipy.stats.skewnorm.stats(skew, loc=adjMean, scale=adjStdev, moments='mvsk')[0])))
print('target stdev={:.4f} actual stdev={:.4f}'.format(stdev, math.sqrt(float(scipy.stats.skewnorm.stats(skew, loc=adjMean, scale=adjStdev, moments='mvsk')[1]))))

Here's the output:

target mean=1.5000 actual mean=1.5000
target stdev=3.0000 actual stdev=3.0000

Upvotes: 1

Related Questions