From d540b8d07829af6317c3992b3b1a9c11fffcabeb Mon Sep 17 00:00:00 2001 From: "Per.Andreas.Brodtkorb" Date: Thu, 4 Oct 2012 02:20:21 +0000 Subject: [PATCH] Updated distributions.py --- pywafo/src/wafo/integrate.py | 60 +-- pywafo/src/wafo/stats/distributions.py | 598 ++++++++++++++++--------- 2 files changed, 410 insertions(+), 248 deletions(-) diff --git a/pywafo/src/wafo/integrate.py b/pywafo/src/wafo/integrate.py index 58aead7..0bd336a 100644 --- a/pywafo/src/wafo/integrate.py +++ b/pywafo/src/wafo/integrate.py @@ -1136,7 +1136,7 @@ def richardson(Q, k): R = Q[k] + (Q[k] - Q[k - 1]) / c return R -def quadgr(fun, a, b, abseps=1e-5, maxiter=17): +def quadgr(fun, a, b, abseps=1e-5, max_iter=17): ''' Gauss-Legendre quadrature with Richardson extrapolation. @@ -1446,35 +1446,35 @@ def qdemo(f, a, b): def main(): - val, err = clencurt(np.exp, 0, 2) - valt = np.exp(2) - np.exp(0) - [Q, err] = quadgr(lambda x: x ** 2, 1, 4, 1e-9) - [Q, err] = quadgr(humps, 1, 4, 1e-9) - - [x, w] = h_roots(11, 'newton') - sum(w) - [x2, w2] = la_roots(11, 1, 't') - - from scitools import numpyutils as npu #@UnresolvedImport - fun = npu.wrap2callable('x**2') - p0 = fun(0) - A = [0, 1, 1]; B = [2, 4, 3] - area, err = gaussq(fun, A, B) - - fun = npu.wrap2callable('x**2') - [val1, err1] = gaussq(fun, A, B) - - - #Integration of x^2*exp(-x) from zero to infinity: - fun2 = npu.wrap2callable('1') - [val2, err2] = gaussq(fun2, 0, np.inf, wfun=3, alpha=2) - [val2, err2] = gaussq(lambda x: x ** 2, 0, np.inf, wfun=3, alpha=0) - - #Integrate humps from 0 to 2 and from 1 to 4 - [val3, err3] = gaussq(humps, A, B) - - [x, w] = p_roots(11, 'newton', 1, 3) - y = np.sum(x ** 2 * w) +# val, err = clencurt(np.exp, 0, 2) +# valt = np.exp(2) - np.exp(0) +# [Q, err] = quadgr(lambda x: x ** 2, 1, 4, 1e-9) +# [Q, err] = quadgr(humps, 1, 4, 1e-9) +# +# [x, w] = h_roots(11, 'newton') +# sum(w) +# [x2, w2] = la_roots(11, 1, 't') +# +# from scitools import numpyutils as npu #@UnresolvedImport +# fun = npu.wrap2callable('x**2') +# p0 = fun(0) +# A = [0, 1, 1]; B = [2, 4, 3] +# area, err = gaussq(fun, A, B) +# +# fun = npu.wrap2callable('x**2') +# [val1, err1] = gaussq(fun, A, B) +# +# +# #Integration of x^2*exp(-x) from zero to infinity: +# fun2 = npu.wrap2callable('1') +# [val2, err2] = gaussq(fun2, 0, np.inf, wfun=3, alpha=2) +# [val2, err2] = gaussq(lambda x: x ** 2, 0, np.inf, wfun=3, alpha=0) +# +# #Integrate humps from 0 to 2 and from 1 to 4 +# [val3, err3] = gaussq(humps, A, B) +# +# [x, w] = p_roots(11, 'newton', 1, 3) +# y = np.sum(x ** 2 * w) x = np.linspace(0, np.pi / 2) q0 = np.trapz(humps(x), x) diff --git a/pywafo/src/wafo/stats/distributions.py b/pywafo/src/wafo/stats/distributions.py index d1d4e89..da9f189 100644 --- a/pywafo/src/wafo/stats/distributions.py +++ b/pywafo/src/wafo/stats/distributions.py @@ -58,10 +58,9 @@ __all__ = [ 'rdist', 'rayleigh', 'truncrayleigh','reciprocal', 'rice', 'recipinvgauss', 'semicircular', 'triang', 'truncexpon', 'truncnorm', 'tukeylambda', 'uniform', 'vonmises', 'wald', 'wrapcauchy', - 'entropy', 'rv_discrete', - 'binom', 'bernoulli', 'nbinom', 'geom', 'hypergeom', 'logser', - 'poisson', 'planck', 'boltzmann', 'randint', 'zipf', 'dlaplace', - 'skellam' + 'entropy', 'rv_discrete', 'binom', 'bernoulli', 'nbinom', 'geom', + 'hypergeom', 'logser', 'poisson', 'planck', 'boltzmann', 'randint', + 'zipf', 'dlaplace', 'skellam' ] floatinfo = numpy.finfo(float) @@ -241,6 +240,8 @@ Display frozen pdf >>> x = np.linspace(0, np.minimum(rv.dist.b, 3)) >>> h = plt.plot(x, rv.pdf(x)) +Here, ``rv.dist.b`` is the right endpoint of the support of ``rv.dist``. + Check accuracy of cdf and ppf >>> prb = %(name)s.cdf(x, %(shapes)s) @@ -1165,8 +1166,8 @@ class rv_generic(object): arg1, arg2, ... : array_like The shape parameter(s) for the distribution (see docstring of the instance object for more information) - loc: array_like, optioal - location parameter (deafult = 0) + loc: array_like, optional + location parameter (default = 0) scale : array_like, optional scale paramter (default = 1) @@ -1204,9 +1205,9 @@ class rv_continuous(rv_generic): Upper bound of the support of the distribution, default is plus infinity. xa : float, optional - Lower bound for fixed point calculation for generic ppf. + DEPRECATED xb : float, optional - Upper bound for fixed point calculation for generic ppf. + DEPRECATED xtol : float, optional The tolerance for fixed point calculation for generic ppf. badvalue : object, optional @@ -1438,7 +1439,7 @@ class rv_continuous(rv_generic): """ - def __init__(self, momtype=1, a=None, b=None, xa=-10.0, xb=10.0, + def __init__(self, momtype=1, a=None, b=None, xa=None, xb=None, xtol=1e-14, badvalue=None, name=None, longname=None, shapes=None, extradoc=None): @@ -1457,6 +1458,12 @@ class rv_continuous(rv_generic): self.a = -inf if b is None: self.b = inf + if xa is not None: + warnings.warn("The `xa` parameter is deprecated and will be " + "removed in scipy 0.12", DeprecationWarning) + if xb is not None: + warnings.warn("The `xb` parameter is deprecated and will be " + "removed in scipy 0.12", DeprecationWarning) self.xa = xa self.xb = xb self.xtol = xtol @@ -1539,7 +1546,28 @@ class rv_continuous(rv_generic): return apply(self.cdf, (x, )+args)-q def _ppf_single_call(self, q, *args): - return optimize.brentq(self._ppf_to_solve, self.xa, self.xb, args=(q,)+args, xtol=self.xtol) + left = right = None + if self.a > -np.inf: + left = self.a + if self.b < np.inf: + right = self.b + + factor = 10. + if not left: # i.e. self.a = -inf + left = -1.*factor + while self._ppf_to_solve(left, q,*args) > 0.: + right = left + left *= factor + # left is now such that cdf(left) < q + if not right: # i.e. self.b = inf + right = factor + while self._ppf_to_solve(right, q,*args) < 0.: + left = right + right *= factor + # right is now such that cdf(right) > q + + return optimize.brentq(self._ppf_to_solve, \ + left, right, args=(q,)+args, xtol=self.xtol) # moment from definition def _mom_integ0(self, x,m,*args): @@ -1727,7 +1755,6 @@ class rv_continuous(rv_generic): cond2 = (x >= self.b) & cond0 cond = cond0 & cond1 output = zeros(shape(cond),'d') - #place(output,(1-cond0)*(cond1==cond1),self.badvalue) place(output,(1-cond0)+np.isnan(x),self.badvalue) place(output,cond2,1.0) if any(cond): #call only if at least 1 entry @@ -2396,14 +2423,15 @@ class rv_continuous(rv_generic): args = list(args) Nargs = len(args) fixedn = [] - index = range(Nargs) + index = range(Nargs) names = ['f%d' % n for n in range(Nargs - 2)] + ['floc', 'fscale'] - x0 = args[:] - for n, key in zip(index[::-1], names[::-1]): + x0 = [] + for n, key in zip(index, names): if kwds.has_key(key): fixedn.append(n) args[n] = kwds[key] - del x0[n] + else: + x0.append(args[n]) method = kwds.get('method', 'ml').lower() if method.startswith('mps'): fitfun = self.nlogps @@ -2737,12 +2765,7 @@ def _norm_logpdf(x): def _norm_cdf(x): return special.ndtr(x) def _norm_logcdf(x): - logcdf = log(special.ndtr(x)) - large_x = (5 0) & (floor(n)==n) - def _pdf(self, x, n): - Px = (x)**(n-1.0)*exp(-x)/special.gamma(n) + def _rvs(self, a): + return gamma.rvs(a, size=self._size) + def _arg_check(self, a): + return (a > 0) & (floor(a)==a) + def _pdf(self, x, a): + Px = (x)**(a-1.0)*exp(-x)/special.gamma(a) return Px - def _logpdf(self, x, n): - return (n-1.0)*log(x) - x - gamln(n) - def _cdf(self, x, n): - return special.gdtr(1.0,n,x) - def _sf(self, x, n): - return special.gdtrc(1.0,n,x) - def _ppf(self, q, n): - return special.gdtrix(1.0, n, q) - def _stats(self, n): - n = n*1.0 - return n, n, 2/sqrt(n), 6/n - def _entropy(self, n): - return special.psi(n)*(1-n) + 1 + gamln(n) -erlang = erlang_gen(a=0.0, name='erlang', shapes='n') + def _logpdf(self, x, a): + return (a-1.0)*log(x) - x - gamln(a) + def _cdf(self, x, a): + return special.gdtr(1.0,a,x) + def _sf(self, x, a): + return special.gdtrc(1.0,a,x) + def _ppf(self, q, a): + return special.gdtrix(1.0, a, q) + def _stats(self, a): + a = a*1.0 + return a, a, 2/sqrt(a), 6/a + def _entropy(self, a): + return special.psi(a)*(1-a) + 1 + gamln(a) +erlang = erlang_gen(a=0.0, name='erlang', shapes='a') ## Exponential (gamma distributed with a=1.0, loc=loc and scale=scale) @@ -3405,12 +3432,14 @@ class expon_gen(rv_continuous): ----- The probability density function for `expon` is:: - expon.pdf(x) = exp(-x) + expon.pdf(x) = lambda * exp(- lambda*x) for ``x >= 0``. The scale parameter is equal to ``scale = 1.0 / lambda``. + `expon` does not have shape parameters. + %(example)s """ @@ -3601,8 +3630,7 @@ class foldcauchy_gen(rv_continuous): return 1.0/pi*(arctan(x-c) + arctan(x+c)) def _stats(self, c): return inf, inf, nan, nan -# setting xb=1000 allows to calculate ppf for up to q=0.9993 -foldcauchy = foldcauchy_gen(a=0.0, name='foldcauchy', xb=1000, shapes='c') +foldcauchy = foldcauchy_gen(a=0.0, name='foldcauchy', shapes='c') ## F @@ -4218,14 +4246,25 @@ class gamma_gen(rv_continuous): Notes ----- - When ``a`` is an integer, this is the Erlang distribution, and for ``a=1`` - it is the exponential distribution. The probability density function for `gamma` is:: - gamma.pdf(x, a) = x**(a-1) * exp(-x) / gamma(a) + gamma.pdf(x, a) = (lambda*x)**(a-1) * exp(-lambda*x) / gamma(a) + + for ``x >= 0``, ``a > 0``. Here ``gamma(a)`` refers to the gamma function. + + The scale parameter is equal to ``scale = 1.0 / lambda``. + + `gamma` has a shape parameter `a` which needs to be set explicitly. For instance: + + >>> from scipy.stats import gamma + >>> rv = gamma(3., loc = 0., scale = 2.) + + produces a frozen form of `gamma` with shape ``a = 3.``, ``loc = + 0.`` and ``lambda = 1./scale = 1./2.``. - for ``x >= 0``, ``a > 0``. + When ``a`` is an integer, `gamma` reduces to the Erlang + distribution, and when ``a=1`` to the exponential distribution. %(example)s @@ -5590,9 +5629,10 @@ class powerlaw_gen(rv_continuous): def _ppf(self, q, a): return pow(q, 1.0/a) def _stats(self, a): - return a/(a+1.0), a*(a+2.0)/(a+1.0)**2, \ - 2*(1.0-a)*sqrt((a+2.0)/(a*(a+3.0))), \ - 6*polyval([1,-1,-6,2],a)/(a*(a+3.0)*(a+4)) + return (a / (a + 1.0), + a / (a + 2.0) / (a + 1.0) ** 2, + -2.0 * ((a - 1.0) / (a + 3.0)) * sqrt((a + 2.0) / a), + 6 * polyval([1, -1, -6, 2], a) / (a * (a + 3.0) * (a + 4))) def _entropy(self, a): return 1 - 1.0/a - log(a) powerlaw = powerlaw_gen(a=0.0, b=1.0, name="powerlaw", shapes="a") @@ -5894,9 +5934,8 @@ class recipinvgauss_gen(rv_continuous): trm2 = 1.0/mu + x isqx = 1.0/sqrt(x) return 1.0-_norm_cdf(isqx*trm1)-exp(2.0/mu)*_norm_cdf(-isqx*trm2) - # xb=50 or something large is necessary for stats to converge without exception -recipinvgauss = recipinvgauss_gen(a=0.0, xb=50, name='recipinvgauss', - shapes="mu") +recipinvgauss = recipinvgauss_gen(a=0.0, name='recipinvgauss', shapes="mu") + # Semicircular @@ -6096,16 +6135,8 @@ class tukeylambda_gen(rv_continuous): vals2 = log(q/(1-q)) return where((lam == 0)&(q==q), vals2, vals1) def _stats(self, lam): - mu2 = 2*gam(lam+1.5)-lam*pow(4,-lam)*sqrt(pi)*gam(lam)*(1-2*lam) - mu2 /= lam*lam*(1+2*lam)*gam(1+1.5) - mu4 = 3*gam(lam)*gam(lam+0.5)*pow(2,-2*lam) / lam**3 / gam(2*lam+1.5) - mu4 += 2.0/lam**4 / (1+4*lam) - mu4 -= 2*sqrt(3)*gam(lam)*pow(2,-6*lam)*pow(3,3*lam) * \ - gam(lam+1.0/3)*gam(lam+2.0/3) / (lam**3.0 * gam(2*lam+1.5) * \ - gam(lam+0.5)) - g2 = mu4 / mu2 / mu2 - 3.0 - - return 0, mu2, 0, g2 + return 0, _tlvar(lam), 0, _tlkurt(lam) + def _entropy(self, lam): def integ(p): return log(pow(p,lam-1)+pow(1-p,lam-1)) @@ -6118,7 +6149,7 @@ tukeylambda = tukeylambda_gen(name='tukeylambda', shapes="lam") class uniform_gen(rv_continuous): """A uniform continuous random variable. - This distribution is constant between `loc` and ``loc = scale``. + This distribution is constant between `loc` and ``loc + scale``. %(before_notes)s @@ -6263,32 +6294,52 @@ wrapcauchy = wrapcauchy_gen(a=0.0, b=2*pi, name='wrapcauchy', shapes="c") ### DISCRETE DISTRIBUTIONS ### -def entropy(pk,qk=None): - """S = entropy(pk,qk=None) +def entropy(pk, qk=None, base=None): + """ Calculate the entropy of a distribution for given probability values. + + If only probabilities `pk` are given, the entropy is calculated as + ``S = -sum(pk * log(pk), axis=0)``. - calculate the entropy of a distribution given the p_k values - S = -sum(pk * log(pk), axis=0) + If `qk` is not None, then compute a relative entropy + ``S = sum(pk * log(pk / qk), axis=0)``. - If qk is not None, then compute a relative entropy - S = sum(pk * log(pk / qk), axis=0) + This routine will normalize `pk` and `qk` if they don't sum to 1. + + Parameters + ---------- + pk : sequence + Defines the (discrete) distribution. ``pk[i]`` is the (possibly + unnormalized) probability of event ``i``. + qk : sequence, optional + Sequence against which the relative entropy is computed. Should be in + the same format as `pk`. + base : float, optional + The logarithmic base to use, defaults to ``e`` (natural logarithm). + + Returns + ------- + S : float + The calculated entropy. - Routine will normalize pk and qk if they don't sum to 1 """ pk = arr(pk) - pk = 1.0* pk / sum(pk,axis=0) + pk = 1.0* pk / sum(pk, axis=0) if qk is None: vec = where(pk == 0, 0.0, pk*log(pk)) else: qk = arr(qk) if len(qk) != len(pk): raise ValueError("qk and pk must have same length.") - qk = 1.0*qk / sum(qk,axis=0) + qk = 1.0*qk / sum(qk, axis=0) # If qk is zero anywhere, then unless pk is zero at those places # too, the relative entropy is infinite. - if any(take(pk,nonzero(qk==0.0),axis=0)!=0.0, 0): + if any(take(pk, nonzero(qk == 0.0), axis=0) != 0.0, 0): return inf vec = where (pk == 0, 0.0, -pk*log(pk / qk)) - return -sum(vec,axis=0) + S = -sum(vec, axis=0) + if base is not None: + S /= log(base) + return S ## Handlers for generic case where xk and pk are given @@ -6568,6 +6619,8 @@ class rv_discrete(rv_generic): >>> x = np.arange(0, np.min(rv.dist.b, 3)+1) >>> h = plt.plot(x, rv.pmf(x)) + Here, ``rv.dist.b`` is the right endpoint of the support of ``rv.dist``. + Check accuracy of cdf and ppf: >>> prb = generic.cdf(x, ) @@ -6588,8 +6641,8 @@ class rv_discrete(rv_generic): moment_tol=1e-8,values=None,inc=1,longname=None, shapes=None, extradoc=None): - super(rv_generic, self).__init__() - + super(rv_generic,self).__init__() + self.fix_loc = self._fix_loc if badvalue is None: @@ -7517,87 +7570,96 @@ class nbinom_gen(rv_discrete): Notes ----- - Probability mass function, given by - ``np.choose(k+n-1, n-1) * p**n * (1-p)**k`` for ``k >= 0``. + The probability mass function for `nbinom` is:: + + nbinom.pmf(k) = choose(k+n-1, n-1) * p**n * (1-p)**k + + for ``k >= 0``. + + `nbinom` takes ``n`` and ``p`` as shape parameters. %(example)s """ - def _rvs(self, n, pr): - return mtrand.negative_binomial(n, pr, self._size) - def _argcheck(self, n, pr): - return (n >= 0) & (pr >= 0) & (pr <= 1) - def _pmf(self, x, n, pr): + def _rvs(self, n, p): + return mtrand.negative_binomial(n, p, self._size) + def _argcheck(self, n, p): + return (n >= 0) & (p >= 0) & (p <= 1) + def _pmf(self, x, n, p): coeff = exp(gamln(n+x) - gamln(x+1) - gamln(n)) - return coeff * power(pr,n) * power(1-pr,x) - def _logpmf(self, x, n, pr): + return coeff * power(p,n) * power(1-p,x) + def _logpmf(self, x, n, p): coeff = gamln(n+x) - gamln(x+1) - gamln(n) - return coeff + n*log(pr) + x*log1p(-pr) - def _cdf(self, x, n, pr): + return coeff + n*log(p) + x*log1p(-p) + def _cdf(self, x, n, p): k = floor(x) - return special.betainc(n, k+1, pr) - def _sf_skip(self, x, n, pr): + return special.betainc(n, k+1, p) + def _sf_skip(self, x, n, p): #skip because special.nbdtrc doesn't work for 0= q, vals1, vals) - def _stats(self, n, pr): - Q = 1.0 / pr + def _stats(self, n, p): + Q = 1.0 / p P = Q - 1.0 mu = n*P var = n*P*Q g1 = (Q+P)/sqrt(n*P*Q) g2 = (1.0 + 6*P*Q) / (n*P*Q) return mu, var, g1, g2 -nbinom = nbinom_gen(name='nbinom', shapes="n, pr", extradoc=""" - -Negative binomial distribution - -nbinom.pmf(k,n,p) = choose(k+n-1,n-1) * p**n * (1-p)**k -for k >= 0. -""" - ) +nbinom = nbinom_gen(name='nbinom', shapes="n, p") ## Geometric distribution class geom_gen(rv_discrete): - def _rvs(self, pr): - return mtrand.geometric(pr,size=self._size) - def _argcheck(self, pr): - return (pr<=1) & (pr >= 0) - def _pmf(self, k, pr): - return (1-pr)**(k-1) * pr - def _logpmf(self, k, pr): - return (k-1)*log1p(-pr) + pr - def _cdf(self, x, pr): + """A geometric discrete random variable. + + %(before_notes)s + + Notes + ----- + The probability mass function for `geom` is:: + + geom.pmf(k) = (1-p)**(k-1)*p + + for ``k >= 1``. + + `geom` takes ``p`` as shape parameter. + + %(example)s + + """ + def _rvs(self, p): + return mtrand.geometric(p,size=self._size) + def _argcheck(self, p): + return (p<=1) & (p >= 0) + def _pmf(self, k, p): + return (1-p)**(k-1) * p + def _logpmf(self, k, p): + return (k-1)*log1p(-p) + p + def _cdf(self, x, p): k = floor(x) - return (1.0-(1.0-pr)**k) - def _sf(self, x, pr): + return (1.0-(1.0-p)**k) + def _sf(self, x, p): k = floor(x) - return (1.0-pr)**k - def _ppf(self, q, pr): - vals = ceil(log1p(-q)/log1p(-pr)) - temp = 1.0-(1.0-pr)**(vals-1) + return (1.0-p)**k + def _ppf(self, q, p): + vals = ceil(log1p(-q)/log1p(-p)) + temp = 1.0-(1.0-p)**(vals-1) return where((temp >= q) & (vals > 0), vals-1, vals) - def _stats(self, pr): - mu = 1.0/pr - qr = 1.0-pr - var = qr / pr / pr - g1 = (2.0-pr) / sqrt(qr) - g2 = numpy.polyval([1,-6,6],pr)/(1.0-pr) + def _stats(self, p): + mu = 1.0/p + qr = 1.0-p + var = qr / p / p + g1 = (2.0-p) / sqrt(qr) + g2 = numpy.polyval([1,-6,6],p)/(1.0-p) return mu, var, g1, g2 geom = geom_gen(a=1,name='geom', longname="A geometric", - shapes="pr", extradoc=""" + shapes="p") -Geometric distribution - -geom.pmf(k,p) = (1-p)**(k-1)*p -for k >= 1 -""" - ) ## Hypergeometric distribution @@ -7618,7 +7680,37 @@ class hypergeom_gen(rv_discrete): pmf(k, M, n, N) = choose(n, k) * choose(M - n, N - k) / choose(M, N), for N - (M-n) <= k <= min(m,N) - %(example)s + Examples + -------- + >>> from scipy.stats import hypergeom + + Suppose we have a collection of 20 animals, of which 7 are dogs. Then if + we want to know the probability of finding a given number of dogs if we + choose at random 12 of the 20 animals, we can initialize a frozen + distribution and plot the probability mass function: + + >>> [M, n, N] = [20, 7, 12] + >>> rv = hypergeom(M, n, N) + >>> x = np.arange(0, n+1) + >>> pmf_dogs = rv.pmf(x) + + >>> fig = plt.figure() + >>> ax = fig.add_subplot(111) + >>> ax.plot(x, pmf_dogs, 'bo') + >>> ax.vlines(x, 0, pmf_dogs, lw=2) + >>> ax.set_xlabel('# of dogs in our group of chosen animals') + >>> ax.set_ylabel('hypergeom PMF') + >>> plt.show() + + Instead of using a frozen distribution we can also use `hypergeom` + methods directly. To for example obtain the cumulative distribution + function, use: + + >>> prb = hypergeom.cdf(x, M, n, N) + + And to generate random numbers: + + >>> R = hypergeom.rvs(M, n, N, size=10) """ def _rvs(self, M, n, N): @@ -7681,6 +7773,23 @@ hypergeom = hypergeom_gen(name='hypergeom', shapes="M, n, N") ## Logarithmic (Log-Series), (Series) distribution # FIXME: Fails _cdfvec class logser_gen(rv_discrete): + """A Logarithmic (Log-Series, Series) discrete random variable. + + %(before_notes)s + + Notes + ----- + The probability mass function for `logser` is:: + + logser.pmf(k) = - p**k / (k*log(1-p)) + + for ``k >= 1``. + + `logser` takes ``p`` as shape parameter. + + %(example)s + + """ def _rvs(self, pr): # looks wrong for pr>0.5, too few k=1 # trying to use generic is worse, no k=1 at all @@ -7704,23 +7813,35 @@ class logser_gen(rv_discrete): g2 = mu4 / var**2 - 3.0 return mu, var, g1, g2 logser = logser_gen(a=1,name='logser', longname='A logarithmic', - shapes='pr', extradoc=""" - -Logarithmic (Log-Series, Series) distribution - -logser.pmf(k,p) = - p**k / (k*log(1-p)) -for k >= 1 -""" - ) + shapes='p') ## Poisson distribution class poisson_gen(rv_discrete): + """A Poisson discrete random variable. + + %(before_notes)s + + Notes + ----- + The probability mass function for `poisson` is:: + + poisson.pmf(k) = exp(-mu) * mu**k / k! + + for ``k >= 0``. + + `poisson` takes ``mu`` as shape parameter. + + %(example)s + + """ def _rvs(self, mu): return mtrand.poisson(mu, self._size) - def _pmf(self, k, mu): + def _logpmf(self, k, mu): Pk = k*log(mu)-gamln(k+1) - mu - return exp(Pk) + return Pk + def _pmf(self, k, mu): + return exp(self._logpmf(k, mu)) def _cdf(self, x, mu): k = floor(x) return special.pdtr(k,mu) @@ -7737,19 +7858,27 @@ class poisson_gen(rv_discrete): g1 = 1.0/arr(sqrt(mu)) g2 = 1.0 / arr(mu) return mu, var, g1, g2 -poisson = poisson_gen(name="poisson", longname='A Poisson', - shapes="mu", extradoc=""" +poisson = poisson_gen(name="poisson", longname='A Poisson', shapes="mu") -Poisson distribution +## (Planck) Discrete Exponential +class planck_gen(rv_discrete): + """A Planck discrete exponential random variable. -poisson.pmf(k, mu) = exp(-mu) * mu**k / k! -for k >= 0 -""" - ) + %(before_notes)s -## (Planck) Discrete Exponential + Notes + ----- + The probability mass function for `planck` is:: -class planck_gen(rv_discrete): + planck.pmf(k) = (1-exp(-lambda))*exp(-lambda*k) + + for ``k*lambda >= 0``. + + `planck` takes ``lambda`` as shape parameter. + + %(example)s + + """ def _argcheck(self, lambda_): if (lambda_ > 0): self.a = 0 @@ -7781,18 +7910,28 @@ class planck_gen(rv_discrete): l = lambda_ C = -expm1(-l) return l * exp(-l) / C - log(C) -planck = planck_gen(name='planck', longname='A discrete exponential ', - shapes="lamda", - extradoc=""" - -Planck (Discrete Exponential) +planck = planck_gen(name='planck',longname='A discrete exponential ', + shapes="lamda") -planck.pmf(k,b) = (1-exp(-b))*exp(-b*k) -for k*b >= 0 -""" - ) class boltzmann_gen(rv_discrete): + """A Boltzmann (Truncated Discrete Exponential) random variable. + + %(before_notes)s + + Notes + ----- + The probability mass function for `boltzmann` is:: + + boltzmann.pmf(k) = (1-exp(-lambda)*exp(-lambda*k)/(1-exp(-lambda*N)) + + for ``k = 0,...,N-1``. + + `boltzmann` takes ``lambda`` and ``N`` as shape parameters. + + %(example)s + + """ def _pmf(self, k, lambda_, N): fact = (expm1(-lambda_))/(expm1(-lambda_*N)) return fact*exp(-lambda_*k) @@ -7819,22 +7958,28 @@ class boltzmann_gen(rv_discrete): return mu, var, g1, g2 boltzmann = boltzmann_gen(name='boltzmann',longname='A truncated discrete exponential ', - shapes="lamda, N", - extradoc=""" + shapes="lamda, N") -Boltzmann (Truncated Discrete Exponential) +## Discrete Uniform -boltzmann.pmf(k,b,N) = (1-exp(-b))*exp(-b*k)/(1-exp(-b*N)) -for k=0,..,N-1 -""" - ) +class randint_gen(rv_discrete): + """A uniform discrete random variable. + + %(before_notes)s + + Notes + ----- + The probability mass function for `randint` is:: + randint.pmf(k) = 1./(max- min) + for ``k = min,...,max``. + `randint` takes ``min`` and ``max`` as shape parameters. -## Discrete Uniform + %(example)s -class randint_gen(rv_discrete): + """ def _argcheck(self, min, max): self.a = min self.b = max-1 @@ -7868,23 +8013,30 @@ class randint_gen(rv_discrete): def _entropy(self, min, max): return log(max-min) randint = randint_gen(name='randint',longname='A discrete uniform '\ - '(random integer)', shapes="min, max", - extradoc=""" - -Discrete Uniform - - Random integers >=min and = 1``. + + `zipf` takes ``a`` as shape parameter. + + %(example)s + + """ def _rvs(self, a): return mtrand.zipf(a, size=self._size) def _argcheck(self, a): @@ -7910,19 +8062,28 @@ class zipf_gen(rv_discrete): g2 = mu4 / arr(var**2) - 3.0 return mu, var, g1, g2 zipf = zipf_gen(a=1,name='zipf', longname='A Zipf', - shapes="a", extradoc=""" + shapes="a") -Zipf distribution -zipf.pmf(k,a) = 1/(zeta(a)*k**a) -for k >= 1 -""" - ) +# Discrete Laplacian +class dlaplace_gen(rv_discrete): + """A Laplacian discrete random variable. + %(before_notes)s -# Discrete Laplacian + Notes + ----- + The probability mass function for `dlaplace` is:: -class dlaplace_gen(rv_discrete): + dlaplace.pmf(k) = tanh(a/2) * exp(-a*abs(k)) + + for ``a >0``. + + `dlaplace` takes ``a`` as shape parameter. + + %(example)s + + """ def _pmf(self, k, a): return tanh(a/2.0)*exp(-a*abs(k)) def _cdf(self, x, a): @@ -7955,17 +8116,35 @@ class dlaplace_gen(rv_discrete): return a / sinh(a) - log(tanh(a/2.0)) dlaplace = dlaplace_gen(a=-inf, name='dlaplace', longname='A discrete Laplacian', - shapes="a", extradoc=""" + shapes="a") -Discrete Laplacian distribution. -dlaplace.pmf(k,a) = tanh(a/2) * exp(-a*abs(k)) -for a > 0. -""" - ) +class skellam_gen(rv_discrete): + """A Skellam discrete random variable. + %(before_notes)s -class skellam_gen(rv_discrete): + Notes + ----- + Probability distribution of the difference of two correlated or + uncorrelated Poisson random variables. + + Let k1 and k2 be two Poisson-distributed r.v. with expected values + lam1 and lam2. Then, ``k1 - k2`` follows a Skellam distribution with + parameters ``mu1 = lam1 - rho*sqrt(lam1*lam2)`` and + ``mu2 = lam2 - rho*sqrt(lam1*lam2)``, where rho is the correlation + coefficient between k1 and k2. If the two Poisson-distributed r.v. + are independent then ``rho = 0``. + + Parameters mu1 and mu2 must be strictly positive. + + For details see: http://en.wikipedia.org/wiki/Skellam_distribution + + `skellam` takes ``mu1`` and ``mu2`` as shape parameters. + + %(example)s + + """ def _rvs(self, mu1, mu2): n = self._size return np.random.poisson(mu1, n)-np.random.poisson(mu2, n) @@ -7993,26 +8172,9 @@ class skellam_gen(rv_discrete): g2 = 1 / var return mean, var, g1, g2 skellam = skellam_gen(a=-np.inf, name="skellam", longname='A Skellam', - shapes="mu1,mu2", extradoc=""" - -Skellam distribution - - Probability distribution of the difference of two correlated or - uncorrelated Poisson random variables. + shapes="mu1,mu2") - Let k1 and k2 be two Poisson-distributed r.v. with expected values - lam1 and lam2. Then, k1-k2 follows a Skellam distribution with - parameters mu1 = lam1 - rho*sqrt(lam1*lam2) and - mu2 = lam2 - rho*sqrt(lam1*lam2), where rho is the correlation - coefficient between k1 and k2. If the two Poisson-distributed r.v. - are independent then rho = 0. - Parameters mu1 and mu2 must be strictly positive. - - For details see: http://en.wikipedia.org/wiki/Skellam_distribution - -""" - ) def test_lognorm():