You cannot select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
2080 lines
76 KiB
Python
2080 lines
76 KiB
Python
""" Test functions for stats module
|
|
|
|
"""
|
|
from __future__ import division, print_function, absolute_import
|
|
|
|
import warnings
|
|
import re
|
|
import sys
|
|
|
|
from numpy.testing import (TestCase, run_module_suite, assert_equal,
|
|
assert_array_equal, assert_almost_equal, assert_array_almost_equal,
|
|
assert_allclose, assert_, assert_raises, rand, dec)
|
|
from nose import SkipTest
|
|
|
|
import numpy
|
|
import numpy as np
|
|
from numpy import typecodes, array
|
|
from scipy.lib._version import NumpyVersion
|
|
from scipy import special
|
|
import wafo.stats as stats
|
|
from wafo.stats._distn_infrastructure import argsreduce
|
|
import wafo.stats.distributions
|
|
|
|
from scipy.special import xlogy
|
|
|
|
|
|
# python -OO strips docstrings
|
|
DOCSTRINGS_STRIPPED = sys.flags.optimize > 1
|
|
|
|
|
|
# Generate test cases to test cdf and distribution consistency.
|
|
# Note that this list does not include all distributions.
|
|
dists = ['uniform','norm','lognorm','expon','beta',
|
|
'powerlaw','bradford','burr','fisk','cauchy','halfcauchy',
|
|
'foldcauchy','gamma','gengamma','loggamma',
|
|
'alpha','anglit','arcsine','betaprime',
|
|
'dgamma','exponweib','exponpow','frechet_l','frechet_r',
|
|
'gilbrat','f','ncf','chi2','chi','nakagami','genpareto',
|
|
'genextreme','genhalflogistic','pareto','lomax','halfnorm',
|
|
'halflogistic','fatiguelife','foldnorm','ncx2','t','nct',
|
|
'weibull_min','weibull_max','dweibull','maxwell','rayleigh',
|
|
'genlogistic', 'logistic','gumbel_l','gumbel_r','gompertz',
|
|
'hypsecant', 'laplace', 'reciprocal','triang','tukeylambda',
|
|
'vonmises', 'vonmises_line', 'pearson3']
|
|
|
|
|
|
def _assert_hasattr(a, b, msg=None):
|
|
if msg is None:
|
|
msg = '%s does not have attribute %s' % (a, b)
|
|
assert_(hasattr(a, b), msg=msg)
|
|
|
|
|
|
def test_api_regression():
|
|
# https://github.com/scipy/scipy/issues/3802
|
|
_assert_hasattr(stats.distributions, 'f_gen')
|
|
|
|
|
|
# check function for test generator
|
|
|
|
|
|
def check_distribution(dist, args, alpha):
|
|
D,pval = stats.kstest(dist,'', args=args, N=1000)
|
|
if (pval < alpha):
|
|
D,pval = stats.kstest(dist,'',args=args, N=1000)
|
|
# if (pval < alpha):
|
|
# D,pval = stats.kstest(dist,'',args=args, N=1000)
|
|
assert_(pval > alpha, msg="D = " + str(D) + "; pval = " + str(pval) +
|
|
"; alpha = " + str(alpha) + "\nargs = " + str(args))
|
|
|
|
# nose test generator
|
|
|
|
|
|
def test_all_distributions():
|
|
for dist in dists:
|
|
distfunc = getattr(stats, dist)
|
|
nargs = distfunc.numargs
|
|
alpha = 0.01
|
|
if dist == 'fatiguelife':
|
|
alpha = 0.001
|
|
|
|
if dist == 'frechet':
|
|
args = tuple(2*rand(1))+(0,)+tuple(2*rand(2))
|
|
elif dist == 'triang':
|
|
args = tuple(rand(nargs))
|
|
elif dist == 'reciprocal':
|
|
vals = rand(nargs)
|
|
vals[1] = vals[0] + 1.0
|
|
args = tuple(vals)
|
|
elif dist == 'vonmises':
|
|
yield check_distribution, dist, (10,), alpha
|
|
yield check_distribution, dist, (101,), alpha
|
|
args = tuple(1.0+rand(nargs))
|
|
else:
|
|
args = tuple(1.0+rand(nargs))
|
|
|
|
yield check_distribution, dist, args, alpha
|
|
|
|
|
|
def check_vonmises_pdf_periodic(k,l,s,x):
|
|
vm = stats.vonmises(k,loc=l,scale=s)
|
|
assert_almost_equal(vm.pdf(x),vm.pdf(x % (2*numpy.pi*s)))
|
|
|
|
|
|
def check_vonmises_cdf_periodic(k,l,s,x):
|
|
vm = stats.vonmises(k,loc=l,scale=s)
|
|
assert_almost_equal(vm.cdf(x) % 1,vm.cdf(x % (2*numpy.pi*s)) % 1)
|
|
|
|
|
|
def test_vonmises_pdf_periodic():
|
|
for k in [0.1, 1, 101]:
|
|
for x in [0,1,numpy.pi,10,100]:
|
|
yield check_vonmises_pdf_periodic, k, 0, 1, x
|
|
yield check_vonmises_pdf_periodic, k, 1, 1, x
|
|
yield check_vonmises_pdf_periodic, k, 0, 10, x
|
|
|
|
yield check_vonmises_cdf_periodic, k, 0, 1, x
|
|
yield check_vonmises_cdf_periodic, k, 1, 1, x
|
|
yield check_vonmises_cdf_periodic, k, 0, 10, x
|
|
|
|
|
|
def test_vonmises_line_support():
|
|
assert_equal(stats.vonmises_line.a, -np.pi)
|
|
assert_equal(stats.vonmises_line.b, np.pi)
|
|
|
|
|
|
class TestRandInt(TestCase):
|
|
def test_rvs(self):
|
|
vals = stats.randint.rvs(5,30,size=100)
|
|
assert_(numpy.all(vals < 30) & numpy.all(vals >= 5))
|
|
assert_(len(vals) == 100)
|
|
vals = stats.randint.rvs(5,30,size=(2,50))
|
|
assert_(numpy.shape(vals) == (2,50))
|
|
assert_(vals.dtype.char in typecodes['AllInteger'])
|
|
val = stats.randint.rvs(15,46)
|
|
assert_((val >= 15) & (val < 46))
|
|
assert_(isinstance(val, numpy.ScalarType), msg=repr(type(val)))
|
|
val = stats.randint(15,46).rvs(3)
|
|
assert_(val.dtype.char in typecodes['AllInteger'])
|
|
|
|
def test_pdf(self):
|
|
k = numpy.r_[0:36]
|
|
out = numpy.where((k >= 5) & (k < 30), 1.0/(30-5), 0)
|
|
vals = stats.randint.pmf(k,5,30)
|
|
assert_array_almost_equal(vals,out)
|
|
|
|
def test_cdf(self):
|
|
x = numpy.r_[0:36:100j]
|
|
k = numpy.floor(x)
|
|
out = numpy.select([k >= 30,k >= 5],[1.0,(k-5.0+1)/(30-5.0)],0)
|
|
vals = stats.randint.cdf(x,5,30)
|
|
assert_array_almost_equal(vals, out, decimal=12)
|
|
|
|
|
|
class TestBinom(TestCase):
|
|
def test_rvs(self):
|
|
vals = stats.binom.rvs(10, 0.75, size=(2, 50))
|
|
assert_(numpy.all(vals >= 0) & numpy.all(vals <= 10))
|
|
assert_(numpy.shape(vals) == (2, 50))
|
|
assert_(vals.dtype.char in typecodes['AllInteger'])
|
|
val = stats.binom.rvs(10, 0.75)
|
|
assert_(isinstance(val, int))
|
|
val = stats.binom(10, 0.75).rvs(3)
|
|
assert_(isinstance(val, numpy.ndarray))
|
|
assert_(val.dtype.char in typecodes['AllInteger'])
|
|
|
|
def test_pmf(self):
|
|
# regression test for Ticket #1842
|
|
vals1 = stats.binom.pmf(100, 100,1)
|
|
vals2 = stats.binom.pmf(0, 100,0)
|
|
assert_allclose(vals1, 1.0, rtol=1e-15, atol=0)
|
|
assert_allclose(vals2, 1.0, rtol=1e-15, atol=0)
|
|
|
|
def test_entropy(self):
|
|
# Basic entropy tests.
|
|
b = stats.binom(2, 0.5)
|
|
expected_p = np.array([0.25, 0.5, 0.25])
|
|
expected_h = -sum(xlogy(expected_p, expected_p))
|
|
h = b.entropy()
|
|
assert_allclose(h, expected_h)
|
|
|
|
b = stats.binom(2, 0.0)
|
|
h = b.entropy()
|
|
assert_equal(h, 0.0)
|
|
|
|
b = stats.binom(2, 1.0)
|
|
h = b.entropy()
|
|
assert_equal(h, 0.0)
|
|
|
|
def test_warns_p0(self):
|
|
# no spurious warnigns are generated for p=0; gh-3817
|
|
with warnings.catch_warnings():
|
|
warnings.simplefilter("error", RuntimeWarning)
|
|
assert_equal(stats.binom(n=2, p=0).mean(), 0)
|
|
assert_equal(stats.binom(n=2, p=0).std(), 0)
|
|
|
|
|
|
class TestBernoulli(TestCase):
|
|
def test_rvs(self):
|
|
vals = stats.bernoulli.rvs(0.75, size=(2, 50))
|
|
assert_(numpy.all(vals >= 0) & numpy.all(vals <= 1))
|
|
assert_(numpy.shape(vals) == (2, 50))
|
|
assert_(vals.dtype.char in typecodes['AllInteger'])
|
|
val = stats.bernoulli.rvs(0.75)
|
|
assert_(isinstance(val, int))
|
|
val = stats.bernoulli(0.75).rvs(3)
|
|
assert_(isinstance(val, numpy.ndarray))
|
|
assert_(val.dtype.char in typecodes['AllInteger'])
|
|
|
|
def test_entropy(self):
|
|
# Simple tests of entropy.
|
|
b = stats.bernoulli(0.25)
|
|
expected_h = -0.25*np.log(0.25) - 0.75*np.log(0.75)
|
|
h = b.entropy()
|
|
assert_allclose(h, expected_h)
|
|
|
|
b = stats.bernoulli(0.0)
|
|
h = b.entropy()
|
|
assert_equal(h, 0.0)
|
|
|
|
b = stats.bernoulli(1.0)
|
|
h = b.entropy()
|
|
assert_equal(h, 0.0)
|
|
|
|
|
|
class TestNBinom(TestCase):
|
|
def test_rvs(self):
|
|
vals = stats.nbinom.rvs(10, 0.75, size=(2, 50))
|
|
assert_(numpy.all(vals >= 0))
|
|
assert_(numpy.shape(vals) == (2, 50))
|
|
assert_(vals.dtype.char in typecodes['AllInteger'])
|
|
val = stats.nbinom.rvs(10, 0.75)
|
|
assert_(isinstance(val, int))
|
|
val = stats.nbinom(10, 0.75).rvs(3)
|
|
assert_(isinstance(val, numpy.ndarray))
|
|
assert_(val.dtype.char in typecodes['AllInteger'])
|
|
|
|
def test_pmf(self):
|
|
# regression test for ticket 1779
|
|
assert_allclose(np.exp(stats.nbinom.logpmf(700, 721, 0.52)),
|
|
stats.nbinom.pmf(700, 721, 0.52))
|
|
|
|
|
|
class TestGeom(TestCase):
|
|
def test_rvs(self):
|
|
vals = stats.geom.rvs(0.75, size=(2, 50))
|
|
assert_(numpy.all(vals >= 0))
|
|
assert_(numpy.shape(vals) == (2, 50))
|
|
assert_(vals.dtype.char in typecodes['AllInteger'])
|
|
val = stats.geom.rvs(0.75)
|
|
assert_(isinstance(val, int))
|
|
val = stats.geom(0.75).rvs(3)
|
|
assert_(isinstance(val, numpy.ndarray))
|
|
assert_(val.dtype.char in typecodes['AllInteger'])
|
|
|
|
def test_pmf(self):
|
|
vals = stats.geom.pmf([1,2,3],0.5)
|
|
assert_array_almost_equal(vals,[0.5,0.25,0.125])
|
|
|
|
def test_logpmf(self):
|
|
# regression test for ticket 1793
|
|
vals1 = np.log(stats.geom.pmf([1,2,3], 0.5))
|
|
vals2 = stats.geom.logpmf([1,2,3], 0.5)
|
|
assert_allclose(vals1, vals2, rtol=1e-15, atol=0)
|
|
|
|
def test_cdf_sf(self):
|
|
vals = stats.geom.cdf([1, 2, 3], 0.5)
|
|
vals_sf = stats.geom.sf([1, 2, 3], 0.5)
|
|
expected = array([0.5, 0.75, 0.875])
|
|
assert_array_almost_equal(vals, expected)
|
|
assert_array_almost_equal(vals_sf, 1-expected)
|
|
|
|
def test_logcdf_logsf(self):
|
|
vals = stats.geom.logcdf([1, 2, 3], 0.5)
|
|
vals_sf = stats.geom.logsf([1, 2, 3], 0.5)
|
|
expected = array([0.5, 0.75, 0.875])
|
|
assert_array_almost_equal(vals, np.log(expected))
|
|
assert_array_almost_equal(vals_sf, np.log1p(-expected))
|
|
|
|
def test_ppf(self):
|
|
vals = stats.geom.ppf([0.5, 0.75, 0.875], 0.5)
|
|
expected = array([1.0, 2.0, 3.0])
|
|
assert_array_almost_equal(vals, expected)
|
|
|
|
|
|
class TestTruncnorm(TestCase):
|
|
def test_ppf_ticket1131(self):
|
|
vals = stats.truncnorm.ppf([-0.5,0,1e-4,0.5, 1-1e-4,1,2], -1., 1.,
|
|
loc=[3]*7, scale=2)
|
|
expected = np.array([np.nan, 1, 1.00056419, 3, 4.99943581, 5, np.nan])
|
|
assert_array_almost_equal(vals, expected)
|
|
|
|
def test_isf_ticket1131(self):
|
|
vals = stats.truncnorm.isf([-0.5,0,1e-4,0.5, 1-1e-4,1,2], -1., 1.,
|
|
loc=[3]*7, scale=2)
|
|
expected = np.array([np.nan, 5, 4.99943581, 3, 1.00056419, 1, np.nan])
|
|
assert_array_almost_equal(vals, expected)
|
|
|
|
def test_gh_2477_small_values(self):
|
|
# Check a case that worked in the original issue.
|
|
low, high = -11, -10
|
|
x = stats.truncnorm.rvs(low, high, 0, 1, size=10)
|
|
assert_(low < x.min() < x.max() < high)
|
|
# Check a case that failed in the original issue.
|
|
low, high = 10, 11
|
|
x = stats.truncnorm.rvs(low, high, 0, 1, size=10)
|
|
assert_(low < x.min() < x.max() < high)
|
|
|
|
def test_gh_2477_large_values(self):
|
|
# Check a case that fails because of extreme tailness.
|
|
raise SkipTest('truncnorm rvs is know to fail at extreme tails')
|
|
low, high = 100, 101
|
|
x = stats.truncnorm.rvs(low, high, 0, 1, size=10)
|
|
assert_(low < x.min() < x.max() < high)
|
|
|
|
def test_gh_1489_trac_962_rvs(self):
|
|
# Check the original example.
|
|
low, high = 10, 15
|
|
x = stats.truncnorm.rvs(low, high, 0, 1, size=10)
|
|
assert_(low < x.min() < x.max() < high)
|
|
|
|
|
|
class TestHypergeom(TestCase):
|
|
def test_rvs(self):
|
|
vals = stats.hypergeom.rvs(20, 10, 3, size=(2, 50))
|
|
assert_(numpy.all(vals >= 0) &
|
|
numpy.all(vals <= 3))
|
|
assert_(numpy.shape(vals) == (2, 50))
|
|
assert_(vals.dtype.char in typecodes['AllInteger'])
|
|
val = stats.hypergeom.rvs(20, 3, 10)
|
|
assert_(isinstance(val, int))
|
|
val = stats.hypergeom(20, 3, 10).rvs(3)
|
|
assert_(isinstance(val, numpy.ndarray))
|
|
assert_(val.dtype.char in typecodes['AllInteger'])
|
|
|
|
def test_precision(self):
|
|
# comparison number from mpmath
|
|
M = 2500
|
|
n = 50
|
|
N = 500
|
|
tot = M
|
|
good = n
|
|
hgpmf = stats.hypergeom.pmf(2, tot, good, N)
|
|
assert_almost_equal(hgpmf, 0.0010114963068932233, 11)
|
|
|
|
def test_cdf_above_one(self):
|
|
# for some values of parameters, hypergeom cdf was >1, see gh-2238
|
|
assert_(0 <= stats.hypergeom.cdf(30, 13397950, 4363, 12390) <= 1.0)
|
|
|
|
def test_precision2(self):
|
|
# Test hypergeom precision for large numbers. See #1218.
|
|
# Results compared with those from R.
|
|
oranges = 9.9e4
|
|
pears = 1.1e5
|
|
fruits_eaten = np.array([3, 3.8, 3.9, 4, 4.1, 4.2, 5]) * 1e4
|
|
quantile = 2e4
|
|
res = []
|
|
for eaten in fruits_eaten:
|
|
res.append(stats.hypergeom.sf(quantile, oranges + pears, oranges, eaten))
|
|
expected = np.array([0, 1.904153e-114, 2.752693e-66, 4.931217e-32,
|
|
8.265601e-11, 0.1237904, 1])
|
|
assert_allclose(res, expected, atol=0, rtol=5e-7)
|
|
|
|
# Test with array_like first argument
|
|
quantiles = [1.9e4, 2e4, 2.1e4, 2.15e4]
|
|
res2 = stats.hypergeom.sf(quantiles, oranges + pears, oranges, 4.2e4)
|
|
expected2 = [1, 0.1237904, 6.511452e-34, 3.277667e-69]
|
|
assert_allclose(res2, expected2, atol=0, rtol=5e-7)
|
|
|
|
def test_entropy(self):
|
|
# Simple tests of entropy.
|
|
hg = stats.hypergeom(4, 1, 1)
|
|
h = hg.entropy()
|
|
expected_p = np.array([0.75, 0.25])
|
|
expected_h = -np.sum(xlogy(expected_p, expected_p))
|
|
assert_allclose(h, expected_h)
|
|
|
|
hg = stats.hypergeom(1, 1, 1)
|
|
h = hg.entropy()
|
|
assert_equal(h, 0.0)
|
|
|
|
|
|
class TestLoggamma(TestCase):
|
|
|
|
def test_stats(self):
|
|
# The following precomputed values are from the table in section 2.2
|
|
# of "A Statistical Study of Log-Gamma Distribution", by Ping Shing
|
|
# Chan (thesis, McMaster University, 1993).
|
|
table = np.array([
|
|
# c, mean, var, skew, exc. kurt.
|
|
0.5, -1.9635, 4.9348, -1.5351, 4.0000,
|
|
1.0, -0.5772, 1.6449, -1.1395, 2.4000,
|
|
12.0, 2.4427, 0.0869, -0.2946, 0.1735,
|
|
]).reshape(-1, 5)
|
|
for c, mean, var, skew, kurt in table:
|
|
computed = stats.loggamma.stats(c, moments='msvk')
|
|
assert_array_almost_equal(computed, [mean, var, skew, kurt],
|
|
decimal=4)
|
|
|
|
|
|
class TestLogser(TestCase):
|
|
def test_rvs(self):
|
|
vals = stats.logser.rvs(0.75, size=(2, 50))
|
|
assert_(numpy.all(vals >= 1))
|
|
assert_(numpy.shape(vals) == (2, 50))
|
|
assert_(vals.dtype.char in typecodes['AllInteger'])
|
|
val = stats.logser.rvs(0.75)
|
|
assert_(isinstance(val, int))
|
|
val = stats.logser(0.75).rvs(3)
|
|
assert_(isinstance(val, numpy.ndarray))
|
|
assert_(val.dtype.char in typecodes['AllInteger'])
|
|
|
|
|
|
class TestPareto(TestCase):
|
|
def test_stats(self):
|
|
# Check the stats() method with some simple values. Also check
|
|
# that the calculations do not trigger RuntimeWarnings.
|
|
with warnings.catch_warnings():
|
|
warnings.simplefilter("error", RuntimeWarning)
|
|
|
|
m, v, s, k = stats.pareto.stats(0.5, moments='mvsk')
|
|
assert_equal(m, np.inf)
|
|
assert_equal(v, np.inf)
|
|
assert_equal(s, np.nan)
|
|
assert_equal(k, np.nan)
|
|
|
|
m, v, s, k = stats.pareto.stats(1.0, moments='mvsk')
|
|
assert_equal(m, np.inf)
|
|
assert_equal(v, np.inf)
|
|
assert_equal(s, np.nan)
|
|
assert_equal(k, np.nan)
|
|
|
|
m, v, s, k = stats.pareto.stats(1.5, moments='mvsk')
|
|
assert_equal(m, 3.0)
|
|
assert_equal(v, np.inf)
|
|
assert_equal(s, np.nan)
|
|
assert_equal(k, np.nan)
|
|
|
|
m, v, s, k = stats.pareto.stats(2.0, moments='mvsk')
|
|
assert_equal(m, 2.0)
|
|
assert_equal(v, np.inf)
|
|
assert_equal(s, np.nan)
|
|
assert_equal(k, np.nan)
|
|
|
|
m, v, s, k = stats.pareto.stats(2.5, moments='mvsk')
|
|
assert_allclose(m, 2.5 / 1.5)
|
|
assert_allclose(v, 2.5 / (1.5*1.5*0.5))
|
|
assert_equal(s, np.nan)
|
|
assert_equal(k, np.nan)
|
|
|
|
m, v, s, k = stats.pareto.stats(3.0, moments='mvsk')
|
|
assert_allclose(m, 1.5)
|
|
assert_allclose(v, 0.75)
|
|
assert_equal(s, np.nan)
|
|
assert_equal(k, np.nan)
|
|
|
|
m, v, s, k = stats.pareto.stats(3.5, moments='mvsk')
|
|
assert_allclose(m, 3.5 / 2.5)
|
|
assert_allclose(v, 3.5 / (2.5*2.5*1.5))
|
|
assert_allclose(s, (2*4.5/0.5)*np.sqrt(1.5/3.5))
|
|
assert_equal(k, np.nan)
|
|
|
|
m, v, s, k = stats.pareto.stats(4.0, moments='mvsk')
|
|
assert_allclose(m, 4.0 / 3.0)
|
|
assert_allclose(v, 4.0 / 18.0)
|
|
assert_allclose(s, 2*(1+4.0)/(4.0-3) * np.sqrt((4.0-2)/4.0))
|
|
assert_equal(k, np.nan)
|
|
|
|
m, v, s, k = stats.pareto.stats(4.5, moments='mvsk')
|
|
assert_allclose(m, 4.5 / 3.5)
|
|
assert_allclose(v, 4.5 / (3.5*3.5*2.5))
|
|
assert_allclose(s, (2*5.5/1.5) * np.sqrt(2.5/4.5))
|
|
assert_allclose(k, 6*(4.5**3 + 4.5**2 - 6*4.5 - 2)/(4.5*1.5*0.5))
|
|
|
|
|
|
class TestGenpareto(TestCase):
|
|
def test_ab(self):
|
|
# c >= 0: a, b = [0, inf]
|
|
for c in [1., 0.]:
|
|
c = np.asarray(c)
|
|
stats.genpareto._argcheck(c) # ugh
|
|
assert_equal(stats.genpareto.a, 0.)
|
|
assert_(np.isposinf(stats.genpareto.b))
|
|
|
|
# c < 0: a=0, b=1/|c|
|
|
c = np.asarray(-2.)
|
|
stats.genpareto._argcheck(c)
|
|
assert_allclose([stats.genpareto.a, stats.genpareto.b], [0., 0.5])
|
|
|
|
def test_c0(self):
|
|
# with c=0, genpareto reduces to the exponential distribution
|
|
rv = stats.genpareto(c=0.)
|
|
x = np.linspace(0, 10., 30)
|
|
assert_allclose(rv.pdf(x), stats.expon.pdf(x))
|
|
assert_allclose(rv.cdf(x), stats.expon.cdf(x))
|
|
assert_allclose(rv.sf(x), stats.expon.sf(x))
|
|
|
|
q = np.linspace(0., 1., 10)
|
|
assert_allclose(rv.ppf(q), stats.expon.ppf(q))
|
|
|
|
def test_cm1(self):
|
|
# with c=-1, genpareto reduces to the uniform distr on [0, 1]
|
|
rv = stats.genpareto(c=-1.)
|
|
x = np.linspace(0, 10., 30)
|
|
assert_allclose(rv.pdf(x), stats.uniform.pdf(x))
|
|
assert_allclose(rv.cdf(x), stats.uniform.cdf(x))
|
|
assert_allclose(rv.sf(x), stats.uniform.sf(x))
|
|
|
|
q = np.linspace(0., 1., 10)
|
|
assert_allclose(rv.ppf(q), stats.uniform.ppf(q))
|
|
|
|
# logpdf(1., c=-1) should be zero
|
|
assert_allclose(rv.logpdf(1), 0)
|
|
|
|
def test_x_inf(self):
|
|
# make sure x=inf is handled gracefully
|
|
rv = stats.genpareto(c=0.1)
|
|
assert_allclose([rv.pdf(np.inf), rv.cdf(np.inf)], [0., 1.])
|
|
assert_(np.isneginf(rv.logpdf(np.inf)))
|
|
|
|
rv = stats.genpareto(c=0.)
|
|
assert_allclose([rv.pdf(np.inf), rv.cdf(np.inf)], [0., 1.])
|
|
assert_(np.isneginf(rv.logpdf(np.inf)))
|
|
|
|
rv = stats.genpareto(c=-1.)
|
|
assert_allclose([rv.pdf(np.inf), rv.cdf(np.inf)], [0., 1.])
|
|
assert_(np.isneginf(rv.logpdf(np.inf)))
|
|
|
|
def test_c_continuity(self):
|
|
# pdf is continuous at c=0, -1
|
|
x = np.linspace(0, 10, 30)
|
|
for c in [0, -1]:
|
|
pdf0 = stats.genpareto.pdf(x, c)
|
|
for dc in [1e-14, -1e-14]:
|
|
pdfc = stats.genpareto.pdf(x, c + dc)
|
|
assert_allclose(pdf0, pdfc, atol=1e-12)
|
|
|
|
cdf0 = stats.genpareto.cdf(x, c)
|
|
for dc in [1e-14, 1e-14]:
|
|
cdfc = stats.genpareto.cdf(x, c + dc)
|
|
assert_allclose(cdf0, cdfc, atol=1e-12)
|
|
|
|
def test_c_continuity_ppf(self):
|
|
q = np.r_[np.logspace(1e-12, 0.01, base=0.1),
|
|
np.linspace(0.01, 1, 30, endpoint=False),
|
|
1. - np.logspace(1e-12, 0.01, base=0.1)]
|
|
for c in [0., -1.]:
|
|
ppf0 = stats.genpareto.ppf(q, c)
|
|
for dc in [1e-14, -1e-14]:
|
|
ppfc = stats.genpareto.ppf(q, c + dc)
|
|
assert_allclose(ppf0, ppfc, atol=1e-12)
|
|
|
|
def test_c_continuity_isf(self):
|
|
q = np.r_[np.logspace(1e-12, 0.01, base=0.1),
|
|
np.linspace(0.01, 1, 30, endpoint=False),
|
|
1. - np.logspace(1e-12, 0.01, base=0.1)]
|
|
for c in [0., -1.]:
|
|
isf0 = stats.genpareto.isf(q, c)
|
|
for dc in [1e-14, -1e-14]:
|
|
isfc = stats.genpareto.isf(q, c + dc)
|
|
assert_allclose(isf0, isfc, atol=1e-12)
|
|
|
|
def test_cdf_ppf_roundtrip(self):
|
|
# this should pass with machine precision. hat tip @pbrod
|
|
q = np.r_[np.logspace(1e-12, 0.01, base=0.1),
|
|
np.linspace(0.01, 1, 30, endpoint=False),
|
|
1. - np.logspace(1e-12, 0.01, base=0.1)]
|
|
for c in [1e-8, -1e-18, 1e-15, -1e-15]:
|
|
assert_allclose(stats.genpareto.cdf(stats.genpareto.ppf(q, c), c),
|
|
q, atol=1e-15)
|
|
|
|
|
|
class TestPearson3(TestCase):
|
|
def test_rvs(self):
|
|
vals = stats.pearson3.rvs(0.1, size=(2, 50))
|
|
assert_(numpy.shape(vals) == (2, 50))
|
|
assert_(vals.dtype.char in typecodes['AllFloat'])
|
|
val = stats.pearson3.rvs(0.5)
|
|
assert_(isinstance(val, float))
|
|
val = stats.pearson3(0.5).rvs(3)
|
|
assert_(isinstance(val, numpy.ndarray))
|
|
assert_(val.dtype.char in typecodes['AllFloat'])
|
|
assert_(len(val) == 3)
|
|
|
|
def test_pdf(self):
|
|
vals = stats.pearson3.pdf(2, [0.0, 0.1, 0.2])
|
|
assert_allclose(vals, np.array([0.05399097, 0.05555481, 0.05670246]),
|
|
atol=1e-6)
|
|
vals = stats.pearson3.pdf(-3, 0.1)
|
|
assert_allclose(vals, np.array([0.00313791]), atol=1e-6)
|
|
vals = stats.pearson3.pdf([-3,-2,-1,0,1], 0.1)
|
|
assert_allclose(vals, np.array([0.00313791, 0.05192304, 0.25028092,
|
|
0.39885918, 0.23413173]), atol=1e-6)
|
|
|
|
def test_cdf(self):
|
|
vals = stats.pearson3.cdf(2, [0.0, 0.1, 0.2])
|
|
assert_allclose(vals, np.array([0.97724987, 0.97462004, 0.97213626]),
|
|
atol=1e-6)
|
|
vals = stats.pearson3.cdf(-3, 0.1)
|
|
assert_allclose(vals, [0.00082256], atol=1e-6)
|
|
vals = stats.pearson3.cdf([-3,-2,-1,0,1], 0.1)
|
|
assert_allclose(vals, [8.22563821e-04, 1.99860448e-02, 1.58550710e-01,
|
|
5.06649130e-01, 8.41442111e-01], atol=1e-6)
|
|
|
|
|
|
class TestPoisson(TestCase):
|
|
def test_rvs(self):
|
|
vals = stats.poisson.rvs(0.5, size=(2, 50))
|
|
assert_(numpy.all(vals >= 0))
|
|
assert_(numpy.shape(vals) == (2, 50))
|
|
assert_(vals.dtype.char in typecodes['AllInteger'])
|
|
val = stats.poisson.rvs(0.5)
|
|
assert_(isinstance(val, int))
|
|
val = stats.poisson(0.5).rvs(3)
|
|
assert_(isinstance(val, numpy.ndarray))
|
|
assert_(val.dtype.char in typecodes['AllInteger'])
|
|
|
|
def test_stats(self):
|
|
mu = 16.0
|
|
result = stats.poisson.stats(mu, moments='mvsk')
|
|
assert_allclose(result, [mu, mu, np.sqrt(1.0/mu), 1.0/mu])
|
|
|
|
|
|
class TestZipf(TestCase):
|
|
def test_rvs(self):
|
|
vals = stats.zipf.rvs(1.5, size=(2, 50))
|
|
assert_(numpy.all(vals >= 1))
|
|
assert_(numpy.shape(vals) == (2, 50))
|
|
assert_(vals.dtype.char in typecodes['AllInteger'])
|
|
val = stats.zipf.rvs(1.5)
|
|
assert_(isinstance(val, int))
|
|
val = stats.zipf(1.5).rvs(3)
|
|
assert_(isinstance(val, numpy.ndarray))
|
|
assert_(val.dtype.char in typecodes['AllInteger'])
|
|
|
|
def test_moments(self):
|
|
# n-th moment is finite iff a > n + 1
|
|
m, v = stats.zipf.stats(a=2.8)
|
|
assert_(np.isfinite(m))
|
|
assert_equal(v, np.inf)
|
|
|
|
s, k = stats.zipf.stats(a=4.8, moments='sk')
|
|
assert_(not np.isfinite([s, k]).all())
|
|
|
|
|
|
class TestDLaplace(TestCase):
|
|
def test_rvs(self):
|
|
vals = stats.dlaplace.rvs(1.5, size=(2, 50))
|
|
assert_(numpy.shape(vals) == (2, 50))
|
|
assert_(vals.dtype.char in typecodes['AllInteger'])
|
|
val = stats.dlaplace.rvs(1.5)
|
|
assert_(isinstance(val, int))
|
|
val = stats.dlaplace(1.5).rvs(3)
|
|
assert_(isinstance(val, numpy.ndarray))
|
|
assert_(val.dtype.char in typecodes['AllInteger'])
|
|
assert_(stats.dlaplace.rvs(0.8) is not None)
|
|
|
|
def test_stats(self):
|
|
# compare the explicit formulas w/ direct summation using pmf
|
|
a = 1.
|
|
dl = stats.dlaplace(a)
|
|
m, v, s, k = dl.stats('mvsk')
|
|
|
|
N = 37
|
|
xx = np.arange(-N, N+1)
|
|
pp = dl.pmf(xx)
|
|
m2, m4 = np.sum(pp*xx**2), np.sum(pp*xx**4)
|
|
assert_equal((m, s), (0,0))
|
|
assert_allclose((v, k), (m2, m4/m2**2 - 3.), atol=1e-14, rtol=1e-8)
|
|
|
|
def test_stats2(self):
|
|
a = np.log(2.)
|
|
dl = stats.dlaplace(a)
|
|
m, v, s, k = dl.stats('mvsk')
|
|
assert_equal((m, s), (0.,0.))
|
|
assert_allclose((v, k), (4., 3.25))
|
|
|
|
|
|
class TestInvGamma(TestCase):
|
|
@dec.skipif(NumpyVersion(np.__version__) < '1.7.0',
|
|
"assert_* funcs broken with inf/nan")
|
|
def test_invgamma_inf_gh_1866(self):
|
|
# invgamma's moments are only finite for a>n
|
|
# specific numbers checked w/ boost 1.54
|
|
with warnings.catch_warnings():
|
|
warnings.simplefilter('error', RuntimeWarning)
|
|
mvsk = stats.invgamma.stats(a=19.31, moments='mvsk')
|
|
assert_allclose(mvsk,
|
|
[0.05461496450, 0.0001723162534, 1.020362676, 2.055616582])
|
|
|
|
a = [1.1, 3.1, 5.6]
|
|
mvsk = stats.invgamma.stats(a=a, moments='mvsk')
|
|
expected = ([10., 0.476190476, 0.2173913043], # mmm
|
|
[np.inf, 0.2061430632, 0.01312749422], # vvv
|
|
[np.nan, 41.95235392, 2.919025532], # sss
|
|
[np.nan, np.nan, 24.51923076]) # kkk
|
|
for x, y in zip(mvsk, expected):
|
|
assert_almost_equal(x, y)
|
|
|
|
|
|
class TestF(TestCase):
|
|
def test_f_moments(self):
|
|
# n-th moment of F distributions is only finite for n < dfd / 2
|
|
m, v, s, k = stats.f.stats(11, 6.5, moments='mvsk')
|
|
assert_(np.isfinite(m))
|
|
assert_(np.isfinite(v))
|
|
assert_(np.isfinite(s))
|
|
assert_(not np.isfinite(k))
|
|
|
|
def test_moments_warnings(self):
|
|
# no warnings should be generated for dfd = 2, 4, 6, 8 (div by zero)
|
|
with warnings.catch_warnings():
|
|
warnings.simplefilter('error', RuntimeWarning)
|
|
stats.f.stats(dfn=[11]*4, dfd=[2, 4, 6, 8], moments='mvsk')
|
|
|
|
@dec.knownfailureif(True, 'f stats does not properly broadcast')
|
|
def test_stats_broadcast(self):
|
|
# stats do not fully broadcast just yet
|
|
mv = stats.f.stats(dfn=11, dfd=[11, 12])
|
|
|
|
|
|
def test_rvgeneric_std():
|
|
# Regression test for #1191
|
|
assert_array_almost_equal(stats.t.std([5, 6]), [1.29099445, 1.22474487])
|
|
|
|
|
|
class TestRvDiscrete(TestCase):
|
|
def test_rvs(self):
|
|
states = [-1,0,1,2,3,4]
|
|
probability = [0.0,0.3,0.4,0.0,0.3,0.0]
|
|
samples = 1000
|
|
r = stats.rv_discrete(name='sample',values=(states,probability))
|
|
x = r.rvs(size=samples)
|
|
assert_(isinstance(x, numpy.ndarray))
|
|
|
|
for s,p in zip(states,probability):
|
|
assert_(abs(sum(x == s)/float(samples) - p) < 0.05)
|
|
|
|
x = r.rvs()
|
|
assert_(isinstance(x, int))
|
|
|
|
def test_entropy(self):
|
|
# Basic tests of entropy.
|
|
pvals = np.array([0.25, 0.45, 0.3])
|
|
p = stats.rv_discrete(values=([0, 1, 2], pvals))
|
|
expected_h = -sum(xlogy(pvals, pvals))
|
|
h = p.entropy()
|
|
assert_allclose(h, expected_h)
|
|
|
|
p = stats.rv_discrete(values=([0, 1, 2], [1.0, 0, 0]))
|
|
h = p.entropy()
|
|
assert_equal(h, 0.0)
|
|
|
|
|
|
class TestExpon(TestCase):
|
|
def test_zero(self):
|
|
assert_equal(stats.expon.pdf(0),1)
|
|
|
|
def test_tail(self): # Regression test for ticket 807
|
|
assert_equal(stats.expon.cdf(1e-18), 1e-18)
|
|
assert_equal(stats.expon.isf(stats.expon.sf(40)), 40)
|
|
|
|
|
|
class TestGenExpon(TestCase):
|
|
def test_pdf_unity_area(self):
|
|
from scipy.integrate import simps
|
|
# PDF should integrate to one
|
|
assert_almost_equal(simps(stats.genexpon.pdf(numpy.arange(0,10,0.01),
|
|
0.5, 0.5, 2.0),
|
|
dx=0.01), 1, 1)
|
|
|
|
def test_cdf_bounds(self):
|
|
# CDF should always be positive
|
|
cdf = stats.genexpon.cdf(numpy.arange(0, 10, 0.01), 0.5, 0.5, 2.0)
|
|
assert_(numpy.all((0 <= cdf) & (cdf <= 1)))
|
|
|
|
|
|
class TestExponpow(TestCase):
|
|
def test_tail(self):
|
|
assert_almost_equal(stats.exponpow.cdf(1e-10, 2.), 1e-20)
|
|
assert_almost_equal(stats.exponpow.isf(stats.exponpow.sf(5, .8), .8), 5)
|
|
|
|
|
|
class TestSkellam(TestCase):
|
|
def test_pmf(self):
|
|
# comparison to R
|
|
k = numpy.arange(-10, 15)
|
|
mu1, mu2 = 10, 5
|
|
skpmfR = numpy.array(
|
|
[4.2254582961926893e-005, 1.1404838449648488e-004,
|
|
2.8979625801752660e-004, 6.9177078182101231e-004,
|
|
1.5480716105844708e-003, 3.2412274963433889e-003,
|
|
6.3373707175123292e-003, 1.1552351566696643e-002,
|
|
1.9606152375042644e-002, 3.0947164083410337e-002,
|
|
4.5401737566767360e-002, 6.1894328166820688e-002,
|
|
7.8424609500170578e-002, 9.2418812533573133e-002,
|
|
1.0139793148019728e-001, 1.0371927988298846e-001,
|
|
9.9076583077406091e-002, 8.8546660073089561e-002,
|
|
7.4187842052486810e-002, 5.8392772862200251e-002,
|
|
4.3268692953013159e-002, 3.0248159818374226e-002,
|
|
1.9991434305603021e-002, 1.2516877303301180e-002,
|
|
7.4389876226229707e-003])
|
|
|
|
assert_almost_equal(stats.skellam.pmf(k, mu1, mu2), skpmfR, decimal=15)
|
|
|
|
def test_cdf(self):
|
|
# comparison to R, only 5 decimals
|
|
k = numpy.arange(-10, 15)
|
|
mu1, mu2 = 10, 5
|
|
skcdfR = numpy.array(
|
|
[6.4061475386192104e-005, 1.7810985988267694e-004,
|
|
4.6790611790020336e-004, 1.1596768997212152e-003,
|
|
2.7077485103056847e-003, 5.9489760066490718e-003,
|
|
1.2286346724161398e-002, 2.3838698290858034e-002,
|
|
4.3444850665900668e-002, 7.4392014749310995e-002,
|
|
1.1979375231607835e-001, 1.8168808048289900e-001,
|
|
2.6011268998306952e-001, 3.5253150251664261e-001,
|
|
4.5392943399683988e-001, 5.5764871387982828e-001,
|
|
6.5672529695723436e-001, 7.4527195703032389e-001,
|
|
8.1945979908281064e-001, 8.7785257194501087e-001,
|
|
9.2112126489802404e-001, 9.5136942471639818e-001,
|
|
9.7136085902200120e-001, 9.8387773632530240e-001,
|
|
9.9131672394792536e-001])
|
|
|
|
assert_almost_equal(stats.skellam.cdf(k, mu1, mu2), skcdfR, decimal=5)
|
|
|
|
|
|
class TestLognorm(TestCase):
|
|
def test_pdf(self):
|
|
# Regression test for Ticket #1471: avoid nan with 0/0 situation
|
|
with np.errstate(divide='ignore'):
|
|
pdf = stats.lognorm.pdf([0, 0.5, 1], 1)
|
|
assert_array_almost_equal(pdf, [0.0, 0.62749608, 0.39894228])
|
|
|
|
|
|
class TestBeta(TestCase):
|
|
def test_logpdf(self):
|
|
# Regression test for Ticket #1326: avoid nan with 0*log(0) situation
|
|
logpdf = stats.beta.logpdf(0,1,0.5)
|
|
assert_almost_equal(logpdf, -0.69314718056)
|
|
logpdf = stats.beta.logpdf(0,0.5,1)
|
|
assert_almost_equal(logpdf, np.inf)
|
|
|
|
def test_logpdf_ticket_1866(self):
|
|
alpha, beta = 267, 1472
|
|
x = np.array([0.2, 0.5, 0.6])
|
|
b = stats.beta(alpha, beta)
|
|
assert_allclose(b.logpdf(x).sum(), -1201.699061824062)
|
|
assert_allclose(b.pdf(x), np.exp(b.logpdf(x)))
|
|
|
|
|
|
class TestBetaPrime(TestCase):
|
|
def test_logpdf(self):
|
|
alpha, beta = 267, 1472
|
|
x = np.array([0.2, 0.5, 0.6])
|
|
b = stats.betaprime(alpha, beta)
|
|
assert_(np.isfinite(b.logpdf(x)).all())
|
|
assert_allclose(b.pdf(x), np.exp(b.logpdf(x)))
|
|
|
|
|
|
class TestGamma(TestCase):
|
|
def test_pdf(self):
|
|
# a few test cases to compare with R
|
|
pdf = stats.gamma.pdf(90, 394, scale=1./5)
|
|
assert_almost_equal(pdf, 0.002312341)
|
|
|
|
pdf = stats.gamma.pdf(3, 10, scale=1./5)
|
|
assert_almost_equal(pdf, 0.1620358)
|
|
|
|
def test_logpdf(self):
|
|
# Regression test for Ticket #1326: cornercase avoid nan with 0*log(0)
|
|
# situation
|
|
logpdf = stats.gamma.logpdf(0,1)
|
|
assert_almost_equal(logpdf, 0)
|
|
|
|
|
|
class TestChi2(TestCase):
|
|
# regression tests after precision improvements, ticket:1041, not verified
|
|
def test_precision(self):
|
|
assert_almost_equal(stats.chi2.pdf(1000, 1000), 8.919133934753128e-003, 14)
|
|
assert_almost_equal(stats.chi2.pdf(100, 100), 0.028162503162596778, 14)
|
|
|
|
|
|
class TestArrayArgument(TestCase): # test for ticket:992
|
|
def test_noexception(self):
|
|
rvs = stats.norm.rvs(loc=(np.arange(5)), scale=np.ones(5), size=(10,5))
|
|
assert_equal(rvs.shape, (10,5))
|
|
|
|
|
|
class TestDocstring(TestCase):
|
|
def test_docstrings(self):
|
|
# See ticket #761
|
|
if stats.rayleigh.__doc__ is not None:
|
|
self.assertTrue("rayleigh" in stats.rayleigh.__doc__.lower())
|
|
if stats.bernoulli.__doc__ is not None:
|
|
self.assertTrue("bernoulli" in stats.bernoulli.__doc__.lower())
|
|
|
|
def test_no_name_arg(self):
|
|
# If name is not given, construction shouldn't fail. See #1508.
|
|
stats.rv_continuous()
|
|
stats.rv_discrete()
|
|
|
|
|
|
class TestEntropy(TestCase):
|
|
def test_entropy_positive(self):
|
|
# See ticket #497
|
|
pk = [0.5,0.2,0.3]
|
|
qk = [0.1,0.25,0.65]
|
|
eself = stats.entropy(pk,pk)
|
|
edouble = stats.entropy(pk,qk)
|
|
assert_(0.0 == eself)
|
|
assert_(edouble >= 0.0)
|
|
|
|
def test_entropy_base(self):
|
|
pk = np.ones(16, float)
|
|
S = stats.entropy(pk, base=2.)
|
|
assert_(abs(S - 4.) < 1.e-5)
|
|
|
|
qk = np.ones(16, float)
|
|
qk[:8] = 2.
|
|
S = stats.entropy(pk, qk)
|
|
S2 = stats.entropy(pk, qk, base=2.)
|
|
assert_(abs(S/S2 - np.log(2.)) < 1.e-5)
|
|
|
|
def test_entropy_zero(self):
|
|
# Test for PR-479
|
|
assert_almost_equal(stats.entropy([0, 1, 2]), 0.63651416829481278,
|
|
decimal=12)
|
|
|
|
def test_entropy_2d(self):
|
|
pk = [[0.1, 0.2], [0.6, 0.3], [0.3, 0.5]]
|
|
qk = [[0.2, 0.1], [0.3, 0.6], [0.5, 0.3]]
|
|
assert_array_almost_equal(stats.entropy(pk, qk),
|
|
[0.1933259, 0.18609809])
|
|
|
|
@dec.skipif(NumpyVersion(np.__version__) < '1.7.0',
|
|
"assert_* funcs broken with inf/nan")
|
|
def test_entropy_2d_zero(self):
|
|
pk = [[0.1, 0.2], [0.6, 0.3], [0.3, 0.5]]
|
|
qk = [[0.0, 0.1], [0.3, 0.6], [0.5, 0.3]]
|
|
assert_array_almost_equal(stats.entropy(pk, qk),
|
|
[np.inf, 0.18609809])
|
|
|
|
pk[0][0] = 0.0
|
|
assert_array_almost_equal(stats.entropy(pk, qk),
|
|
[0.17403988, 0.18609809])
|
|
|
|
|
|
def TestArgsreduce():
|
|
a = array([1,3,2,1,2,3,3])
|
|
b,c = argsreduce(a > 1, a, 2)
|
|
|
|
assert_array_equal(b, [3,2,2,3,3])
|
|
assert_array_equal(c, [2,2,2,2,2])
|
|
|
|
b,c = argsreduce(2 > 1, a, 2)
|
|
assert_array_equal(b, a[0])
|
|
assert_array_equal(c, [2])
|
|
|
|
b,c = argsreduce(a > 0, a, 2)
|
|
assert_array_equal(b, a)
|
|
assert_array_equal(c, [2] * numpy.size(a))
|
|
|
|
|
|
class TestFitMethod(object):
|
|
skip = ['ncf']
|
|
|
|
@dec.slow
|
|
def test_fit(self):
|
|
def check(func, dist, args, alpha):
|
|
if dist in self.skip:
|
|
raise SkipTest("%s fit known to fail" % dist)
|
|
distfunc = getattr(stats, dist)
|
|
with np.errstate(all='ignore'):
|
|
res = distfunc.rvs(*args, **{'size':200})
|
|
vals = distfunc.fit(res)
|
|
vals2 = distfunc.fit(res, optimizer='powell')
|
|
# Only check the length of the return
|
|
# FIXME: should check the actual results to see if we are 'close'
|
|
# to what was created --- but what is 'close' enough
|
|
if dist == 'frechet':
|
|
assert_(len(vals) == len(args))
|
|
assert_(len(vals2) == len(args))
|
|
else:
|
|
assert_(len(vals) == 2+len(args))
|
|
assert_(len(vals2) == 2+len(args))
|
|
|
|
for func, dist, args, alpha in test_all_distributions():
|
|
yield check, func, dist, args, alpha
|
|
|
|
@dec.slow
|
|
def test_fix_fit(self):
|
|
def check(func, dist, args, alpha):
|
|
# Not sure why 'ncf', and 'beta' are failing
|
|
# frechet has different len(args) than distfunc.numargs
|
|
if dist in self.skip + ['frechet']:
|
|
raise SkipTest("%s fit known to fail" % dist)
|
|
distfunc = getattr(stats, dist)
|
|
with np.errstate(all='ignore'):
|
|
res = distfunc.rvs(*args, **{'size':200})
|
|
vals = distfunc.fit(res,floc=0)
|
|
vals2 = distfunc.fit(res,fscale=1)
|
|
assert_(len(vals) == 2+len(args))
|
|
assert_(vals[-2] == 0)
|
|
assert_(vals2[-1] == 1)
|
|
assert_(len(vals2) == 2+len(args))
|
|
if len(args) > 0:
|
|
vals3 = distfunc.fit(res, f0=args[0])
|
|
assert_(len(vals3) == 2+len(args))
|
|
assert_(vals3[0] == args[0])
|
|
if len(args) > 1:
|
|
vals4 = distfunc.fit(res, f1=args[1])
|
|
assert_(len(vals4) == 2+len(args))
|
|
assert_(vals4[1] == args[1])
|
|
if len(args) > 2:
|
|
vals5 = distfunc.fit(res, f2=args[2])
|
|
assert_(len(vals5) == 2+len(args))
|
|
assert_(vals5[2] == args[2])
|
|
|
|
for func, dist, args, alpha in test_all_distributions():
|
|
yield check, func, dist, args, alpha
|
|
|
|
def test_fix_fit_2args_lognorm(self):
|
|
# Regression test for #1551.
|
|
np.random.seed(12345)
|
|
with np.errstate(all='ignore'):
|
|
x = stats.lognorm.rvs(0.25, 0., 20.0, size=20)
|
|
assert_allclose(np.array(stats.lognorm.fit(x, floc=0, fscale=20)),
|
|
[0.25888672, 0, 20], atol=1e-5)
|
|
|
|
def test_fix_fit_norm(self):
|
|
x = np.arange(1, 6)
|
|
|
|
loc, scale = stats.norm.fit(x)
|
|
assert_almost_equal(loc, 3)
|
|
assert_almost_equal(scale, np.sqrt(2))
|
|
|
|
loc, scale = stats.norm.fit(x, floc=2)
|
|
assert_equal(loc, 2)
|
|
assert_equal(scale, np.sqrt(3))
|
|
|
|
loc, scale = stats.norm.fit(x, fscale=2)
|
|
assert_almost_equal(loc, 3)
|
|
assert_equal(scale, 2)
|
|
|
|
def test_fix_fit_gamma(self):
|
|
x = np.arange(1, 6)
|
|
meanlog = np.log(x).mean()
|
|
|
|
# A basic test of gamma.fit with floc=0.
|
|
floc = 0
|
|
a, loc, scale = stats.gamma.fit(x, floc=floc)
|
|
s = np.log(x.mean()) - meanlog
|
|
assert_almost_equal(np.log(a) - special.digamma(a), s, decimal=5)
|
|
assert_equal(loc, floc)
|
|
assert_almost_equal(scale, x.mean()/a, decimal=8)
|
|
|
|
# Regression tests for gh-2514.
|
|
# The problem was that if `floc=0` was given, any other fixed
|
|
# parameters were ignored.
|
|
f0 = 1
|
|
floc = 0
|
|
a, loc, scale = stats.gamma.fit(x, f0=f0, floc=floc)
|
|
assert_equal(a, f0)
|
|
assert_equal(loc, floc)
|
|
assert_almost_equal(scale, x.mean()/a, decimal=8)
|
|
|
|
f0 = 2
|
|
floc = 0
|
|
a, loc, scale = stats.gamma.fit(x, f0=f0, floc=floc)
|
|
assert_equal(a, f0)
|
|
assert_equal(loc, floc)
|
|
assert_almost_equal(scale, x.mean()/a, decimal=8)
|
|
|
|
# loc and scale fixed.
|
|
floc = 0
|
|
fscale = 2
|
|
a, loc, scale = stats.gamma.fit(x, floc=floc, fscale=fscale)
|
|
assert_equal(loc, floc)
|
|
assert_equal(scale, fscale)
|
|
c = meanlog - np.log(fscale)
|
|
assert_almost_equal(special.digamma(a), c)
|
|
|
|
def test_fix_fit_beta(self):
|
|
# Test beta.fit when both floc and fscale are given.
|
|
|
|
def mlefunc(a, b, x):
|
|
# Zeros of this function are critical points of
|
|
# the maximum likelihood function.
|
|
n = len(x)
|
|
s1 = np.log(x).sum()
|
|
s2 = np.log(1-x).sum()
|
|
psiab = special.psi(a + b)
|
|
func = [s1 - n * (-psiab + special.psi(a)),
|
|
s2 - n * (-psiab + special.psi(b))]
|
|
return func
|
|
|
|
# Basic test with floc and fscale given.
|
|
x = np.array([0.125, 0.25, 0.5])
|
|
a, b, loc, scale = stats.beta.fit(x, floc=0, fscale=1)
|
|
assert_equal(loc, 0)
|
|
assert_equal(scale, 1)
|
|
assert_allclose(mlefunc(a, b, x), [0,0], atol=1e-6)
|
|
|
|
# Basic test with f0, floc and fscale given.
|
|
# This is also a regression test for gh-2514.
|
|
x = np.array([0.125, 0.25, 0.5])
|
|
a, b, loc, scale = stats.beta.fit(x, f0=2, floc=0, fscale=1)
|
|
assert_equal(a, 2)
|
|
assert_equal(loc, 0)
|
|
assert_equal(scale, 1)
|
|
da, db = mlefunc(a, b, x)
|
|
assert_allclose(db, 0, atol=1e-5)
|
|
|
|
# Same floc and fscale values as above, but reverse the data
|
|
# and fix b (f1).
|
|
x2 = 1 - x
|
|
a2, b2, loc2, scale2 = stats.beta.fit(x2, f1=2, floc=0, fscale=1)
|
|
assert_equal(b2, 2)
|
|
assert_equal(loc2, 0)
|
|
assert_equal(scale2, 1)
|
|
da, db = mlefunc(a2, b2, x2)
|
|
assert_allclose(da, 0, atol=1e-5)
|
|
# a2 of this test should equal b from above.
|
|
assert_almost_equal(a2, b)
|
|
|
|
# Check for detection of data out of bounds when floc and fscale
|
|
# are given.
|
|
assert_raises(ValueError, stats.beta.fit, x, floc=0.5, fscale=1)
|
|
y = np.array([0, .5, 1])
|
|
assert_raises(ValueError, stats.beta.fit, y, floc=0, fscale=1)
|
|
assert_raises(ValueError, stats.beta.fit, y, floc=0, fscale=1, f0=2)
|
|
assert_raises(ValueError, stats.beta.fit, y, floc=0, fscale=1, f1=2)
|
|
|
|
# Check that attempting to fix all the parameters raises a ValueError.
|
|
assert_raises(ValueError, stats.beta.fit, y, f0=0, f1=1,
|
|
floc=2, fscale=3)
|
|
|
|
|
|
class TestFrozen(TestCase):
|
|
# Test that a frozen distribution gives the same results as the original object.
|
|
#
|
|
# Only tested for the normal distribution (with loc and scale specified)
|
|
# and for the gamma distribution (with a shape parameter specified).
|
|
def test_norm(self):
|
|
dist = stats.norm
|
|
frozen = stats.norm(loc=10.0, scale=3.0)
|
|
|
|
result_f = frozen.pdf(20.0)
|
|
result = dist.pdf(20.0, loc=10.0, scale=3.0)
|
|
assert_equal(result_f, result)
|
|
|
|
result_f = frozen.cdf(20.0)
|
|
result = dist.cdf(20.0, loc=10.0, scale=3.0)
|
|
assert_equal(result_f, result)
|
|
|
|
result_f = frozen.ppf(0.25)
|
|
result = dist.ppf(0.25, loc=10.0, scale=3.0)
|
|
assert_equal(result_f, result)
|
|
|
|
result_f = frozen.isf(0.25)
|
|
result = dist.isf(0.25, loc=10.0, scale=3.0)
|
|
assert_equal(result_f, result)
|
|
|
|
result_f = frozen.sf(10.0)
|
|
result = dist.sf(10.0, loc=10.0, scale=3.0)
|
|
assert_equal(result_f, result)
|
|
|
|
result_f = frozen.median()
|
|
result = dist.median(loc=10.0, scale=3.0)
|
|
assert_equal(result_f, result)
|
|
|
|
result_f = frozen.mean()
|
|
result = dist.mean(loc=10.0, scale=3.0)
|
|
assert_equal(result_f, result)
|
|
|
|
result_f = frozen.var()
|
|
result = dist.var(loc=10.0, scale=3.0)
|
|
assert_equal(result_f, result)
|
|
|
|
result_f = frozen.std()
|
|
result = dist.std(loc=10.0, scale=3.0)
|
|
assert_equal(result_f, result)
|
|
|
|
result_f = frozen.entropy()
|
|
result = dist.entropy(loc=10.0, scale=3.0)
|
|
assert_equal(result_f, result)
|
|
|
|
result_f = frozen.moment(2)
|
|
result = dist.moment(2,loc=10.0, scale=3.0)
|
|
assert_equal(result_f, result)
|
|
|
|
def test_gamma(self):
|
|
a = 2.0
|
|
dist = stats.gamma
|
|
frozen = stats.gamma(a)
|
|
|
|
result_f = frozen.pdf(20.0)
|
|
result = dist.pdf(20.0, a)
|
|
assert_equal(result_f, result)
|
|
|
|
result_f = frozen.cdf(20.0)
|
|
result = dist.cdf(20.0, a)
|
|
assert_equal(result_f, result)
|
|
|
|
result_f = frozen.ppf(0.25)
|
|
result = dist.ppf(0.25, a)
|
|
assert_equal(result_f, result)
|
|
|
|
result_f = frozen.isf(0.25)
|
|
result = dist.isf(0.25, a)
|
|
assert_equal(result_f, result)
|
|
|
|
result_f = frozen.sf(10.0)
|
|
result = dist.sf(10.0, a)
|
|
assert_equal(result_f, result)
|
|
|
|
result_f = frozen.median()
|
|
result = dist.median(a)
|
|
assert_equal(result_f, result)
|
|
|
|
result_f = frozen.mean()
|
|
result = dist.mean(a)
|
|
assert_equal(result_f, result)
|
|
|
|
result_f = frozen.var()
|
|
result = dist.var(a)
|
|
assert_equal(result_f, result)
|
|
|
|
result_f = frozen.std()
|
|
result = dist.std(a)
|
|
assert_equal(result_f, result)
|
|
|
|
result_f = frozen.entropy()
|
|
result = dist.entropy(a)
|
|
assert_equal(result_f, result)
|
|
|
|
result_f = frozen.moment(2)
|
|
result = dist.moment(2, a)
|
|
assert_equal(result_f, result)
|
|
|
|
def test_regression_ticket_1293(self):
|
|
# Create a frozen distribution.
|
|
frozen = stats.lognorm(1)
|
|
# Call one of its methods that does not take any keyword arguments.
|
|
m1 = frozen.moment(2)
|
|
# Now call a method that takes a keyword argument.
|
|
frozen.stats(moments='mvsk')
|
|
# Call moment(2) again.
|
|
# After calling stats(), the following was raising an exception.
|
|
# So this test passes if the following does not raise an exception.
|
|
m2 = frozen.moment(2)
|
|
# The following should also be true, of course. But it is not
|
|
# the focus of this test.
|
|
assert_equal(m1, m2)
|
|
|
|
def test_ab(self):
|
|
# test that the support of a frozen distribution
|
|
# (i) remains frozen even if it changes for the original one
|
|
# (ii) is actually correct if the shape parameters are such that
|
|
# the values of [a, b] are not the default [0, inf]
|
|
# take a genpareto as an example where the support
|
|
# depends on the value of the shape parameter:
|
|
# for c > 0: a, b = 0, inf
|
|
# for c < 0: a, b = 0, -1/c
|
|
rv = stats.genpareto(c=-0.1)
|
|
a, b = rv.dist.a, rv.dist.b
|
|
assert_equal([a, b], [0., 10.])
|
|
|
|
stats.genpareto.pdf(0, c=0.1) # this changes genpareto.b
|
|
assert_equal([rv.dist.a, rv.dist.b], [a, b])
|
|
|
|
rv1 = stats.genpareto(c=0.1)
|
|
assert_(rv1.dist is not rv.dist)
|
|
|
|
def test_rv_frozen_in_namespace(self):
|
|
# Regression test for gh-3522
|
|
assert_(hasattr(stats.distributions, 'rv_frozen'))
|
|
|
|
|
|
class TestExpect(TestCase):
|
|
# Test for expect method.
|
|
#
|
|
# Uses normal distribution and beta distribution for finite bounds, and
|
|
# hypergeom for discrete distribution with finite support
|
|
def test_norm(self):
|
|
v = stats.norm.expect(lambda x: (x-5)*(x-5), loc=5, scale=2)
|
|
assert_almost_equal(v, 4, decimal=14)
|
|
|
|
m = stats.norm.expect(lambda x: (x), loc=5, scale=2)
|
|
assert_almost_equal(m, 5, decimal=14)
|
|
|
|
lb = stats.norm.ppf(0.05, loc=5, scale=2)
|
|
ub = stats.norm.ppf(0.95, loc=5, scale=2)
|
|
prob90 = stats.norm.expect(lambda x: 1, loc=5, scale=2, lb=lb, ub=ub)
|
|
assert_almost_equal(prob90, 0.9, decimal=14)
|
|
|
|
prob90c = stats.norm.expect(lambda x: 1, loc=5, scale=2, lb=lb, ub=ub,
|
|
conditional=True)
|
|
assert_almost_equal(prob90c, 1., decimal=14)
|
|
|
|
def test_beta(self):
|
|
# case with finite support interval
|
|
v = stats.beta.expect(lambda x: (x-19/3.)*(x-19/3.), args=(10,5),
|
|
loc=5, scale=2)
|
|
assert_almost_equal(v, 1./18., decimal=13)
|
|
|
|
m = stats.beta.expect(lambda x: x, args=(10,5), loc=5., scale=2.)
|
|
assert_almost_equal(m, 19/3., decimal=13)
|
|
|
|
ub = stats.beta.ppf(0.95, 10, 10, loc=5, scale=2)
|
|
lb = stats.beta.ppf(0.05, 10, 10, loc=5, scale=2)
|
|
prob90 = stats.beta.expect(lambda x: 1., args=(10,10), loc=5.,
|
|
scale=2.,lb=lb, ub=ub, conditional=False)
|
|
assert_almost_equal(prob90, 0.9, decimal=13)
|
|
|
|
prob90c = stats.beta.expect(lambda x: 1, args=(10,10), loc=5,
|
|
scale=2, lb=lb, ub=ub, conditional=True)
|
|
assert_almost_equal(prob90c, 1., decimal=13)
|
|
|
|
def test_hypergeom(self):
|
|
# test case with finite bounds
|
|
|
|
# without specifying bounds
|
|
m_true, v_true = stats.hypergeom.stats(20, 10, 8, loc=5.)
|
|
m = stats.hypergeom.expect(lambda x: x, args=(20, 10, 8), loc=5.)
|
|
assert_almost_equal(m, m_true, decimal=13)
|
|
|
|
v = stats.hypergeom.expect(lambda x: (x-9.)**2, args=(20, 10, 8),
|
|
loc=5.)
|
|
assert_almost_equal(v, v_true, decimal=14)
|
|
|
|
# with bounds, bounds equal to shifted support
|
|
v_bounds = stats.hypergeom.expect(lambda x: (x-9.)**2, args=(20, 10, 8),
|
|
loc=5., lb=5, ub=13)
|
|
assert_almost_equal(v_bounds, v_true, decimal=14)
|
|
|
|
# drop boundary points
|
|
prob_true = 1-stats.hypergeom.pmf([5, 13], 20, 10, 8, loc=5).sum()
|
|
prob_bounds = stats.hypergeom.expect(lambda x: 1, args=(20, 10, 8),
|
|
loc=5., lb=6, ub=12)
|
|
assert_almost_equal(prob_bounds, prob_true, decimal=13)
|
|
|
|
# conditional
|
|
prob_bc = stats.hypergeom.expect(lambda x: 1, args=(20, 10, 8), loc=5.,
|
|
lb=6, ub=12, conditional=True)
|
|
assert_almost_equal(prob_bc, 1, decimal=14)
|
|
|
|
# check simple integral
|
|
prob_b = stats.hypergeom.expect(lambda x: 1, args=(20, 10, 8),
|
|
lb=0, ub=8)
|
|
assert_almost_equal(prob_b, 1, decimal=13)
|
|
|
|
def test_poisson(self):
|
|
# poisson, use lower bound only
|
|
prob_bounds = stats.poisson.expect(lambda x: 1, args=(2,), lb=3,
|
|
conditional=False)
|
|
prob_b_true = 1-stats.poisson.cdf(2,2)
|
|
assert_almost_equal(prob_bounds, prob_b_true, decimal=14)
|
|
|
|
prob_lb = stats.poisson.expect(lambda x: 1, args=(2,), lb=2,
|
|
conditional=True)
|
|
assert_almost_equal(prob_lb, 1, decimal=14)
|
|
|
|
def test_genhalflogistic(self):
|
|
# genhalflogistic, changes upper bound of support in _argcheck
|
|
# regression test for gh-2622
|
|
halflog = stats.genhalflogistic
|
|
# check consistency when calling expect twice with the same input
|
|
res1 = halflog.expect(args=(1.5,))
|
|
halflog.expect(args=(0.5,))
|
|
res2 = halflog.expect(args=(1.5,))
|
|
assert_almost_equal(res1, res2, decimal=14)
|
|
|
|
def test_rice_overflow(self):
|
|
# rice.pdf(999, 0.74) was inf since special.i0 silentyly overflows
|
|
# check that using i0e fixes it
|
|
assert_(np.isfinite(stats.rice.pdf(999, 0.74)))
|
|
|
|
assert_(np.isfinite(stats.rice.expect(lambda x: 1, args=(0.74,))))
|
|
assert_(np.isfinite(stats.rice.expect(lambda x: 2, args=(0.74,))))
|
|
assert_(np.isfinite(stats.rice.expect(lambda x: 3, args=(0.74,))))
|
|
|
|
|
|
class TestNct(TestCase):
|
|
def test_nc_parameter(self):
|
|
# Parameter values c<=0 were not enabled (gh-2402).
|
|
# For negative values c and for c=0 results of rv.cdf(0) below were nan
|
|
rv = stats.nct(5, 0)
|
|
assert_equal(rv.cdf(0), 0.5)
|
|
rv = stats.nct(5, -1)
|
|
assert_almost_equal(rv.cdf(0), 0.841344746069, decimal=10)
|
|
|
|
def test_broadcasting(self):
|
|
res = stats.nct.pdf(5, np.arange(4,7)[:,None], np.linspace(0.1, 1, 4))
|
|
expected = array([[0.00321886, 0.00557466, 0.00918418, 0.01442997],
|
|
[0.00217142, 0.00395366, 0.00683888, 0.01126276],
|
|
[0.00153078, 0.00291093, 0.00525206, 0.00900815]])
|
|
assert_allclose(res, expected, rtol=1e-5)
|
|
|
|
def text_variance_gh_issue_2401(self):
|
|
# Computation of the variance of a non-central t-distribution resulted
|
|
# in a TypeError: ufunc 'isinf' not supported for the input types,
|
|
# and the inputs could not be safely coerced to any supported types
|
|
# according to the casting rule 'safe'
|
|
rv = stats.nct(4, 0)
|
|
assert_equal(rv.var(), 2.0)
|
|
|
|
def test_nct_inf_moments(self):
|
|
# n-th moment of nct only exists for df > n
|
|
m, v, s, k = stats.nct.stats(df=1.9, nc=0.3, moments='mvsk')
|
|
assert_(np.isfinite(m))
|
|
assert_equal([v, s, k], [np.inf, np.nan, np.nan])
|
|
|
|
m, v, s, k = stats.nct.stats(df=3.1, nc=0.3, moments='mvsk')
|
|
assert_(np.isfinite([m, v, s]).all())
|
|
assert_equal(k, np.nan)
|
|
|
|
|
|
class TestRice(TestCase):
|
|
def test_rice_zero_b(self):
|
|
# rice distribution should work with b=0, cf gh-2164
|
|
x = [0.2, 1., 5.]
|
|
assert_(np.isfinite(stats.rice.pdf(x, b=0.)).all())
|
|
assert_(np.isfinite(stats.rice.logpdf(x, b=0.)).all())
|
|
assert_(np.isfinite(stats.rice.cdf(x, b=0.)).all())
|
|
assert_(np.isfinite(stats.rice.logcdf(x, b=0.)).all())
|
|
|
|
q = [0.1, 0.1, 0.5, 0.9]
|
|
assert_(np.isfinite(stats.rice.ppf(q, b=0.)).all())
|
|
|
|
mvsk = stats.rice.stats(0, moments='mvsk')
|
|
assert_(np.isfinite(mvsk).all())
|
|
|
|
# furthermore, pdf is continuous as b\to 0
|
|
# rice.pdf(x, b\to 0) = x exp(-x^2/2) + O(b^2)
|
|
# see e.g. Abramovich & Stegun 9.6.7 & 9.6.10
|
|
b = 1e-8
|
|
assert_allclose(stats.rice.pdf(x, 0), stats.rice.pdf(x, b),
|
|
atol=b, rtol=0)
|
|
|
|
def test_rice_rvs(self):
|
|
rvs = stats.rice.rvs
|
|
assert_equal(rvs(b=3.).size, 1)
|
|
assert_equal(rvs(b=3., size=(3, 5)).shape, (3, 5))
|
|
|
|
|
|
class TestErlang(TestCase):
|
|
def test_erlang_runtimewarning(self):
|
|
# erlang should generate a RuntimeWarning if a non-integer
|
|
# shape parameter is used.
|
|
with warnings.catch_warnings():
|
|
warnings.simplefilter("error", RuntimeWarning)
|
|
|
|
# The non-integer shape parameter 1.3 should trigger a RuntimeWarning
|
|
assert_raises(RuntimeWarning,
|
|
stats.erlang.rvs, 1.3, loc=0, scale=1, size=4)
|
|
|
|
# Calling the fit method with `f0` set to an integer should
|
|
# *not* trigger a RuntimeWarning. It should return the same
|
|
# values as gamma.fit(...).
|
|
data = [0.5, 1.0, 2.0, 4.0]
|
|
result_erlang = stats.erlang.fit(data, f0=1)
|
|
result_gamma = stats.gamma.fit(data, f0=1)
|
|
assert_allclose(result_erlang, result_gamma, rtol=1e-3)
|
|
|
|
|
|
class TestExponWeib(TestCase):
|
|
|
|
def test_pdf_logpdf(self):
|
|
# Regression test for gh-3508.
|
|
x = 0.1
|
|
a = 1.0
|
|
c = 100.0
|
|
p = stats.exponweib.pdf(x, a, c)
|
|
logp = stats.exponweib.logpdf(x, a, c)
|
|
# Expected values were computed with mpmath.
|
|
assert_allclose([p, logp],
|
|
[1.0000000000000054e-97, -223.35075402042244])
|
|
|
|
def test_a_is_1(self):
|
|
# For issue gh-3508.
|
|
# Check that when a=1, the pdf and logpdf methods of exponweib are the
|
|
# same as those of weibull_min.
|
|
x = np.logspace(-4, -1, 4)
|
|
a = 1
|
|
c = 100
|
|
|
|
p = stats.exponweib.pdf(x, a, c)
|
|
expected = stats.weibull_min.pdf(x, c)
|
|
assert_allclose(p, expected)
|
|
|
|
logp = stats.exponweib.logpdf(x, a, c)
|
|
expected = stats.weibull_min.logpdf(x, c)
|
|
assert_allclose(logp, expected)
|
|
|
|
def test_a_is_1_c_is_1(self):
|
|
# When a = 1 and c = 1, the distribution is exponential.
|
|
x = np.logspace(-8, 1, 10)
|
|
a = 1
|
|
c = 1
|
|
|
|
p = stats.exponweib.pdf(x, a, c)
|
|
expected = stats.expon.pdf(x)
|
|
assert_allclose(p, expected)
|
|
|
|
logp = stats.exponweib.logpdf(x, a, c)
|
|
expected = stats.expon.logpdf(x)
|
|
assert_allclose(logp, expected)
|
|
|
|
|
|
class TestRdist(TestCase):
|
|
@dec.slow
|
|
def test_rdist_cdf_gh1285(self):
|
|
# check workaround in rdist._cdf for issue gh-1285.
|
|
distfn = stats.rdist
|
|
values = [0.001, 0.5, 0.999]
|
|
assert_almost_equal(distfn.cdf(distfn.ppf(values, 541.0), 541.0),
|
|
values, decimal=5)
|
|
|
|
|
|
def test_540_567():
|
|
# test for nan returned in tickets 540, 567
|
|
assert_almost_equal(stats.norm.cdf(-1.7624320982),0.03899815971089126,
|
|
decimal=10, err_msg='test_540_567')
|
|
assert_almost_equal(stats.norm.cdf(-1.7624320983),0.038998159702449846,
|
|
decimal=10, err_msg='test_540_567')
|
|
assert_almost_equal(stats.norm.cdf(1.38629436112, loc=0.950273420309,
|
|
scale=0.204423758009),0.98353464004309321,
|
|
decimal=10, err_msg='test_540_567')
|
|
|
|
|
|
def test_regression_ticket_1316():
|
|
# The following was raising an exception, because _construct_default_doc()
|
|
# did not handle the default keyword extradoc=None. See ticket #1316.
|
|
g = stats._continuous_distns.gamma_gen(name='gamma')
|
|
|
|
|
|
def test_regression_ticket_1326():
|
|
# adjust to avoid nan with 0*log(0)
|
|
assert_almost_equal(stats.chi2.pdf(0.0, 2), 0.5, 14)
|
|
|
|
|
|
def test_regression_tukey_lambda():
|
|
# Make sure that Tukey-Lambda distribution correctly handles non-positive lambdas.
|
|
x = np.linspace(-5.0, 5.0, 101)
|
|
|
|
olderr = np.seterr(divide='ignore')
|
|
try:
|
|
for lam in [0.0, -1.0, -2.0, np.array([[-1.0], [0.0], [-2.0]])]:
|
|
p = stats.tukeylambda.pdf(x, lam)
|
|
assert_((p != 0.0).all())
|
|
assert_(~np.isnan(p).all())
|
|
|
|
lam = np.array([[-1.0], [0.0], [2.0]])
|
|
p = stats.tukeylambda.pdf(x, lam)
|
|
finally:
|
|
np.seterr(**olderr)
|
|
|
|
assert_(~np.isnan(p).all())
|
|
assert_((p[0] != 0.0).all())
|
|
assert_((p[1] != 0.0).all())
|
|
assert_((p[2] != 0.0).any())
|
|
assert_((p[2] == 0.0).any())
|
|
|
|
|
|
@dec.skipif(DOCSTRINGS_STRIPPED)
|
|
def test_regression_ticket_1421():
|
|
assert_('pdf(x, mu, loc=0, scale=1)' not in stats.poisson.__doc__)
|
|
assert_('pmf(x,' in stats.poisson.__doc__)
|
|
|
|
|
|
def test_nan_arguments_gh_issue_1362():
|
|
with np.errstate(invalid='ignore'):
|
|
assert_(np.isnan(stats.t.logcdf(1, np.nan)))
|
|
assert_(np.isnan(stats.t.cdf(1, np.nan)))
|
|
assert_(np.isnan(stats.t.logsf(1, np.nan)))
|
|
assert_(np.isnan(stats.t.sf(1, np.nan)))
|
|
assert_(np.isnan(stats.t.pdf(1, np.nan)))
|
|
assert_(np.isnan(stats.t.logpdf(1, np.nan)))
|
|
assert_(np.isnan(stats.t.ppf(1, np.nan)))
|
|
assert_(np.isnan(stats.t.isf(1, np.nan)))
|
|
|
|
assert_(np.isnan(stats.bernoulli.logcdf(np.nan, 0.5)))
|
|
assert_(np.isnan(stats.bernoulli.cdf(np.nan, 0.5)))
|
|
assert_(np.isnan(stats.bernoulli.logsf(np.nan, 0.5)))
|
|
assert_(np.isnan(stats.bernoulli.sf(np.nan, 0.5)))
|
|
assert_(np.isnan(stats.bernoulli.pmf(np.nan, 0.5)))
|
|
assert_(np.isnan(stats.bernoulli.logpmf(np.nan, 0.5)))
|
|
assert_(np.isnan(stats.bernoulli.ppf(np.nan, 0.5)))
|
|
assert_(np.isnan(stats.bernoulli.isf(np.nan, 0.5)))
|
|
|
|
|
|
def test_frozen_fit_ticket_1536():
|
|
np.random.seed(5678)
|
|
true = np.array([0.25, 0., 0.5])
|
|
x = stats.lognorm.rvs(true[0], true[1], true[2], size=100)
|
|
|
|
olderr = np.seterr(divide='ignore')
|
|
try:
|
|
params = np.array(stats.lognorm.fit(x, floc=0.))
|
|
finally:
|
|
np.seterr(**olderr)
|
|
|
|
assert_almost_equal(params, true, decimal=2)
|
|
|
|
params = np.array(stats.lognorm.fit(x, fscale=0.5, loc=0))
|
|
assert_almost_equal(params, true, decimal=2)
|
|
|
|
params = np.array(stats.lognorm.fit(x, f0=0.25, loc=0))
|
|
assert_almost_equal(params, true, decimal=2)
|
|
|
|
params = np.array(stats.lognorm.fit(x, f0=0.25, floc=0))
|
|
assert_almost_equal(params, true, decimal=2)
|
|
|
|
np.random.seed(5678)
|
|
loc = 1
|
|
floc = 0.9
|
|
x = stats.norm.rvs(loc, 2., size=100)
|
|
params = np.array(stats.norm.fit(x, floc=floc))
|
|
expected = np.array([floc, np.sqrt(((x-floc)**2).mean())])
|
|
assert_almost_equal(params, expected, decimal=4)
|
|
|
|
|
|
def test_regression_ticket_1530():
|
|
# Check the starting value works for Cauchy distribution fit.
|
|
np.random.seed(654321)
|
|
rvs = stats.cauchy.rvs(size=100)
|
|
params = stats.cauchy.fit(rvs)
|
|
expected = (0.045, 1.142)
|
|
assert_almost_equal(params, expected, decimal=1)
|
|
|
|
|
|
def test_tukeylambda_stats_ticket_1545():
|
|
# Some test for the variance and kurtosis of the Tukey Lambda distr.
|
|
# See test_tukeylamdba_stats.py for more tests.
|
|
|
|
mv = stats.tukeylambda.stats(0, moments='mvsk')
|
|
# Known exact values:
|
|
expected = [0, np.pi**2/3, 0, 1.2]
|
|
assert_almost_equal(mv, expected, decimal=10)
|
|
|
|
mv = stats.tukeylambda.stats(3.13, moments='mvsk')
|
|
# 'expected' computed with mpmath.
|
|
expected = [0, 0.0269220858861465102, 0, -0.898062386219224104]
|
|
assert_almost_equal(mv, expected, decimal=10)
|
|
|
|
mv = stats.tukeylambda.stats(0.14, moments='mvsk')
|
|
# 'expected' computed with mpmath.
|
|
expected = [0, 2.11029702221450250, 0, -0.02708377353223019456]
|
|
assert_almost_equal(mv, expected, decimal=10)
|
|
|
|
|
|
def test_poisson_logpmf_ticket_1436():
|
|
assert_(np.isfinite(stats.poisson.logpmf(1500, 200)))
|
|
|
|
|
|
def test_powerlaw_stats():
|
|
"""Test the powerlaw stats function.
|
|
|
|
This unit test is also a regression test for ticket 1548.
|
|
|
|
The exact values are:
|
|
mean:
|
|
mu = a / (a + 1)
|
|
variance:
|
|
sigma**2 = a / ((a + 2) * (a + 1) ** 2)
|
|
skewness:
|
|
One formula (see http://en.wikipedia.org/wiki/Skewness) is
|
|
gamma_1 = (E[X**3] - 3*mu*E[X**2] + 2*mu**3) / sigma**3
|
|
A short calculation shows that E[X**k] is a / (a + k), so gamma_1
|
|
can be implemented as
|
|
n = a/(a+3) - 3*(a/(a+1))*a/(a+2) + 2*(a/(a+1))**3
|
|
d = sqrt(a/((a+2)*(a+1)**2)) ** 3
|
|
gamma_1 = n/d
|
|
Either by simplifying, or by a direct calculation of mu_3 / sigma**3,
|
|
one gets the more concise formula:
|
|
gamma_1 = -2.0 * ((a - 1) / (a + 3)) * sqrt((a + 2) / a)
|
|
kurtosis: (See http://en.wikipedia.org/wiki/Kurtosis)
|
|
The excess kurtosis is
|
|
gamma_2 = mu_4 / sigma**4 - 3
|
|
A bit of calculus and algebra (sympy helps) shows that
|
|
mu_4 = 3*a*(3*a**2 - a + 2) / ((a+1)**4 * (a+2) * (a+3) * (a+4))
|
|
so
|
|
gamma_2 = 3*(3*a**2 - a + 2) * (a+2) / (a*(a+3)*(a+4)) - 3
|
|
which can be rearranged to
|
|
gamma_2 = 6 * (a**3 - a**2 - 6*a + 2) / (a*(a+3)*(a+4))
|
|
"""
|
|
cases = [(1.0, (0.5, 1./12, 0.0, -1.2)),
|
|
(2.0, (2./3, 2./36, -0.56568542494924734, -0.6))]
|
|
for a, exact_mvsk in cases:
|
|
mvsk = stats.powerlaw.stats(a, moments="mvsk")
|
|
assert_array_almost_equal(mvsk, exact_mvsk)
|
|
|
|
|
|
def test_ksone_fit_freeze():
|
|
# Regression test for ticket #1638.
|
|
d = np.array(
|
|
[-0.18879233, 0.15734249, 0.18695107, 0.27908787, -0.248649,
|
|
-0.2171497, 0.12233512, 0.15126419, 0.03119282, 0.4365294,
|
|
0.08930393, -0.23509903, 0.28231224, -0.09974875, -0.25196048,
|
|
0.11102028, 0.1427649, 0.10176452, 0.18754054, 0.25826724,
|
|
0.05988819, 0.0531668, 0.21906056, 0.32106729, 0.2117662,
|
|
0.10886442, 0.09375789, 0.24583286, -0.22968366, -0.07842391,
|
|
-0.31195432, -0.21271196, 0.1114243, -0.13293002, 0.01331725,
|
|
-0.04330977, -0.09485776, -0.28434547, 0.22245721, -0.18518199,
|
|
-0.10943985, -0.35243174, 0.06897665, -0.03553363, -0.0701746,
|
|
-0.06037974, 0.37670779, -0.21684405])
|
|
|
|
try:
|
|
olderr = np.seterr(invalid='ignore')
|
|
with warnings.catch_warnings():
|
|
warnings.simplefilter('ignore', UserWarning)
|
|
warnings.simplefilter('ignore', RuntimeWarning)
|
|
stats.ksone.fit(d)
|
|
finally:
|
|
np.seterr(**olderr)
|
|
|
|
|
|
def test_norm_logcdf():
|
|
# Test precision of the logcdf of the normal distribution.
|
|
# This precision was enhanced in ticket 1614.
|
|
x = -np.asarray(list(range(0, 120, 4)))
|
|
# Values from R
|
|
expected = [-0.69314718, -10.36010149, -35.01343716, -75.41067300,
|
|
-131.69539607, -203.91715537, -292.09872100, -396.25241451,
|
|
-516.38564863, -652.50322759, -804.60844201, -972.70364403,
|
|
-1156.79057310, -1356.87055173, -1572.94460885, -1805.01356068,
|
|
-2053.07806561, -2317.13866238, -2597.19579746, -2893.24984493,
|
|
-3205.30112136, -3533.34989701, -3877.39640444, -4237.44084522,
|
|
-4613.48339520, -5005.52420869, -5413.56342187, -5837.60115548,
|
|
-6277.63751711, -6733.67260303]
|
|
|
|
olderr = np.seterr(divide='ignore')
|
|
try:
|
|
assert_allclose(stats.norm().logcdf(x), expected, atol=1e-8)
|
|
finally:
|
|
np.seterr(**olderr)
|
|
|
|
|
|
def test_levy_cdf_ppf():
|
|
# Test levy.cdf, including small arguments.
|
|
x = np.array([1000, 1.0, 0.5, 0.1, 0.01, 0.001])
|
|
|
|
# Expected values were calculated separately with mpmath.
|
|
# E.g.
|
|
# >>> mpmath.mp.dps = 100
|
|
# >>> x = mpmath.mp.mpf('0.01')
|
|
# >>> cdf = mpmath.erfc(mpmath.sqrt(1/(2*x)))
|
|
expected = np.array([0.9747728793699604,
|
|
0.3173105078629141,
|
|
0.1572992070502851,
|
|
0.0015654022580025495,
|
|
1.523970604832105e-23,
|
|
1.795832784800726e-219])
|
|
|
|
y = stats.levy.cdf(x)
|
|
assert_allclose(y, expected, rtol=1e-10)
|
|
|
|
# ppf(expected) should get us back to x.
|
|
xx = stats.levy.ppf(expected)
|
|
assert_allclose(xx, x, rtol=1e-13)
|
|
|
|
|
|
def test_hypergeom_interval_1802():
|
|
# these two had endless loops
|
|
assert_equal(stats.hypergeom.interval(.95, 187601, 43192, 757),
|
|
(152.0, 197.0))
|
|
assert_equal(stats.hypergeom.interval(.945, 187601, 43192, 757),
|
|
(152.0, 197.0))
|
|
# this was working also before
|
|
assert_equal(stats.hypergeom.interval(.94, 187601, 43192, 757),
|
|
(153.0, 196.0))
|
|
|
|
# degenerate case .a == .b
|
|
assert_equal(stats.hypergeom.ppf(0.02, 100, 100, 8), 8)
|
|
assert_equal(stats.hypergeom.ppf(1, 100, 100, 8), 8)
|
|
|
|
|
|
def test_distribution_too_many_args():
|
|
# Check that a TypeError is raised when too many args are given to a method
|
|
# Regression test for ticket 1815.
|
|
x = np.linspace(0.1, 0.7, num=5)
|
|
assert_raises(TypeError, stats.gamma.pdf, x, 2, 3, loc=1.0)
|
|
assert_raises(TypeError, stats.gamma.pdf, x, 2, 3, 4, loc=1.0)
|
|
assert_raises(TypeError, stats.gamma.pdf, x, 2, 3, 4, 5)
|
|
assert_raises(TypeError, stats.gamma.pdf, x, 2, 3, loc=1.0, scale=0.5)
|
|
assert_raises(TypeError, stats.gamma.rvs, 2., 3, loc=1.0, scale=0.5)
|
|
assert_raises(TypeError, stats.gamma.cdf, x, 2., 3, loc=1.0, scale=0.5)
|
|
assert_raises(TypeError, stats.gamma.ppf, x, 2., 3, loc=1.0, scale=0.5)
|
|
assert_raises(TypeError, stats.gamma.stats, 2., 3, loc=1.0, scale=0.5)
|
|
assert_raises(TypeError, stats.gamma.entropy, 2., 3, loc=1.0, scale=0.5)
|
|
assert_raises(TypeError, stats.gamma.fit, x, 2., 3, loc=1.0, scale=0.5)
|
|
|
|
# These should not give errors
|
|
stats.gamma.pdf(x, 2, 3) # loc=3
|
|
stats.gamma.pdf(x, 2, 3, 4) # loc=3, scale=4
|
|
stats.gamma.stats(2., 3)
|
|
stats.gamma.stats(2., 3, 4)
|
|
stats.gamma.stats(2., 3, 4, 'mv')
|
|
stats.gamma.rvs(2., 3, 4, 5)
|
|
stats.gamma.fit(stats.gamma.rvs(2., size=7), 2.)
|
|
|
|
# Also for a discrete distribution
|
|
stats.geom.pmf(x, 2, loc=3) # no error, loc=3
|
|
assert_raises(TypeError, stats.geom.pmf, x, 2, 3, 4)
|
|
assert_raises(TypeError, stats.geom.pmf, x, 2, 3, loc=4)
|
|
|
|
# And for distributions with 0, 2 and 3 args respectively
|
|
assert_raises(TypeError, stats.expon.pdf, x, 3, loc=1.0)
|
|
assert_raises(TypeError, stats.exponweib.pdf, x, 3, 4, 5, loc=1.0)
|
|
assert_raises(TypeError, stats.exponweib.pdf, x, 3, 4, 5, 0.1, 0.1)
|
|
assert_raises(TypeError, stats.ncf.pdf, x, 3, 4, 5, 6, loc=1.0)
|
|
assert_raises(TypeError, stats.ncf.pdf, x, 3, 4, 5, 6, 1.0, scale=0.5)
|
|
stats.ncf.pdf(x, 3, 4, 5, 6, 1.0) # 3 args, plus loc/scale
|
|
|
|
|
|
def test_ncx2_tails_ticket_955():
|
|
# Trac #955 -- check that the cdf computed by special functions
|
|
# matches the integrated pdf
|
|
a = stats.ncx2.cdf(np.arange(20, 25, 0.2), 2, 1.07458615e+02)
|
|
b = stats.ncx2._cdfvec(np.arange(20, 25, 0.2), 2, 1.07458615e+02)
|
|
assert_allclose(a, b, rtol=1e-3, atol=0)
|
|
|
|
|
|
def test_foldnorm_zero():
|
|
# Parameter value c=0 was not enabled, see gh-2399.
|
|
rv = stats.foldnorm(0, scale=1)
|
|
assert_equal(rv.cdf(0), 0) # rv.cdf(0) previously resulted in: nan
|
|
|
|
|
|
def test_stats_shapes_argcheck():
|
|
# stats method was failing for vector shapes if some of the values
|
|
# were outside of the allowed range, see gh-2678
|
|
mv3 = stats.invgamma.stats([0.0, 0.5, 1.0], 1, 0.5) # 0 is not a legal `a`
|
|
mv2 = stats.invgamma.stats([0.5, 1.0], 1, 0.5)
|
|
mv2_augmented = tuple(np.r_[np.nan, _] for _ in mv2)
|
|
assert_equal(mv2_augmented, mv3)
|
|
|
|
mv3 = stats.lognorm.stats([2, 2.4, -1]) # -1 is not a legal shape parameter
|
|
mv2 = stats.lognorm.stats([2, 2.4])
|
|
mv2_augmented = tuple(np.r_[_, np.nan] for _ in mv2)
|
|
assert_equal(mv2_augmented, mv3)
|
|
|
|
# FIXME: this is only a quick-and-dirty test of a quick-and-dirty bugfix.
|
|
# stats method with multiple shape parameters is not properly vectorized
|
|
# anyway, so some distributions may or may not fail.
|
|
|
|
|
|
## Test subclassing distributions w/ explicit shapes
|
|
|
|
class _distr_gen(stats.rv_continuous):
|
|
def _pdf(self, x, a):
|
|
return 42
|
|
|
|
|
|
class _distr2_gen(stats.rv_continuous):
|
|
def _cdf(self, x, a):
|
|
return 42 * a + x
|
|
|
|
|
|
class _distr3_gen(stats.rv_continuous):
|
|
def _pdf(self, x, a, b):
|
|
return a + b
|
|
|
|
def _cdf(self, x, a):
|
|
# Different # of shape params from _pdf, to be able to check that
|
|
# inspection catches the inconsistency."""
|
|
return 42 * a + x
|
|
|
|
|
|
class _distr6_gen(stats.rv_continuous):
|
|
# Two shape parameters (both _pdf and _cdf defined, consistent shapes.)
|
|
def _pdf(self, x, a, b):
|
|
return a*x + b
|
|
|
|
def _cdf(self, x, a, b):
|
|
return 42 * a + x
|
|
|
|
|
|
class TestSubclassingExplicitShapes(TestCase):
|
|
# Construct a distribution w/ explicit shapes parameter and test it.
|
|
|
|
def test_correct_shapes(self):
|
|
dummy_distr = _distr_gen(name='dummy', shapes='a')
|
|
assert_equal(dummy_distr.pdf(1, a=1), 42)
|
|
|
|
def test_wrong_shapes_1(self):
|
|
dummy_distr = _distr_gen(name='dummy', shapes='A')
|
|
assert_raises(TypeError, dummy_distr.pdf, 1, **dict(a=1))
|
|
|
|
def test_wrong_shapes_2(self):
|
|
dummy_distr = _distr_gen(name='dummy', shapes='a, b, c')
|
|
dct = dict(a=1, b=2, c=3)
|
|
assert_raises(TypeError, dummy_distr.pdf, 1, **dct)
|
|
|
|
def test_shapes_string(self):
|
|
# shapes must be a string
|
|
dct = dict(name='dummy', shapes=42)
|
|
assert_raises(TypeError, _distr_gen, **dct)
|
|
|
|
def test_shapes_identifiers_1(self):
|
|
# shapes must be a comma-separated list of valid python identifiers
|
|
dct = dict(name='dummy', shapes='(!)')
|
|
assert_raises(SyntaxError, _distr_gen, **dct)
|
|
|
|
def test_shapes_identifiers_2(self):
|
|
dct = dict(name='dummy', shapes='4chan')
|
|
assert_raises(SyntaxError, _distr_gen, **dct)
|
|
|
|
def test_shapes_identifiers_3(self):
|
|
dct = dict(name='dummy', shapes='m(fti)')
|
|
assert_raises(SyntaxError, _distr_gen, **dct)
|
|
|
|
def test_shapes_identifiers_nodefaults(self):
|
|
dct = dict(name='dummy', shapes='a=2')
|
|
assert_raises(SyntaxError, _distr_gen, **dct)
|
|
|
|
def test_shapes_args(self):
|
|
dct = dict(name='dummy', shapes='*args')
|
|
assert_raises(SyntaxError, _distr_gen, **dct)
|
|
|
|
def test_shapes_kwargs(self):
|
|
dct = dict(name='dummy', shapes='**kwargs')
|
|
assert_raises(SyntaxError, _distr_gen, **dct)
|
|
|
|
def test_shapes_keywords(self):
|
|
# python keywords cannot be used for shape parameters
|
|
dct = dict(name='dummy', shapes='a, b, c, lambda')
|
|
assert_raises(SyntaxError, _distr_gen, **dct)
|
|
|
|
def test_shapes_signature(self):
|
|
# test explicit shapes which agree w/ the signature of _pdf
|
|
class _dist_gen(stats.rv_continuous):
|
|
def _pdf(self, x, a):
|
|
return stats.norm._pdf(x) * a
|
|
|
|
dist = _dist_gen(shapes='a')
|
|
assert_equal(dist.pdf(0.5, a=2), stats.norm.pdf(0.5)*2)
|
|
|
|
def test_shapes_signature_inconsistent(self):
|
|
# test explicit shapes which do not agree w/ the signature of _pdf
|
|
class _dist_gen(stats.rv_continuous):
|
|
def _pdf(self, x, a):
|
|
return stats.norm._pdf(x) * a
|
|
|
|
dist = _dist_gen(shapes='a, b')
|
|
assert_raises(TypeError, dist.pdf, 0.5, **dict(a=1, b=2))
|
|
|
|
def test_star_args(self):
|
|
# test _pdf with only starargs
|
|
# NB: **kwargs of pdf will never reach _pdf
|
|
class _dist_gen(stats.rv_continuous):
|
|
def _pdf(self, x, *args):
|
|
extra_kwarg = args[0]
|
|
return stats.norm._pdf(x) * extra_kwarg
|
|
|
|
dist = _dist_gen(shapes='extra_kwarg')
|
|
assert_equal(dist.pdf(0.5, extra_kwarg=33), stats.norm.pdf(0.5)*33)
|
|
assert_equal(dist.pdf(0.5, 33), stats.norm.pdf(0.5)*33)
|
|
assert_raises(TypeError, dist.pdf, 0.5, **dict(xxx=33))
|
|
|
|
def test_star_args_2(self):
|
|
# test _pdf with named & starargs
|
|
# NB: **kwargs of pdf will never reach _pdf
|
|
class _dist_gen(stats.rv_continuous):
|
|
def _pdf(self, x, offset, *args):
|
|
extra_kwarg = args[0]
|
|
return stats.norm._pdf(x) * extra_kwarg + offset
|
|
|
|
dist = _dist_gen(shapes='offset, extra_kwarg')
|
|
assert_equal(dist.pdf(0.5, offset=111, extra_kwarg=33),
|
|
stats.norm.pdf(0.5)*33 + 111)
|
|
assert_equal(dist.pdf(0.5, 111, 33),
|
|
stats.norm.pdf(0.5)*33 + 111)
|
|
|
|
def test_extra_kwarg(self):
|
|
# **kwargs to _pdf are ignored.
|
|
# this is a limitation of the framework (_pdf(x, *goodargs))
|
|
class _distr_gen(stats.rv_continuous):
|
|
def _pdf(self, x, *args, **kwargs):
|
|
# _pdf should handle *args, **kwargs itself. Here "handling" is
|
|
# ignoring *args and looking for ``extra_kwarg`` and using that.
|
|
extra_kwarg = kwargs.pop('extra_kwarg', 1)
|
|
return stats.norm._pdf(x) * extra_kwarg
|
|
|
|
dist = _distr_gen(shapes='extra_kwarg')
|
|
assert_equal(dist.pdf(1, extra_kwarg=3), stats.norm.pdf(1))
|
|
|
|
def shapes_empty_string(self):
|
|
# shapes='' is equivalent to shapes=None
|
|
class _dist_gen(stats.rv_continuous):
|
|
def _pdf(self, x):
|
|
return stats.norm.pdf(x)
|
|
|
|
dist = _dist_gen(shapes='')
|
|
assert_equal(dist.pdf(0.5), stats.norm.pdf(0.5))
|
|
|
|
|
|
class TestSubclassingNoShapes(TestCase):
|
|
# Construct a distribution w/o explicit shapes parameter and test it.
|
|
|
|
def test_only__pdf(self):
|
|
dummy_distr = _distr_gen(name='dummy')
|
|
assert_equal(dummy_distr.pdf(1, a=1), 42)
|
|
|
|
def test_only__cdf(self):
|
|
# _pdf is determined from _cdf by taking numerical derivative
|
|
dummy_distr = _distr2_gen(name='dummy')
|
|
assert_almost_equal(dummy_distr.pdf(1, a=1), 1)
|
|
|
|
@dec.skipif(DOCSTRINGS_STRIPPED)
|
|
def test_signature_inspection(self):
|
|
# check that _pdf signature inspection works correctly, and is used in
|
|
# the class docstring
|
|
dummy_distr = _distr_gen(name='dummy')
|
|
assert_equal(dummy_distr.numargs, 1)
|
|
assert_equal(dummy_distr.shapes, 'a')
|
|
res = re.findall('logpdf\(x, a, loc=0, scale=1\)',
|
|
dummy_distr.__doc__)
|
|
assert_(len(res) == 1)
|
|
|
|
@dec.skipif(DOCSTRINGS_STRIPPED)
|
|
def test_signature_inspection_2args(self):
|
|
# same for 2 shape params and both _pdf and _cdf defined
|
|
dummy_distr = _distr6_gen(name='dummy')
|
|
assert_equal(dummy_distr.numargs, 2)
|
|
assert_equal(dummy_distr.shapes, 'a, b')
|
|
res = re.findall('logpdf\(x, a, b, loc=0, scale=1\)',
|
|
dummy_distr.__doc__)
|
|
assert_(len(res) == 1)
|
|
|
|
def test_signature_inspection_2args_incorrect_shapes(self):
|
|
# both _pdf and _cdf defined, but shapes are inconsistent: raises
|
|
try:
|
|
_distr3_gen(name='dummy')
|
|
except TypeError:
|
|
pass
|
|
else:
|
|
raise AssertionError('TypeError not raised.')
|
|
|
|
def test_defaults_raise(self):
|
|
# default arguments should raise
|
|
class _dist_gen(stats.rv_continuous):
|
|
def _pdf(self, x, a=42):
|
|
return 42
|
|
assert_raises(TypeError, _dist_gen, **dict(name='dummy'))
|
|
|
|
def test_starargs_raise(self):
|
|
# without explicit shapes, *args are not allowed
|
|
class _dist_gen(stats.rv_continuous):
|
|
def _pdf(self, x, a, *args):
|
|
return 42
|
|
assert_raises(TypeError, _dist_gen, **dict(name='dummy'))
|
|
|
|
def test_kwargs_raise(self):
|
|
# without explicit shapes, **kwargs are not allowed
|
|
class _dist_gen(stats.rv_continuous):
|
|
def _pdf(self, x, a, **kwargs):
|
|
return 42
|
|
assert_raises(TypeError, _dist_gen, **dict(name='dummy'))
|
|
|
|
|
|
@dec.skipif(DOCSTRINGS_STRIPPED)
|
|
def test_docstrings():
|
|
badones = [',\s*,', '\(\s*,', '^\s*:']
|
|
for distname in stats.__all__:
|
|
dist = getattr(stats, distname)
|
|
if isinstance(dist, (stats.rv_discrete, stats.rv_continuous)):
|
|
for regex in badones:
|
|
assert_(re.search(regex, dist.__doc__) is None)
|
|
|
|
|
|
def test_infinite_input():
|
|
assert_almost_equal(stats.skellam.sf(np.inf, 10, 11), 0)
|
|
assert_almost_equal(stats.ncx2._cdf(np.inf, 8, 0.1), 1)
|
|
|
|
|
|
if __name__ == "__main__":
|
|
run_module_suite()
|