From 111d6ce8084e93fb8adffc6101a828228aeb3cb5 Mon Sep 17 00:00:00 2001 From: "Per.Andreas.Brodtkorb" Date: Tue, 23 Nov 2010 08:45:25 +0000 Subject: [PATCH] Added test_distributions.py updated test_estimation.py --- .../wafo/stats/tests/test_distributions.py | 567 ++++++++++++++++++ .../src/wafo/stats/tests/test_estimation.py | 59 +- 2 files changed, 588 insertions(+), 38 deletions(-) create mode 100644 pywafo/src/wafo/stats/tests/test_distributions.py diff --git a/pywafo/src/wafo/stats/tests/test_distributions.py b/pywafo/src/wafo/stats/tests/test_distributions.py new file mode 100644 index 0000000..1df98ef --- /dev/null +++ b/pywafo/src/wafo/stats/tests/test_distributions.py @@ -0,0 +1,567 @@ +""" Test functions for stats module + +""" + +from numpy.testing import TestCase, run_module_suite, assert_equal, \ + assert_array_equal, assert_almost_equal, assert_array_almost_equal, \ + assert_, rand, dec + + +import numpy +import numpy as np +from numpy import typecodes, array +import wafo.stats as stats +from wafo.stats.distributions import argsreduce + +def kolmogorov_check(diststr, args=(), N=20, significance=0.01): + qtest = stats.ksoneisf(significance, N) + cdf = eval('stats.'+diststr+'.cdf') + dist = eval('stats.'+diststr) + # Get random numbers + kwds = {'size':N} + vals = numpy.sort(dist.rvs(*args, **kwds)) + cdfvals = cdf(vals, *args) + q = max(abs(cdfvals - np.arange(1.0, N+1)/N)) + assert_(q < qtest, msg="Failed q=%f, bound=%f, alpha=%f" % (q, qtest, significance)) + return + + +# generate test cases to test cdf and distribution consistency +dists = ['uniform','norm','lognorm','expon','beta', + 'powerlaw','bradford','burr','fisk','cauchy','halfcauchy', + 'foldcauchy','gamma','gengamma','loggamma', + 'alpha','anglit','arcsine','betaprime','erlang', + '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'] + +# 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 == 'erlang': + args = (4,)+tuple(rand(2)) + elif 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 + +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=`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']) + + +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']) + +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']) + +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_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) + + +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']) + +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 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']) + +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']) + +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']) + +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)) + +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 TestHypergeom(TestCase): + 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) + + +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()) + +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 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(TestCase): + skip = ['ncf'] + + @dec.slow + def test_fit(self): + for func, dist, args, alpha in test_all_distributions(): + if dist in self.skip: + continue + distfunc = getattr(stats, dist) + 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 in ['erlang', 'frechet']: + assert_(len(vals)==len(args)) + assert_(len(vals2)==len(args)) + else: + assert_(len(vals) == 2+len(args)) + assert_(len(vals2)==2+len(args)) + + @dec.slow + def test_fix_fit(self): + for func, dist, args, alpha in test_all_distributions(): + # Not sure why 'ncf', and 'beta' are failing + # erlang and frechet have different len(args) than distfunc.numargs + if dist in self.skip + ['erlang', 'frechet', 'beta']: + continue + distfunc = getattr(stats, dist) + 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]) + +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) + 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_02(self): + """Regression test for ticket #1293.""" + # 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. + s = 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_regression_ticket_1316(): + """Regression test for 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.distributions.gamma_gen(name='gamma') + + +if __name__ == "__main__": + run_module_suite() diff --git a/pywafo/src/wafo/stats/tests/test_estimation.py b/pywafo/src/wafo/stats/tests/test_estimation.py index 107c5b5..e609532 100644 --- a/pywafo/src/wafo/stats/tests/test_estimation.py +++ b/pywafo/src/wafo/stats/tests/test_estimation.py @@ -1,65 +1,48 @@ -''' +# -*- coding:utf-8 -*- +""" Created on 19. nov. 2010 @author: pab -''' - -import wafo.stats as ws -from wafo.stats.estimation import Profile, FitDistribution -from numpy import log, array -def test_fit_and_profile(): +""" +from numpy import array, log +import wafo.stats as ws +from wafo.stats.estimation import FitDistribution, Profile +def test_profile(): ''' - # MLE + # MLE import wafo.stats as ws - R = ws.weibull_min.rvs(1,size=100) - >>> R = array([ 1.11214852, 0.55730247, 0.88241032, 0.13421021, 0.21939628, - ... 0.93901239, 4.55065051, 0.06728561, 1.65808964, 3.31725493, - ... 1.26022136, 0.48707756, 1.04600817, 3.45064509, 0.23478845, - ... 0.73197928, 0.78702608, 1.10566899, 0.08770387, 0.29773143, - ... 0.442589 , 0.06017045, 0.99086215, 0.36520787, 0.55931192, - ... 0.02220038, 0.09683231, 0.2407798 , 1.5550326 , 1.21834857, - ... 0.51317074, 3.12672982, 2.01984527, 1.42837289, 1.60665276, - ... 0.07643233, 0.97254151, 0.71870616, 0.15669964, 0.7647863 , - ... 2.15586998, 1.03167507, 0.66384501, 0.26824186, 0.25953582, - ... 0.52773641, 0.3371869 , 1.17897776, 0.36094868, 0.49593992, - ... 0.12124638, 0.13178139, 0.28325022, 1.2189952 , 2.85895075, - ... 1.94868469, 1.5216953 , 0.77066953, 0.42610845, 1.52848366, - ... 0.26570905, 0.52072509, 1.03106309, 0.34430637, 1.36386879, - ... 1.05857501, 0.66754409, 3.18836655, 0.40063524, 0.40696207, - ... 1.24167978, 0.27927484, 3.95644924, 1.34793436, 0.66225993, - ... 0.46567866, 0.09325875, 1.90895739, 1.23831713, 0.29819745, - ... 1.10528719, 0.01973962, 0.56370415, 0.66831751, 1.33781636, - ... 0.53985288, 0.85251279, 1.75369891, 0.08066507, 0.29273944, - ... 0.38886539, 0.20961546, 0.93814728, 0.08732068, 1.6336713 , - ... 1.38893923, 2.05882459, 0.51709862, 0.37553027, 0.06829269]) + R = ws.weibull_min.rvs(1,size=20) + >>> R = array([ 0.08018795, 1.09015299, 2.08591271, 0.51542081, 0.75692042, + ... 0.57837017, 0.45419753, 1.1559131 , 0.26582267, 0.51517273, + ... 0.75008492, 0.59404957, 1.33748264, 0.14472142, 0.77459603, + ... 1.77312556, 1.06347991, 0.42007769, 0.71094628, 0.02366977]) >>> phat = FitDistribution(ws.weibull_min, R, 1, scale=1, floc=0.0) >>> phat.par - array([ 1.07424857, 0. , 0.96282996]) + array([ 1.37836487, 0. , 0.82085633]) # Better CI for phat.par[i=0] >>> Lp = Profile(phat, i=0) >>> Lp.get_bounds(alpha=0.1) - array([ 0.94231243, 1.21444458]) + array([ 1.00225064, 1.8159036 ]) >>> SF = 1./990 >>> x = phat.isf(SF) >>> x - 5.8114656626008818 + 3.3323076459875312 # CI for x >>> Lx = phat.profile(i=0, x=x, link=phat.dist.link) >>> Lx.get_bounds(alpha=0.2) - array([ 4.91172017, 7.08081584]) + array([ 2.52211661, 4.98664787]) # CI for logSF=log(SF) - >>> Lsf = phat.profile(i=0, logSF=log(SF), link=phat.dist.link) + >>> logSF = log(SF) + >>> Lsf = phat.profile(i=0, logSF=logSF, link=phat.dist.link, pmin=logSF-10,pmax=logSF+5) >>> Lsf.get_bounds(alpha=0.2) - array([-8.37580767, -5.66897775]) + array([-10.87488318, -4.36225468]) ''' - pass - -if __name__ == "__main__": +if __name__ == '__main__': import doctest doctest.testmod() \ No newline at end of file