#!/usr/bin/env python
# (C) 2000 Huaiyu Zhu <hzhu@users.sourceforge.net>. Licence: GPL
# $Id: distribs.py,v 1.7 2001/08/26 12:40:44 hzhu Exp $
"""
Stats.distribs : Defines statistical distributions
with pdf, cdf, cdfc, cdfi, cdfci, rand.
Two representations: Eg. (sample from Poisson)
<distname>.rand : <distname> could be poisson
rand<dn> : <dn> could be p.
"""
from MatPy.Matrix import to_Matrix,efunc
def rfunc(f, *args): return to_Matrix(apply(f, args))
#------------------------------------------------------------------
import MLab, cephes, RandomArray
_exp = MLab.exp
_log = MLab.log
_sqrt = MLab.sqrt
_bet = cephes.beta
_gam = cephes.gamma
_lgam = cephes.lgam
pi = MLab.pi
_2pi = 2*pi
#_sqrt2 = _sqrt(2)
#==================================================================
class distribution:
"""
Each subclass has pdf, cdf, cdfc, cdfi, cdfci, rand
pdf(x, *a): pdf
cdf(x, *a) : cdf = integral(t=-Inf..x, pdf(t,a,b))
cdfc(x, *a) : complementary cdf = integral(t=x..Inf, pdf(t, *a))
cdfi(y, *a) : x such that cdf(x, a) = y
cdfci(y, *a) : x such that cdfc(x, a) = y
"""
def pdf(self, x): raise NotImplementedError
def cdf(self, x): raise NotImplementedError
def cdfi(self, x): raise NotImplementedError
def cdfc(self, x): return 1. - self.cdf(x)
def cdfci(self, x): return self.cdfi(1.-x)
def rand(self, s, *a): return NotImplementedError
#==================================================================
# Discrete distributions
#------------------------------------------------------------------
# Binomial distribution
def _rand_b(a, b, s): return RandomArray.binomial(a, b, s)
class binomial(distribution):
"""Binomial distribution (n, p):
pdf = nCk p**k (1-p)**(n-k)
"""
def __init__(self, n=1, p=0.5):
self.a = n
self.b = p
#def pdf(self, k): return efunc(_pdf_b, k, self.a, self.b)
def cdf(self, k): return efunc(cephes.bdtr, k, self.a, self.b)
def cdfc(self, k): return efunc(cephes.bdtrc, k, self.a, self.b)
def cdfi(self, y): return efunc(cephes.bdtri, y, self.a, self.b)
def rand(self, s): return rfunc(_rand_b, self.a, self.b, s)
#------------------------------------------------------------------
# Negative binomial distribution
def _rand_nb(a, b, s): return RandomArray.negative_binomial(a, b, s)
class negbinomial(distribution):
"""Negative binomial distribution (n, p):
sum _j=0..k (n+j-1)Cj p**n (1-p)**j,
In a sequence of Bernoulli trials this is the probability that k or
fewer failures precede the nth success.
"""
def __init__(self, n=1, p=0.5):
self.a = n
self.b = p
#def pdf(self, k): return efunc(_pdfnb, k, self.a, self.b)
def cdf(self, k): return efunc(cephes.nbdtr, k, self.a, self.b)
def cdfc(self, k): return efunc(cephes.nbdtrc, k, self.a, self.b)
def cdfi(self, y): return efunc(cephes.nbdtri, y, self.a, self.b)
def rand(self, s): return rfunc(_rand_nb, self.a, self.b, s)
#------------------------------------------------------------------
# Poisson distribution
def _rand_p(a, s): return RandomArray.poisson(a, s)
class poisson(distribution):
"""Binomial distribution (n, p):
sum(exp(-m) * m**j / j!, j=0..k) = igamc( k+1, m).
Arguments must both be positive and k an integer.
"""
def __init__(self, p=1):
self.a = p
#def pdf(self, k): return efunc(_pdf_p, k, self.a)
def cdf(self, k): return efunc(cephes.pdtr, k, self.a)
def cdfc(self, k): return efunc(cephes.pdtrc, k, self.a)
def cdfi(self, y): return efunc(cephes.pdtri, y, self.a)
def rand(self, s): return rfunc(_rand_p, self.a, s)
#==================================================================
# Continuous distributions
#------------------------------------------------------------------
# Beta distribution
def _rand_bt(a, b, s): return RandomArray.beta(a, b, s)
class beta(distribution):
"""Beta distribution (a>0, b>0):
pdf = 1 / beta(a, b)*integral(t**(a-1) (1-t)**(b-1), t=0..x)
"""
def __init__(self, a=1, b=1):
self.a = a
self.b = b
#def pdf(self, x): return efunc(_pdfg, x, self.a, self.b)
def cdf(self, x): return efunc(cephes.btdtr, self.a, self.b, x)
def cdfc(self, x): return efunc(cephes.btdtrc, self.a, self.b, x)
def rand(self, s): return rfunc(_rand_bt, self.a, self.b, s)
#------------------------------------------------------------------
# Gamma distribution
def _pdf_g(x, a, b): return _exp(-_lgam(a) + _log(x*b)*(a-1.) - x*b) * b
def _rand_g(a, b, s): return RandomArray.gamma(a, b, s)
def pdfg(x, a, b):
"pdfg(x,a,b) : gamma pdf = b/gamma(a) * (bx)^(a-1) * exp(-bx)"
return efunc(_pdf_g, x, a, b)
def cdfg(x, a, b):
"cdfg(x,a,b) : gamma cdf = integral(t=0..x, pdfg(x,a,b))"
return efunc(cephes.gdtr, b, a, x)
def cdfcg(x, a, b):
"cdfg(x,a,b) : complemented gamma cdf = integral(t=x..Inf, pdfg(x,a,b))"
return efunc(cephes.gdtrc, b, a, x)
def randg(s, a=1, b=1):
"""randg(s,a,b) : random sample of gamma distribution, size s
pdf = b/gamma(a) * (bx)^(a-1) * exp(-bx)"""
return rfunc(_rand_g, b, a, s)
class gamma(distribution):
"""Gamma distribution (a>0, b>0):
pdf = b/gamma(a) * (b*x)^(a-1) * exp(-b*x)
"""
def __init__(self, a=1, b=1):
self.a = a
self.b = b
def pdf(self, x): return efunc(_pdf_g, x, self.a, self.b)
def cdf(self, x): return efunc(cephes.gdtr, self.b, self.a, x)
def cdfc(self, x): return efunc(cephes.gdtrc, self.b, self.a, x)
def rand(self, s): return rfunc(_rand_g, self.b, self.a, s)
#------------------------------------------------------------------
# Normal distribution
def _pdf_n(x): return _exp(-(x**2/2)) / _sqrt(_2pi)
def _rand_n(s): return RandomArray.standard_normal(s)
def pdfn(x):
"""pdfn(x) : normal pdf = 1/sqrt(2pi) * exp(-(x^2 / 2)"""
return efunc(_pdf_n, x)
def cdfn(x):
"cdfn(x) : normal cdf = integral(t=-Inf..x, pdfn(x))"
return efunc(cephes.ndtr, x)
def cdfcn(x):
"cdfn(x) : complemented normal cdf = integral(t=x..Inf, pdfn(x))"
return 1 - cdfn(x)
def erf(x):
"erf(x) : error function 2/sqrt(pi) * integral(t=0..x, exp(-t**2 / 2))"
return efunc(cephes.erf, x)
def randn(s):
"""randn(s) : random sample of standard normal distribution, size s
"""
return rfunc(_rand_n, s)
class normal(distribution):
"""Normal distribution:
pdf = 1/(b*sqrt(pi)) * exp(-((x-a)/b)^2 / 2)
"""
def pdf(self, x): return efunc(_pdf_n, x)
def cdf(self, x): return efunc(cephes.ndtr, x)
def cdfi(self, x): return efunc(cephes.ndtri, x)
def rand(self, s): return rfunc(_rand_n, s)
#------------------------------------------------------------------
# Chi square distribution
def _pdf_c2(x, n): return _pdf_g(x, n/2., 1/2.)
def _rand_c2(n, s): return RandomArray.chi_square(n, s)
class chi2(distribution):
"""Chi square distribution with n degrees of freedom (n>0):
pdf_t = 1/(2* gam(n/2)) * (t/2)^(n/2-1) * exp(-t/2)
= pdfg(x, n/2, 1/2)
"""
def __init__(self, n=1):
self.a = n
def pdf(self, x): return efunc(_pdf_c2, x, self.a)
def cdf(self, x): return efunc(cephes.chdtr, self.a, x)
def cdfc(self, x): return efunc(cephes.chdtrc, self.a, x)
def cdfi(self, p): return efunc(cephes.chdtri, self.a, p)
def rand(self, s): return rfunc(_rand_c2, self.a, s)
#------------------------------------------------------------------
# Chi distribution
def _rand_c(n, s): return _sqrt(_rand_c2(n,s))
#------------------------------------------------------------------
# Student's t distribution
def _pdf_t(x, n):
n2 = n/2.; n12 = (n + 1)/2.
return _exp(_lgam(n12) - _log(n*pi)/2
- _lgam(n2) + _log(1.+x**2/n) * (-n12))
def _rand_t(n, s):
return _rand_n(s) / _rand_c(n, s)
class t(distribution):
"""Student t distribution (n>0):
pdf = gam((n+1)/2) / (sqrt(n*pi)* gam(n/2)) * (1+x^2/n)^(-(n+1)/2)
"""
def __init__(self, n):
self.a = n
def pdf(self, x): return efunc(_pdf_t, x, self.a)
def cdf(self, x): return efunc(cephes.stdtr, self.a, x)
def cdfi(self, p): return efunc(cephes.stdtri, self.a, p)
def rand(self, s): return rfunc(_rand_t, self.a, s)
#------------------------------------------------------------------
# F distribution
def _pdf_f(x, m, n):
m2 = m/2.; n2 = n/2.;
xmn = 1.*m/n*x;
return _exp(_log(xmn)*(m2-1)
- (m2+n2)*_log(1+xmn)) * (1.*m/n) / _bet(m2,n2)
def _rand_f(m, n, s):
return 1.*n/m * _rand_c2(m, s) / _rand_c2(n, s)
class F(distribution):
"""F distribution (n1>0, n2>0):
(also known as Snedcor's density or the variance ratio density)
pdf = density of (u1/n1)/(u2/n2),
where u1 and u2 are Chi variates with n1 and n2 degrees of freedom.
"""
def __init__(self, n1=1, n2=1):
self.a = n1
self.b = n2
def pdf(self, x): return efunc(_pdf_f, x, self.a, self.b)
def cdf(self, x): return efunc(cephes.fdtr, self.a, self.b, x)
def cdfc(self, x): return efunc(cephes.fdtrc, self.a, self.b, x)
def cdfi(self, p): return efunc(cephes.fdtri, self.a, self.b, p)
def rand(self, s): return rfunc(_rand_f, self.a, self.b, s)
|