1
# Functions to implement several important functions for
2
# various Continous and Discrete Probability Distributions
4
# Author: Travis Oliphant 2002-2003
7
from __future__ import nested_scopes
9
import scipy.special as special
10
import scipy.optimize as optimize
13
from Numeric import alltrue, where, arange, put, putmask, nonzero, \
14
ravel, compress, take, ones, sum, shape, product, repeat, reshape, \
16
from scipy_base.fastumath import *
17
from scipy_base import atleast_1d, polyval, angle, ceil, insert, extract, \
18
any, argsort, argmax, argmin, vectorize, r_, asarray, nan, inf, select
21
errp = special.errprint
36
"""seed(x, y), set the seed using the integers x, y;
37
Set a random one from clock if y == 0
39
if type (x) != types.IntType or type (y) != types.IntType :
40
raise ArgumentError, "seed requires integer arguments."
43
y = int(rv.initial_seed())
44
x = random.randint(1,int(2l**31-2)) # Python 2.1 and 2.2 require this
45
# bit of hackery to escape some
46
# wierd int/longint conversion
53
"Return the current seed pair"
54
return rand.get_seeds()
56
def _build_random_array(fun, args, size=None):
57
# Build an array by applying function fun to
58
# the arguments in args, creating an array with
59
# the specified shape.
60
# Allows an integer shape n as a shorthand for (n,).
61
if isinstance(size, types.IntType):
63
if size is not None and len(size) != 0:
64
n = Num.multiply.reduce(size)
65
s = apply(fun, args + (n,))
70
s = apply(fun, args + (n,))
73
def random(size=None):
74
"Returns array of random numbers between 0 and 1"
75
return _build_random_array(rand.sample, (), size)
77
def random_integers(max, min=1, size=None):
78
"""random integers in range min to max inclusive"""
79
return randint.rvs(min, max+1, size)
82
"""If arg is an integer, a permutation of indices arange(n), otherwise
83
a permuation of the sequence"""
84
if isinstance(arg,types.IntType):
86
return rand.permutation(arg)
89
## Internal class to compute a ppf given a distribution.
90
## (needs cdf function) and uses brentq from scipy.optimize
91
## to compute ppf from cdf.
92
class general_cont_ppf:
93
def __init__(self, dist, xa=-10.0, xb=10.0, xtol=1e-14):
95
self.cdf = eval('%scdf'%dist)
99
self.vecfunc = sgf(self._single_call,otypes='d')
100
def _tosolve(self, x, q, *args):
101
return apply(self.cdf, (x, )+args) - q
102
def _single_call(self, q, *args):
103
return scipy.optimize.brentq(self._tosolve, self.xa, self.xb, args=(q,)+args, xtol=self.xtol)
104
def __call__(self, q, *args):
105
return self.vecfunc(q, *args)
110
def __init__(self, dist, *args, **kwds):
115
return self.dist.pdf(x,*self.args,**self.kwds)
117
return self.dist.cdf(x,*self.args,**self.kwds)
119
return self.dist.ppf(q,*self.args,**self.kwds)
121
return self.dist.isf(q,*self.args,**self.kwds)
122
def rvs(self, size=None):
124
kwds.update({'size':size})
125
return self.dist.rvs(*self.args,**kwds)
127
return self.dist.sf(x,*self.args,**self.kwds)
129
return self.dist.stats(*self.args,**self.kwds)
133
## NANs are returned for unsupported parameters.
134
## location and scale parameters are optional for each distribution.
135
## The shape parameters are generally required
137
## The loc and scale parameters must be given as keyword parameters.
138
## These are related to the common symbols in the .lyx file
140
## skew is third central moment / variance**(1.5)
141
## kurtosis is fourth central moment / variance**2 - 3
146
## Documentation for ranlib, rv2, cdflib and
148
## Eric Wesstein's world of mathematics http://mathworld.wolfram.com/
149
## http://mathworld.wolfram.com/topics/StatisticalDistributions.html
151
## Documentation to Regress+ by Michael McLaughlin
153
## Engineering and Statistics Handbook (NIST)
154
## http://www.itl.nist.gov/div898/handbook/index.htm
156
## Documentation for DATAPLOT from NIST
157
## http://www.itl.nist.gov/div898/software/dataplot/distribu.htm
159
## Norman Johnson, Samuel Kotz, and N. Balakrishnan "Continuous
160
## Univariate Distributions", second edition,
161
## Volumes I and II, Wiley & Sons, 1994.
164
## Each continuous random variable as the following methods
166
## rvs -- Random Variates (alternatively calling the class could produce these)
169
## sf -- Survival Function (1-CDF)
170
## ppf -- Percent Point Function (Inverse of CDF)
171
## isf -- Inverse Survival Function (Inverse of SF)
172
## stats -- Return mean, variance, (Fisher's) skew, or (Fisher's) kurtosis
173
## nnlf -- negative log likelihood function (to minimize)
174
## fit -- Model-fitting
178
## hf --- Hazard Function (PDF / SF)
179
## chf --- Cumulative hazard function (-log(SF))
180
## psf --- Probability sparsity function (reciprocal of the pdf) in
181
## units of percent-point-function (as a function of q).
182
## Also, the derivative of the percent-point function.
184
## To define a new random variable you subclass the rv_continuous class
187
## _pdf method which will be given clean arguments (in between a and b)
188
## and passing the argument check method
190
## If postive argument checking is not correct for your RV
191
## then you will also need to re-define
194
## Correct, but potentially slow defaults exist for the remaining
195
## methods but for speed and/or accuracy you can over-ride
197
## _cdf, _ppf, _rvs, _isf, _sf
199
## Rarely would you override _isf and _sf but you could.
201
## Statistics are computed using numerical integration by default.
202
## For speed you can redefine this using
204
## _stats --- take shape parameters and return mu, mu2, g1, g2
205
## --- If you can't compute one of these return it as None
207
## --- Can also be defined with a keyword argument moments=<str>
208
## where <str> is a string composed of 'm', 'v', 's',
209
## and/or 'k'. Only the components appearing in string
210
## should be computed and returned in the order 'm', 'v',
211
## 's', or 'k' with missing values returned as None
217
## _munp -- takes n and shape parameters and returns
218
## -- then nth non-central moment of the distribution.
221
def valarray(shape,value=nan,typecode=None):
222
"""Return an array of all value.
224
out = reshape(repeat([value],product(shape)),shape)
228
return out.astype(typecode)
230
def argsreduce(cond, *args):
231
"""Return a sequence of arguments converted to the dimensions of cond
234
expand_arr = (cond==cond)
235
for k in range(len(args)):
236
newargs[k] = extract(cond,arr(args[k])*expand_arr)
240
"""A Generic continuous random variable.
242
Continuous random variables are defined from a standard form chosen
243
for simplicity of representation. The standard form may require
244
some shape parameters to complete its specification. The distributions
245
also take optional location and scale parameters using loc= and scale=
246
keywords (defaults: loc=0, scale=1)
248
These shape, scale, and location parameters can be passed to any of the
249
methods of the RV object such as the following:
251
generic.rvs(<shape(s)>,loc=0,scale=1)
254
generic.pdf(x,<shape(s)>,loc=0,scale=1)
255
- probability density function
257
generic.cdf(x,<shape(s)>,loc=0,scale=1)
258
- cumulative density function
260
generic.sf(x,<shape(s)>,loc=0,scale=1)
261
- survival function (1-cdf --- sometimes more accurate)
263
generic.ppf(q,<shape(s)>,loc=0,scale=1)
264
- percent point function (inverse of cdf --- percentiles)
266
generic.isf(q,<shape(s)>,loc=0,scale=1)
267
- inverse survival function (inverse of sf)
269
generic.stats(<shape(s)>,loc=0,scale=1,moments='mv')
270
- mean('m'), variance('v'), skew('s'), and/or kurtosis('k')
272
generic.entropy(<shape(s)>,loc=0,scale=1)
273
- (differential) entropy of the RV.
275
Alternatively, the object may be called (as a function) to fix
276
the shape, location, and scale parameters returning a
277
"frozen" continuous RV object:
279
myrv = generic(<shape(s)>,loc=0,scale=1)
280
- frozen RV object with the same methods but holding the
281
given shape, location, and scale fixed
283
def __init__(self, momtype=1, a=None, b=None, xa=-10.0, xb=10.0, xtol=1e-14, badvalue=None, name=None, longname=None, shapes=None, extradoc=None):
286
self.badvalue = badvalue
299
self.moment_type = momtype
300
self.vecfunc = sgf(self._ppf_single_call,otypes='d')
301
self.vecentropy = sgf(self._entropy,otypes='d')
302
self.veccdf = sgf(self._cdf_single_call,otypes='d')
305
self.generic_moment = sgf(self._mom0_sc,otypes='d')
307
self.generic_moment = sgf(self._mom1_sc,otypes='d')
308
cdf_signature = inspect.getargspec(self._cdf.im_func)
309
numargs1 = len(cdf_signature[0]) - 2
310
pdf_signature = inspect.getargspec(self._pdf.im_func)
311
numargs2 = len(pdf_signature[0]) - 2
312
self.numargs = max(numargs1, numargs2)
315
if name[0] in ['aeiouAEIOU']: hstr = "An "
317
longname = hstr + name
318
if self.__doc__ is None:
319
self.__doc__ = rv_continuous.__doc__
320
if self.__doc__ is not None:
321
self.__doc__ = self.__doc__.replace("A Generic",longname)
323
self.__doc__ = self.__doc__.replace("generic",name)
325
self.__doc__ = self.__doc__.replace("<shape(s)>,","")
327
self.__doc__ = self.__doc__.replace("<shape(s)>",shapes)
328
if extradoc is not None:
329
self.__doc__ = self.__doc__ + extradoc
331
def _ppf_to_solve(self, x, q,*args):
332
return apply(self.cdf, (x, )+args)-q
334
def _ppf_single_call(self, q, *args):
335
return scipy.optimize.brentq(self._ppf_to_solve, self.xa, self.xb, args=(q,)+args, xtol=self.xtol)
337
# moment from definition
338
def _mom_integ0(self, x,m,*args):
339
return x**m * self.pdf(x,*args)
340
def _mom0_sc(self, m,*args):
341
return scipy.integrate.quad(self._mom_integ0, self.a,
342
self.b, args=(m,)+args)[0]
343
# moment calculated using ppf
344
def _mom_integ1(self, q,m,*args):
345
return (self.ppf(q,*args))**m
346
def _mom1_sc(self, m,*args):
347
return scipy.integrate.quad(self._mom_integ1, 0, 1,args=(m,)+args)[0]
349
## These are the methods you must define (standard form functions)
350
def _argcheck(self, *args):
351
# Default check for correct values on args and keywords.
352
# Returns condition array of 1's where arguments are correct and
353
# 0's where they are not.
356
cond = logical_and(cond,(arr(arg) > 0))
359
def _pdf(self,x,*args):
360
return scipy.derivative(self._cdf,x,dx=1e-5,args=args,order=5)
362
## Could also define any of these (return 1-d using self._size to get number)
363
def _rvs(self, *args):
364
## Use basic inverse cdf algorithm for RV generation as default.
365
U = rand.sample(self._size)
366
Y = self._ppf(U,*args)
369
def _cdf_single_call(self, x, *args):
370
return scipy.integrate.quad(self._pdf, self.a, x, args=args)[0]
372
def _cdf(self, x, *args):
373
return self.veccdf(x,*args)
375
def _sf(self, x, *args):
376
return 1.0-self._cdf(x,*args)
378
def _ppf(self, q, *args):
379
return self.vecfunc(q,*args)
381
def _isf(self, q, *args):
382
return self.vecfunc(1.0-q,*args)
384
# The actual cacluation functions (no basic checking need be done)
385
# If these are defined, the others won't be looked at.
386
# Otherwise, the other set can be defined.
387
def _stats(self,*args, **kwds):
388
moments = kwds.get('moments')
389
return None, None, None, None
392
def _munp(self,n,*args):
393
return self.generic_moment(n,*args)
395
def __fix_loc_scale(self, args, loc, scale):
398
if N == self.numargs + 1 and loc is None: # loc is given without keyword
400
if N == self.numargs + 2 and scale is None: # loc and scale given without keyword
401
loc, scale = args[-2:]
402
args = args[:self.numargs]
407
return args, loc, scale
409
# These are actually called, but should not
410
# be overwritten if you want to keep
411
# the error checking.
412
def rvs(self,*args,**kwds):
413
"""Random variates of given type.
417
The shape parameter(s) for the distribution (see docstring of the
418
instance object for more information)
422
size - number of random variates (default=1)
423
loc - location parameter (default=0)
424
scale - scale parameter (default=1)
426
loc,scale,size=map(kwds.get,['loc','scale','size'])
427
args, loc, scale = self.__fix_loc_scale(args, loc, scale)
428
cond = logical_and(self._argcheck(*args),(scale >= 0))
430
raise ValueError, "Domain error in arguments."
435
self._size = product(size)
436
if scipy.isscalar(size):
440
vals = reshape(self._rvs(*args),size)
442
return loc*ones(size,'d')
444
return vals * scale + loc
446
def pdf(self,x,*args,**kwds):
447
"""Probability density function at x of the given RV.
451
The shape parameter(s) for the distribution (see docstring of the
452
instance object for more information)
456
loc - location parameter (default=0)
457
scale - scale parameter (default=1)
459
loc,scale=map(kwds.get,['loc','scale'])
460
args, loc, scale = self.__fix_loc_scale(args, loc, scale)
461
x,loc,scale = map(arr,(x,loc,scale))
462
args = tuple(map(arr,args))
463
x = arr((x-loc)*1.0/scale)
464
cond0 = self._argcheck(*args) & (scale > 0)
465
cond1 = (scale > 0) & (x > self.a) & (x < self.b)
467
output = zeros(shape(cond),'d')
468
insert(output,(1-cond0)*(cond1==cond1),self.badvalue)
469
goodargs = argsreduce(cond, *((x,)+args+(scale,)))
470
scale, goodargs = goodargs[-1], goodargs[:-1]
471
insert(output,cond,self._pdf(*goodargs) / scale)
474
def cdf(self,x,*args,**kwds):
475
"""Cumulative distribution function at x of the given RV.
479
The shape parameter(s) for the distribution (see docstring of the
480
instance object for more information)
484
loc - location parameter (default=0)
485
scale - scale parameter (default=1)
487
loc,scale=map(kwds.get,['loc','scale'])
488
args, loc, scale = self.__fix_loc_scale(args, loc, scale)
489
x,loc,scale = map(arr,(x,loc,scale))
490
args = tuple(map(arr,args))
491
x = (x-loc)*1.0/scale
492
cond0 = self._argcheck(*args) & (scale > 0)
493
cond1 = (scale > 0) & (x > self.a) & (x < self.b)
494
cond2 = (x >= self.b) & cond0
496
output = zeros(shape(cond),'d')
497
insert(output,(1-cond0)*(cond1==cond1),self.badvalue)
498
insert(output,cond2,1.0)
499
goodargs = argsreduce(cond, *((x,)+args))
500
insert(output,cond,self._cdf(*goodargs))
503
def sf(self,x,*args,**kwds):
504
"""Survival function (1-cdf) at x of the given RV.
508
The shape parameter(s) for the distribution (see docstring of the
509
instance object for more information)
513
loc - location parameter (default=0)
514
scale - scale parameter (default=1)
516
loc,scale=map(kwds.get,['loc','scale'])
517
args, loc, scale = self.__fix_loc_scale(args, loc, scale)
518
x,loc,scale = map(arr,(x,loc,scale))
519
args = tuple(map(arr,args))
520
x = (x-loc)*1.0/scale
521
cond0 = self._argcheck(*args) & (scale > 0)
522
cond1 = (scale > 0) & (x > self.a) & (x < self.b)
523
cond2 = cond0 & (x <= self.a)
525
output = zeros(shape(cond),'d')
526
insert(output,(1-cond0)*(cond1==cond1),self.badvalue)
527
insert(output,cond2,1.0)
528
goodargs = argsreduce(cond, *((x,)+args))
529
insert(output,cond,self._sf(*goodargs))
532
def ppf(self,q,*args,**kwds):
533
"""Percent point function (inverse of cdf) at q of the given RV.
537
The shape parameter(s) for the distribution (see docstring of the
538
instance object for more information)
542
loc - location parameter (default=0)
543
scale - scale parameter (default=1)
545
loc,scale=map(kwds.get,['loc','scale'])
546
args, loc, scale = self.__fix_loc_scale(args, loc, scale)
547
q,loc,scale = map(arr,(q,loc,scale))
548
args = tuple(map(arr,args))
549
cond0 = self._argcheck(*args) & (scale > 0) & (loc==loc)
550
cond1 = (q > 0) & (q < 1)
551
cond2 = (q==1) & cond0
553
output = valarray(shape(cond),value=self.a*scale + loc)
554
insert(output,(1-cond0)+(1-cond1)*(q!=0.0), self.badvalue)
555
insert(output,cond2,self.b*scale + loc)
556
goodargs = argsreduce(cond, *((q,)+args+(scale,loc)))
557
scale, loc, goodargs = goodargs[-2], goodargs[-1], goodargs[:-2]
558
insert(output,cond,self._ppf(*goodargs)*scale + loc)
561
def isf(self,q,*args,**kwds):
562
"""Inverse survival function at q of the given RV.
566
The shape parameter(s) for the distribution (see docstring of the
567
instance object for more information)
571
loc - location parameter (default=0)
572
scale - scale parameter (default=1)
574
loc,scale=map(kwds.get,['loc','scale'])
575
args, loc, scale = self.__fix_loc_scale(args, loc, scale)
576
q,loc,scale = map(arr,(q,loc,scale))
577
args = tuple(map(arr,args))
578
cond0 = self._argcheck(*args) & (scale > 0) & (loc==loc)
579
cond1 = (q > 0) & (q < 1)
580
cond2 = (q==1) & cond0
582
output = valarray(shape(cond),value=self.b)
583
insert(output,(1-cond0)*(cond1==cond1), self.badvalue)
584
insert(output,cond2,self.a)
585
goodargs = argsreduce(cond, *((1.0-q,)+args+(scale,loc)))
586
scale, loc, goodargs = goodargs[-2], goodargs[-1], goodargs[:-2]
587
insert(output,cond,self._ppf(*goodargs)*scale + loc)
590
def stats(self,*args,**kwds):
591
"""Some statistics of the given RV
595
The shape parameter(s) for the distribution (see docstring of the
596
instance object for more information)
600
loc - location parameter (default=0)
601
scale - scale parameter (default=1)
602
moments - a string composed of letters ['mvsk'] specifying
603
which moments to compute (default='mv')
606
's' = (Fisher's) skew,
607
'k' = (Fisher's) kurtosis.
609
loc,scale,moments=map(kwds.get,['loc','scale','moments'])
613
if N == self.numargs + 1 and loc is None: # loc is given without keyword
615
if N == self.numargs + 2 and scale is None: # loc and scale given without keyword
616
loc, scale = args[-2:]
617
if N == self.numargs + 3 and moments is None: # loc, scale, and moments
618
loc, scale, moments = args[-3:]
619
args = args[:self.numargs]
620
if scale is None: scale = 1.0
621
if loc is None: loc = 0.0
622
if moments is None: moments = 'mv'
624
loc,scale = map(arr,(loc,scale))
625
args = tuple(map(arr,args))
626
cond = self._argcheck(*args) & (scale > 0) & (loc==loc)
628
signature = inspect.getargspec(self._stats.im_func)
629
if (signature[2] is not None) or ('moments' in signature[0]):
630
mu, mu2, g1, g2 = self._stats(*args,**{'moments':moments})
632
mu, mu2, g1, g2 = self._stats(*args)
637
default = valarray(shape(cond), self.badvalue)
640
# Use only entries that are valid in calculation
641
goodargs = argsreduce(cond, *(args+(scale,loc)))
642
scale, loc, goodargs = goodargs[-2], goodargs[-1], goodargs[:-2]
646
mu = self._munp(1.0,*goodargs)
647
out0 = default.copy()
648
insert(out0,cond,mu*scale+loc)
653
mu2p = self._munp(2.0,*goodargs)
655
mu = self._munp(1.0,*goodargs)
657
out0 = default.copy()
658
insert(out0,cond,mu2*scale*scale)
663
mu3p = self._munp(3.0,*goodargs)
665
mu = self._munp(1.0,*goodargs)
667
mu2p = self._munp(2.0,*goodargs)
669
mu3 = mu3p - 3*mu*mu2 - mu**3
671
out0 = default.copy()
677
mu4p = self._munp(4.0,*goodargs)
679
mu = self._munp(1.0,*goodargs)
681
mu2p = self._munp(2.0,*goodargs)
684
mu3p = self._munp(3.0,*goodargs)
685
mu3 = mu3p - 3*mu*mu2 - mu**3
686
mu4 = mu4p - 4*mu*mu3 - 6*mu*mu*mu2 - mu**4
687
g2 = mu4 / mu2**2.0 - 3.0
688
out0 = default.copy()
697
def moment(self, n, *args):
699
raise ValueError, "Moment must be an integer."
700
if (n < 0): raise ValueError, "Moment must be positive."
701
if (n == 0): return 1.0
702
if (n > 0) and (n < 5):
703
signature = inspect.getargspec(self._stats.im_func)
704
if (signature[2] is not None) or ('moments' in signature[0]):
705
dict = {'moments':{1:'m',2:'v',3:'vs',4:'vk'}[n]}
708
mu, mu2, g1, g2 = self._stats(*args,**dict)
710
if mu is None: return self._munp(1,*args)
713
if mu2 is None: return self._munp(2,*args)
716
if g1 is None or mu2 is None: return self._munp(3,*args)
717
else: return g1*(mu2**1.5)
719
if g2 is None or mu2 is None: return self._munp(4,*args)
720
else: return (g2+3.0)*(mu2**2.0)
722
return self._munp(n,*args)
724
def _nnlf(self, x, *args):
725
return -sum(log(self._pdf(x, *args)))
727
def nnlf(self, *args):
728
# - sum (log pdf(x, theta))
729
# where theta are the parameters (including loc and scale)
737
raise ValueError, "Not enough input arguments."
738
if not self._argcheck(*args) or scale <= 0:
740
x = arr((x-loc) / scale)
741
cond0 = (x <= self.a) | (x >= self.b)
746
return self._nnlf(self, x, *args) + N*log(scale)
748
def fit(self, data, *args, **kwds):
749
loc0, scale0 = map(kwds.get, ['loc', 'scale'],[0.0, 1.0])
751
if Narg != self.numargs:
752
if Narg > self.numargs:
753
raise ValueError, "Too many input arguments."
755
args += (1.0,)*(self.numargs-Narg)
756
# location and scale are at the end
757
x0 = args + (loc0, scale0)
758
return optimize.fmin(self.nnlf,x0,args=(ravel(data),),disp=0)
760
def est_loc_scale(self, data, *args):
761
mu, mu2, g1, g2 = self.stats(*args,**{'moments':'mv'})
762
muhat = stats.nanmean(data)
763
mu2hat = stats.nanstd(data)
764
Shat = sqrt(mu2hat / mu2)
765
Lhat = muhat - Shat*mu
768
def freeze(self,*args,**kwds):
769
return rv_frozen(self,*args,**kwds)
771
def __call__(self, *args, **kwds):
772
return self.freeze(*args, **kwds)
774
def _entropy(self, *args):
776
val = self._pdf(x, *args)
778
return -scipy.integrate.quad(integ,self.a,self.b)[0]
780
def entropy(self, *args, **kwds):
781
loc,scale=map(kwds.get,['loc','scale'])
782
args, loc, scale = self.__fix_loc_scale(args, loc, scale)
784
cond0 = self._argcheck(*args) & (scale > 0) & (loc==loc)
785
output = zeros(shape(cond0),'d')
786
insert(output,(1-cond0),self.badvalue)
787
goodargs = argsreduce(cond0, *args)
788
insert(output,cond0,self.vecentropy(*goodargs)+log(scale))
791
_EULER = 0.577215664901532860606512090082402431042 # -special.psi(1)
792
_ZETA3 = 1.202056903159594285399738161511449990765 # special.zeta(3,1) Apery's constant
794
## Kolmogorov-Smirnov one-sided and two-sided test statistics
796
class ksone_gen(rv_continuous):
798
return 1.0-special.smirnov(n,x)
800
return special.smirnovi(n,1.0-q)
801
ksone = ksone_gen(a=0.0,name='ksone', longname="Kolmogorov-Smirnov "\
802
"A one-sided test statistic.", shapes="n",
805
General Kolmogorov-Smirnov one-sided test.
809
class kstwobign_gen(rv_continuous):
811
return 1.0-special.kolmogorov(x)
813
return special.kolmogi(1.0-q)
814
kstwobign = kstwobign_gen(a=0.0,name='kstwobign', longname='Kolmogorov-Smirnov two-sided (for large N)', extradoc="""
816
Kolmogorov-Smirnov two-sided test for large N
821
## Normal distribution
823
# loc = mu, scale = std
824
class norm_gen(rv_continuous):
826
return rand.standard_normal(self._size)
828
return 1.0/sqrt(2*pi)*exp(-x**2/2.0)
830
return special.ndtr(x)
832
return special.ndtri(q)
834
return 0.0, 1.0, 0.0, 0.0
836
return 0.5*(log(2*pi)+1)
837
norm = norm_gen(name='norm',longname='A normal',extradoc="""
841
The location (loc) keyword specifies the mean.
842
The scale (scale) keyword specifies the standard deviation.
845
## Alpha distribution
847
class alpha_gen(rv_continuous):
848
def _pdf(self, x, a):
849
return 1.0/arr(x**2)/special.ndtr(a)*norm.pdf(a-1.0/x)
850
def _cdf(self, x, a):
851
return special.ndtr(a-1.0/x) / special.ndtr(a)
852
def _ppf(self, q, a):
853
return 1.0/arr(a-special.ndtri(q*special.ndtr(a)))
855
return [inf]*2 + [nan]*2
856
alpha = alpha_gen(a=0.0,name='alpha',shapes='a',extradoc="""
861
## Anglit distribution
863
class anglit_gen(rv_continuous):
867
return sin(x+pi/4)**2.0
869
return (arcsin(sqrt(q))-pi/4)
871
return 0.0, pi*pi/16-0.5, 0.0, -2*(pi**4 - 96)/(pi*pi-8)**2
874
anglit = anglit_gen(a=-pi/4,b=pi/4,name='anglit', extradoc="""
881
## Arcsine distribution
883
class arcsine_gen(rv_continuous):
885
return 1.0/pi/sqrt(x*(1-x))
887
return 2.0/pi*arcsin(sqrt(x))
889
return sin(pi/2.0*q)**2.0
891
mup = 0.5, 3.0/8.0, 15.0/48.0, 35.0/128.0
896
return mu, mu2, g1, g2
898
return -0.24156447527049044468
899
arcsine = arcsine_gen(a=0.0,b=1.0,name='arcsine',extradoc="""
907
class beta_gen(rv_continuous):
908
def _rvs(self, a, b):
909
return rand.beta(a,b,self._size)
910
def _pdf(self, x, a, b):
911
Px = (1.0-x)**(b-1.0) * x**(a-1.0)
912
Px /= special.beta(a,b)
914
def _cdf(self, x, a, b):
915
return special.btdtr(a,b,x)
916
def _ppf(self, q, a, b):
917
return special.btdtri(a,b,q)
918
def _stats(self, a, b):
919
mn = a *1.0 / (a + b)
920
var = (a*b*1.0)*(a+b+1.0)/(a+b)**2.0
921
g1 = 2.0*(b-a)*sqrt((1.0+a+b)/(a*b)) / (2+a+b)
922
g2 = 6.0*(a**3 + a**2*(1-2*b) + b**2*(1+b) - 2*a*b*(2+b))
923
g2 /= a*b*(a+b+2)*(a+b+3)
924
return mn, var, g1, g2
925
beta = beta_gen(a=0.0, b=1.0, name='beta',shapes='a,b',extradoc="""
931
class betaprime_gen(rv_continuous):
932
def _rvs(self, a, b):
933
u1 = gamma.rvs(a,size=self._size)
934
u2 = gamma.rvs(b,size=self._size)
936
def _pdf(self, x, a, b):
937
return 1.0/special.beta(a,b)*x**(a-1.0)/(1+x)**(a+b)
938
def _cdf(self, x, a, b):
939
x = where(x==1.0, 1.0-1e-6,x)
940
return pow(x,a)*special.hyp2f1(a+b,a,1+a,-x)/a/special.beta(a,b)
941
def _munp(self, n, a, b):
943
return where(b > 1, a/(b-1.0), inf)
945
return where(b > 2, a*(a+1.0)/((b-2.0)*(b-1.0)), inf)
947
return where(b > 3, a*(a+1.0)*(a+2.0)/((b-3.0)*(b-2.0)*(b-1.0)),
951
a*(a+1.0)*(a+2.0)*(a+3.0)/((b-4.0)*(b-3.0) \
952
*(b-2.0)*(b-1.0)), inf)
954
raise NotImplementedError
955
betaprime = betaprime_gen(a=0.0, b=500.0, name='betaprime', shapes='a,b',
958
Beta prime distribution
964
class bradford_gen(rv_continuous):
965
def _pdf(self, x, c):
966
return c / (c*x + 1.0) / log(1.0+c)
967
def _cdf(self, x, c):
968
return log(1.0+c*x) / log(c+1.0)
969
def _ppf(self, q, c):
970
return ((1.0+c)**q-1)/c
971
def _stats(self, c, moments='mv'):
974
mu2 = ((c+2.0)*k-2.0*c)/(2*c*k*k)
978
g1 = sqrt(2)*(12*c*c-9*c*k*(c+2)+2*k*k*(c*(c+3)+3))
979
g1 /= sqrt(c*(c*(k-2)+2*k))*(3*c*(k-2)+6*k)
981
g2 = c**3*(k-3)*(k*(3*k-16)+24)+12*k*c*c*(k-4)*(k-3) \
982
+ 6*c*k*k*(3*k-14) + 12*k**3
983
g2 /= 3*c*(c*(k-2)+2*k)**2
984
return mu, mu2, g1, g2
985
def _entropy(self, c):
987
return k/2.0 - log(c/k)
988
bradford = bradford_gen(a=0.0, b=1.0, name='bradford', longname="A Bradford",
989
shapes='c', extradoc="""
991
Bradford distribution
997
# burr with d=1 is called the fisk distribution
998
class burr_gen(rv_continuous):
999
def _pdf(self, x, c, d):
1000
return c*d*(x**(-c-1.0))*((1+x**(-c*1.0))**(-d-1.0))
1001
def _cdf(self, x, c, d):
1002
return (1+x**(-c*1.0))**(-d**1.0)
1003
def _ppf(self, q, c, d):
1004
return (q**(-1.0/d)-1)**(-1.0/c)
1005
def _stats(self, c, d, moments='mv'):
1006
g2c, g2cd = gam(1-2.0/c), gam(2.0/c+d)
1007
g1c, g1cd = gam(1-1.0/c), gam(1.0/c+d)
1009
k = gd*g2c*g2cd - g1c**2 * g1cd**2
1013
g3c, g3cd = None, None
1015
g3c, g3cd = gam(1-3.0/c), gam(3.0/c+d)
1016
g1 = 2*g1c**3 * g1cd**3 + gd*gd*g3c*g3cd - 3*gd*g2c*g1c*g1cd*g2cd
1023
g4c, g4cd = gam(1-4.0/c), gam(4.0/c+d)
1024
g2 = 6*gd*g2c*g2cd * g1c**2 * g1cd**2 + gd**3 * g4c*g4cd
1025
g2 -= 3*g1c**4 * g1cd**4 -4*gd**2*g3c*g1c*g1cd*g3cd
1026
return mu, mu2, g1, g2
1027
burr = burr_gen(a=0.0, name='burr', longname="Burr",
1028
shapes="c,d", extradoc="""
1034
# burr is a generalization
1036
class fisk_gen(burr_gen):
1037
def _pdf(self, x, c):
1038
return burr_gen._pdf(self, x, c, 1.0)
1039
def _cdf(self, x, c):
1040
return burr_gen._cdf(self, x, c, 1.0)
1041
def _ppf(self, x, c):
1042
return burr_gen._ppf(self, x, c, 1.0)
1043
def _stats(self, c):
1044
return burr_gen._stats(self, c, 1.0)
1045
def _entropy(self, c):
1047
fisk = fisk_gen(a=0.0, name='fink', longname="A funk",
1048
shapes='c', extradoc="""
1058
class cauchy_gen(rv_continuous):
1060
return 1.0/pi/(1.0+x*x)
1062
return 0.5 + 1.0/pi*arctan(x)
1064
return tan(pi*q-pi/2.0)
1066
return 0.5 - 1.0/pi*arctan(x)
1068
return tan(pi/2.0-pi*q)
1070
return inf, inf, nan, nan
1073
cauchy = cauchy_gen(name='cauchy',longname='Cauchy',extradoc="""
1080
## (positive square-root of chi-square)
1081
## chi(1, loc, scale) = halfnormal
1082
## chi(2, 0, scale) = Rayleigh
1083
## chi(3, 0, scale) = MaxWell
1085
class chi_gen(rv_continuous):
1087
return sqrt(chi2.rvs(df,size=self._size))
1088
def _pdf(self, x, df):
1089
return x**(df-1.)*exp(-x*x*0.5)/(2.0)**(df*0.5-1)/gam(df*0.5)
1090
def _cdf(self, x, df):
1091
return special.gammainc(df*0.5,0.5*x*x)
1092
def _ppf(self, q, df):
1093
return sqrt(2*special.gammaincinv(df*0.5,q))
1094
def _stats(self, df):
1095
mu = sqrt(2)*special.gamma(df/2.0+0.5)/special.gamma(df/2.0)
1097
g1 = (2*mu**3.0 + mu*(1-2*df))/arr(mu2**1.5)
1098
g2 = 2*df*(1.0-df)-6*mu**4 + 4*mu**2 * (2*df-1)
1100
return mu, mu2, g1, g2
1101
chi = chi_gen(a=0.0,name='chi',shapes='df',extradoc="""
1108
## Chi-squared (gamma-distributed with loc=0 and scale=2 and shape=df/2)
1109
class chi2_gen(rv_continuous):
1111
return rand.chi2(df,self._size)
1112
def _pdf(self, x, df):
1113
Px = x**(df/2.0-1)*exp(-x/2.0)
1114
Px /= special.gamma(df/2.0)* 2**(df/2.0)
1116
def _cdf(self, x, df):
1117
return special.chdtr(df, x)
1118
def _sf(self, x, df):
1119
return special.chdtrc(df, x)
1120
def _isf(self, p, df):
1121
return special.chdtri(df, p)
1122
def _ppf(self, p, df):
1123
return self._isf(1.0-p, df)
1124
def _stats(self, df):
1129
return mu, mu2, g1, g2
1130
chi2 = chi2_gen(a=0.0,name='chi2',longname='A chi-squared',shapes='df',
1133
Chi-squared distribution
1137
## Cosine (Approximation to the Normal)
1138
class cosine_gen(rv_continuous):
1140
return 1.0/2/pi*(1+cos(x))
1142
return 1.0/2/pi*(pi + x + sin(x))
1144
return 0.0, pi*pi/3.0-2.0, 0.0, -6.0*(pi**4-90)/(5.0*(pi*pi-6)**2)
1146
return log(4*pi)-1.0
1147
cosine = cosine_gen(a=-pi,b=pi,name='cosine',extradoc="""
1149
Cosine distribution (Approximation to the normal)
1152
## Double Gamma distribution
1153
class dgamma_gen(rv_continuous):
1155
u = random(size=self._size)
1156
return (gamma.rvs(a,size=self._size)*Num.where(u>=0.5,1,-1))
1157
def _pdf(self, x, a):
1159
return 1.0/(2*special.gamma(a))*ax**(a-1.0) * exp(-ax)
1160
def _cdf(self, x, a):
1161
fac = 0.5*special.gammainc(a,abs(x))
1162
return where(x>0,0.5+fac,0.5-fac)
1163
def _sf(self, x, a):
1164
fac = 0.5*special.gammainc(a,abs(x))
1165
return where(x>0,0.5-0.5*fac,0.5+0.5*fac)
1166
def _ppf(self, q, a):
1167
fac = special.gammainccinv(a,1-abs(2*q-1))
1168
return where(q>0.5, fac, -fac)
1169
def _stats(self, a):
1171
return 0.0, mu2, 0.0, (a+2.0)*(a+3.0)/mu2-3.0
1172
dgamma = dgamma_gen(name='dgamma',longname="A double gamma",
1173
shapes='a',extradoc="""
1175
Double gamma distribution
1179
## Double Weibull distribution
1181
class dweibull_gen(rv_continuous):
1183
u = random(size=self._size)
1184
return weibull_min.rvs(c, size=self._size)*(Num.where(u>=0.5,1,-1))
1185
def _pdf(self, x, c):
1187
Px = c/2.0*ax**(c-1.0)*exp(-ax**c)
1189
def _cdf(self, x, c):
1190
Cx1 = 0.5*exp(-abs(x)**c)
1191
return where(x > 0, 1-Cx1, Cx1)
1192
def _ppf(self, q, c):
1193
fac = where(q<=0.5,2*q,2*q-1)
1194
fac = pow(arr(log(1.0/fac)),1.0/c)
1195
return where(q>0.5,fac,-fac)
1196
def _stats(self, c):
1198
return 0.0, var, 0.0, gam(1+4.0/c)/var
1199
dweibull = dweibull_gen(name='dweibull',longname="A double Weibull",
1200
shapes='c',extradoc="""
1202
Double Weibull distribution
1208
## Special case of the Gamma distribution with shape parameter an integer.
1210
class erlang_gen(rv_continuous):
1212
return gamma.rvs(n,size=self._size)
1213
def _arg_check(self, n):
1214
return (n > 0) & (floor(n)==n)
1215
def _pdf(self, x, n):
1216
Px = (x)**(n-1.0)*exp(-x)/special.gamma(n)
1218
def _cdf(self, x, n):
1219
return special.gdtr(1.0,n,x)
1220
def _sf(self, x, n):
1221
return special.gdtrc(1.0,n,x)
1222
def _ppf(self, q, n):
1223
return special.gdtrix(1.0, n, q)
1224
def _stats(self, n):
1226
return n, n, 2/sqrt(n), 6/n
1227
def _entropy(self, n):
1228
return special.psi(n)*(1-n) + 1 + special.gammaln(n)
1229
erlang = erlang_gen(a=0.0,name='erlang',longname='An Erlang',
1230
shapes='n',extradoc="""
1232
Erlang distribution (Gamma with integer shape parameter)
1236
## Exponential (gamma distributed with a=1.0, loc=loc and scale=scale)
1237
## scale == 1.0 / lambda
1239
class expon_gen(rv_continuous):
1241
return rand.standard_exp(self._size)
1249
return 1.0, 1.0, 2.0, 6.0
1252
expon = expon_gen(a=0.0,name='expon',longname="An exponential",
1255
Exponential distribution
1260
## Exponentiated Weibull
1261
class exponweib_gen(rv_continuous):
1262
def _pdf(self, x, a, c):
1264
return a*c*(1-exc)**arr(a-1) * exc * x**arr(c-1)
1265
def _cdf(self, x, a, c):
1267
return arr((1-exc)**a)
1268
def _ppf(self, q, a, c):
1269
return (-log(1-q**(1.0/a)))**arr(1.0/c)
1270
exponweib = exponweib_gen(a=0.0,name='exponweib',
1271
longname="An exponentiated Weibull",
1272
shapes="a,c",extradoc="""
1274
Exponentiated Weibull distribution
1278
## Exponential Power
1280
class exponpow_gen(rv_continuous):
1281
def _pdf(self, x, b):
1282
xbm1 = arr(x**(b-1.0))
1284
return exp(1)*b*xbm1 * exp(xb - exp(xb))
1285
def _cdf(self, x, b):
1287
return 1.0-exp(1-exp(xb))
1288
def _ppf(self, q, b):
1289
return pow(log(1.0-log(1.0-q)), 1.0/b)
1290
exponpow = exponpow_gen(a=0.0,name='exponpow',longname="An exponential power",
1291
shapes='b',extradoc="""
1293
Exponential Power distribution
1297
## Faigue-Life (Birnbaum-Sanders)
1298
class fatiguelife_gen(rv_continuous):
1300
z = norm.rvs(size=self._size)
1301
U = random(size=self._size)
1303
det = sqrt(fac*fac - 4)
1306
return t1*(U>0.5) + t2*(U<0.5)
1307
def _pdf(self, x, c):
1308
return (x+1)/arr(2*c*sqrt(2*pi*x**3))*exp(-(x-1)**2/arr((2.0*x*c**2)))
1309
def _cdf(self, x, c):
1310
return special.ndtr(1.0/c*(sqrt(x)-1.0/arr(sqrt(x))))
1311
def _ppf(self, q, c):
1312
tmp = c*special.ndtri(q)
1313
return 0.25*(tmp + sqrt(tmp**2 + 4))**2
1314
def _stats(self, c):
1319
g1 = 4*c*sqrt(11*c2+6.0)/den**1.5
1320
g2 = 6*c2*(93*c2+41.0) / den**2.0
1321
return mu, mu2, g1, g2
1322
fatiguelife = fatiguelife_gen(a=0.0,name='fatiguelife',
1323
longname="A fatigue-life (Birnbaum-Sanders)",
1324
shapes='c',extradoc="""
1326
Fatigue-life (Birnbaum-Sanders) distribution
1332
class foldcauchy_gen(rv_continuous):
1334
return abs(cauchy.rvs(loc=c,size=self._size))
1335
def _pdf(self, x, c):
1336
return 1.0/pi*(1.0/(1+(x-c)**2) + 1.0/(1+(x+c)**2))
1337
def _cdf(self, x, c):
1338
return 1.0/pi*(arctan(x-c) + arctan(x+c))
1339
def _stats(self, x, c):
1340
return inf, inf, nan, nan
1341
foldcauchy = foldcauchy_gen(a=0.0, name='foldcauchy',
1342
longname = "A folded Cauchy",
1343
shapes='c',extradoc="""
1345
A folded Cauchy distributions
1351
class f_gen(rv_continuous):
1352
def _rvs(self, dfn, dfd):
1353
return rand.f(dfn, dfd, self._size)
1354
def _pdf(self, x, dfn, dfd):
1357
Px = m**(m/2) * n**(n/2) * x**(n/2-1)
1358
Px /= (m+n*x)**((n+m)/2)*special.beta(n/2,m/2)
1360
def _cdf(self, x, dfn, dfd):
1361
return special.fdtr(dfn, dfd, x)
1362
def _sf(self, x, dfn, dfd):
1363
return special.fdtrc(dfn, dfd, x)
1364
def _ppf(self, q, dfn, dfd):
1365
return special.fdtri(dfn, dfd, q)
1366
def _stats(self, dfn, dfd):
1369
mu = where (v2 > 2, v2 / arr(v2 - 2), inf)
1370
mu2 = 2*v2*v2*(v2+v1-2)/(v1*(v2-2)**2 * (v2-4))
1371
mu2 = where(v2 > 4, mu2, inf)
1372
g1 = 2*(v2+2*v1-2)/(v2-6)*sqrt((2*v2-4)/(v1*(v2+v1-2)))
1373
g1 = where(v2 > 6, g1, nan)
1374
g2 = 3/(2*v2-16)*(8+g1*g1*(v2-6))
1375
g2 = where(v2 > 8, g2, nan)
1376
return mu, mu2, g1, g2
1377
f = f_gen(a=0.0,name='f',longname='An F',shapes="dfn,dfd",
1385
## abs(Z) where (Z is normal with mu=L and std=S so that c=abs(L)/S)
1387
## note: regress docs have scale parameter correct, but first parameter
1388
## he gives is a shape parameter A = c * scale
1390
## Half-normal is folded normal with shape-parameter c=0.
1392
class foldnorm_gen(rv_continuous):
1394
return abs(norm.rvs(loc=c,size=self._size))
1395
def _pdf(self, x, c):
1396
return sqrt(2.0/pi)*cosh(c*x)*exp(-(x*x+c*c)/2.0)
1397
def _cdf(self, x, c,):
1398
return special.ndtr(x-c) + special.ndtr(x+c) - 1.0
1399
def _stats(self, c):
1400
fac = special.erf(c/sqrt(2))
1401
mu = sqrt(2.0/pi)*exp(-0.5*c*c)+c*fac
1402
mu2 = c*c + 1 - mu*mu
1404
g1 = sqrt(2/pi)*exp(-1.5*c2)*(4-pi*exp(c2)*(2*c2+1.0))
1405
g1 += 2*c*fac*(6*exp(-c2) + 3*sqrt(2*pi)*c*exp(-c2/2.0)*fac + \
1409
g2 = c2*c2+6*c2+3+6*(c2+1)*mu*mu - 3*mu**4
1410
g2 -= 4*exp(-c2/2.0)*mu*(sqrt(2.0/pi)*(c2+2)+c*(c2+3)*exp(c2/2.0)*fac)
1412
return mu, mu2, g1, g2
1413
foldnorm = foldnorm_gen(a=0.0,name='foldnorm',longname='A folded normal',
1414
shapes='c',extradoc="""
1416
Folded normal distribution
1420
## Extreme Value Type II or Frechet
1421
## (defined in Regress+ documentation as Extreme LB) as
1422
## a limiting value distribution.
1424
class frechet_r_gen(rv_continuous):
1425
def _pdf(self, x, c):
1426
return c*pow(x,c-1)*exp(-pow(x,c))
1427
def _cdf(self, x, c):
1428
return 1-exp(-pow(x,c))
1429
def _ppf(self, q, c):
1430
return pow(-log(1-q),1.0/c)
1431
def _munp(self, n, c):
1432
return special.gamma(1.0+n*1.0/c)
1433
def _entropy(self, c):
1434
return -_EULER / c - log(c) + _EULER + 1
1435
frechet_r = frechet_r_gen(a=0.0,name='frechet_r',longname="A Frechet right",
1436
shapes='c',extradoc="""
1438
A Frechet (right) distribution (also called Weibull minimum)
1441
weibull_min = frechet_r_gen(a=0.0,name='weibull_min',
1442
longname="A Weibull minimum",
1443
shapes='c',extradoc="""
1445
A Weibull minimum distribution
1449
class frechet_l_gen(rv_continuous):
1450
def _pdf(self, x, c):
1451
return c*pow(-x,c-1)*exp(-pow(-x,c))
1452
def _cdf(self, x, c):
1453
return exp(-pow(-x,c))
1454
def _ppf(self, q, c):
1455
return -pow(-log(q),1.0/c)
1456
def _munp(self, n, c):
1457
val = special.gamma(1.0+n*1.0/c)
1458
if (int(n) % 2): sgn = -1
1461
def _entropy(self, c):
1462
return -_EULER / c - log(c) + _EULER + 1
1463
frechet_l = frechet_l_gen(b=0.0,name='frechet_l',longname="A Frechet left",
1464
shapes='c',extradoc="""
1466
A Frechet (left) distribution (also called Weibull maximum)
1469
weibull_max = frechet_l_gen(b=0.0,name='weibull_max',
1470
longname="A Weibull maximum",
1471
shapes='c',extradoc="""
1473
A Weibull maximum distribution
1478
## Generalized Logistic
1480
class genlogistic_gen(rv_continuous):
1481
def _pdf(self, x, c):
1482
Px = c*exp(-x)/(1+exp(-x))**(c+1.0)
1484
def _cdf(self, x, c):
1485
Cx = (1+exp(-x))**(-c)
1487
def _ppf(self, q, c):
1488
vals = -log(pow(q,-1.0/c)-1)
1490
def _stats(self, c):
1492
mu = _EULER + special.psi(c)
1493
mu2 = pi*pi/6.0 + zeta(2,c)
1494
g1 = -2*zeta(3,c) + 2*_ZETA3
1496
g2 = pi**4/15.0 + 6*zeta(4,c)
1498
return mu, mu2, g1, g2
1499
genlogistic = genlogistic_gen(name='genlogistic',
1500
longname="A generalized logistic",
1501
shapes='c',extradoc="""
1503
Generalized logistic distribution
1507
## Generalized Pareto
1508
class genpareto_gen(rv_continuous):
1509
def _argcheck(self, c):
1511
self.b = where(c < 0, 1.0/abs(c), inf)
1512
return where(c==0, 0, 1)
1513
def _pdf(self, x, c):
1514
Px = pow(1+c*x,arr(-1.0-1.0/c))
1516
def _cdf(self, x, c):
1517
return 1.0 - pow(1+c*x,arr(-1.0/c))
1518
def _ppf(self, q, c):
1519
vals = 1.0/c * (pow(1-q, -c)-1)
1521
def _munp(self, n, c):
1523
val = (-1.0/c)**n * sum(scipy.comb(n,k)*(-1)**k / (1.0-c*k))
1524
return where(c*n < 1, val, inf)
1525
def _entropy(self, c):
1530
return rv_continuous._entropy(self, c)
1531
genpareto = genpareto_gen(a=0.0,name='genpareto',
1532
longname="A generalized Pareto",
1533
shapes='c',extradoc="""
1535
Generalized Pareto distribution
1539
## Generalized Exponential
1541
class genexpon_gen(rv_continuous):
1542
def _pdf(self, x, a, b, c):
1543
return (a+b*(1-exp(-c*x)))*exp((a-b)*x+b*(1-exp(-c*x))/c)
1544
def _cdf(self, x, a, b, c):
1545
return 1.0-exp((a-b)*x + b*(1-exp(-c*x))/c)
1546
genexpon = genexpon_gen(a=0.0,name='genexpon',
1547
longname='A generalized exponential',
1548
shapes='a,b,c',extradoc="""
1550
Generalized exponential distribution
1554
## Generalized Extreme Value
1555
## c=0 is just gumbel distribution.
1556
## This version does not accept c==0
1557
## Use gumbel_r for c==0
1559
class genextreme_gen(rv_continuous):
1560
def _argcheck(self, c):
1561
self.b = where(c > 0, 1.0 / c, inf)
1562
self.a = where(c < 0, 1.0 / c, -inf)
1564
def _pdf(self, x, c):
1566
pex2 = pow(ex2,1.0/c)
1567
p2 = exp(-pex2)*pex2/ex2
1569
def _cdf(self, x, c):
1570
return exp(-pow(1-c*x,1.0/c))
1571
def _ppf(self, q, c):
1572
return 1.0/c*(1-(-log(q))**c)
1573
def _munp(self, n, c):
1575
vals = 1.0/c**n * sum(scipy.comb(n,k) * (-1)**k * special.gamma(c*k + 1))
1576
return where(c*n > -1, vals, inf)
1577
genextreme = genextreme_gen(name='genextreme',
1578
longname="A generalized extreme value",
1579
shapes='c',extradoc="""
1581
Generalized extreme value (see gumbel_r for c=0)
1585
## Gamma (Use MATLAB and MATHEMATICA (b=theta=scale, a=alpha=shape) definition)
1587
## gamma(a, loc, scale) with a an integer is the Erlang distribution
1588
## gamma(1, loc, scale) is the Exponential distribution
1589
## gamma(df/2, 0, 2) is the chi2 distribution with df degrees of freedom.
1591
class gamma_gen(rv_continuous):
1593
return rand.standard_gamma(a, self._size)
1594
def _pdf(self, x, a):
1595
return x**(a-1)*exp(-x)/special.gamma(a)
1596
def _cdf(self, x, a):
1597
return special.gammainc(a, x)
1598
def _ppf(self, q, a):
1599
return special.gammaincinv(a,q)
1600
def _stats(self, a):
1601
return a, a, 2.0/sqrt(a), 6.0/a
1602
def _entropy(self, a):
1603
return special.psi(a)*(1-a) + 1 + special.gammaln(a)
1604
gamma = gamma_gen(a=0.0,name='gamma',longname='A gamma',
1605
shapes='a',extradoc="""
1612
class gengamma_gen(rv_continuous):
1613
def _argcheck(self, a, c):
1614
return (a > 0) & (c != 0)
1615
def _pdf(self, x, a, c):
1616
return abs(c)* x**(c*a-1) / special.gamma(a) * exp(-x**c)
1617
def _cdf(self, x, a, c):
1618
val = special.gammainc(a,x**c)
1620
return where(cond>0,val,1-val)
1621
def _ppf(self, q, a, c):
1622
val1 = special.gammaincinv(a,q)
1623
val2 = special.gammaincinv(a,1.0-q)
1626
return where(cond > 0,val1**ic,val2**ic)
1627
def _munp(self, n, a, c):
1628
return special.gamma(a+n*1.0/c) / special.gamma(a)
1630
val = special.psi(a)
1631
return a*(1-val) + 1.0/c*val + special.gammaln(a)-log(abs(c))
1632
gengamma = gengamma_gen(a=0.0, name='gengamma',
1633
longname='A generalized gamma',
1634
shapes="a,c", extradoc="""
1636
Generalized gamma distribution
1640
## Generalized Half-Logistic
1643
class genhalflogistic_gen(rv_continuous):
1644
def _argcheck(self, c):
1647
def _pdf(self, x, c):
1650
tmp0 = tmp**(limit-1)
1652
return 2*tmp0 / (1+tmp2)**2
1653
def _cdf(self, x, c):
1657
return (1.0-tmp2) / (1+tmp2)
1658
def _ppf(self, q, c):
1659
return 1.0/c*(1-((1.0-q)/(1.0+q))**c)
1660
def _entropy(self,c):
1661
return 2 - (2*c+1)*log(2)
1662
genhalflogistic = genhalflogistic_gen(a=0.0, name='genhalflogistic',
1663
longname="A generalized half-logistic",
1664
shapes='c',extradoc="""
1666
Generalized half-logistic
1670
## Gompertz (Truncated Gumbel)
1673
class gompertz_gen(rv_continuous):
1674
def _pdf(self, x, c):
1676
return c*ex*exp(-c*(ex-1))
1677
def _cdf(self, x, c):
1678
return 1.0-exp(-c*(exp(x)-1))
1679
def _ppf(self, q, c):
1680
return log(1-1.0/c*log(1-q))
1681
def _entropy(self, c):
1682
return 1.0 - log(c) - exp(c)*special.expn(1,c)
1683
gompertz = gompertz_gen(a=0.0, name='gompertz',
1684
longname="A Gompertz (truncated Gumbel) distribution",
1685
shapes='c',extradoc="""
1687
Gompertz (truncated Gumbel) distribution
1691
## Gumbel, Log-Weibull, Fisher-Tippett, Gompertz
1692
## The left-skewed gumbel distribution.
1693
## and right-skewed are available as gumbel_l and gumbel_r
1695
class gumbel_r_gen(rv_continuous):
1700
return exp(-exp(-x))
1702
return -log(-log(q))
1704
return _EULER, pi*pi/6.0, \
1705
12*sqrt(6)/pi**3 * _ZETA3, 12.0/5
1707
return 1.0608407169541684911
1708
gumbel_r = gumbel_r_gen(name='gumbel_r',longname="A (right-skewed) Gumbel",
1711
Right-skewed Gumbel (Log-Weibull, Fisher-Tippett, Gompertz) distribution
1714
class gumbel_l_gen(rv_continuous):
1719
return 1.0-exp(-exp(x))
1721
return log(-log(1-q))
1723
return _EULER, pi*pi/6.0, \
1724
12*sqrt(6)/pi**3 * _ZETA3, 12.0/5
1726
return 1.0608407169541684911
1727
gumbel_l = gumbel_l_gen(name='gumbel_l',longname="A left-skewed Gumbel",
1730
Left-skewed Gumbel distribution
1736
class halfcauchy_gen(rv_continuous):
1738
return 2.0/pi/(1.0+x*x)
1740
return 2.0/pi*arctan(x)
1744
return inf, inf, nan, nan
1747
halfcauchy = halfcauchy_gen(a=0.0,name='halfcauchy',
1748
longname="A Half-Cauchy",extradoc="""
1750
Half-Cauchy distribution
1758
class halflogistic_gen(rv_continuous):
1760
return 0.5/(cosh(x/2.0))**2.0
1766
if n==1: return 2*log(2)
1767
if n==2: return pi*pi/3.0
1768
if n==3: return 9*_ZETA3
1769
if n==4: return 7*pi**4 / 15.0
1770
return 2*(1-pow(2.0,1-n))*special.gamma(n+1)*special.zeta(n,1)
1773
halflogistic = halflogistic_gen(a=0.0, name='halflogistic',
1774
longname="A half-logistic",
1777
Half-logistic distribution
1782
## Half-normal = chi(1, loc, scale)
1784
class halfnorm_gen(rv_continuous):
1786
return abs(norm.rvs(size=self._size))
1788
return sqrt(2.0/pi)*exp(-x*x/2.0)
1790
return special.ndtr(x)*2-1.0
1792
return special.ndtri((1+q)/2.0)
1794
return sqrt(2.0/pi), 1-2.0/pi, sqrt(2)*(4-pi)/(pi-2)**1.5, \
1797
return 0.5*log(pi/2.0)+0.5
1798
halfnorm = halfnorm_gen(a=0.0, name='halfnorm',
1799
longname="A half-normal",
1802
Half-normal distribution
1806
## Hyperbolic Secant
1808
class hypsecant_gen(rv_continuous):
1810
return 1.0/(pi*cosh(x))
1812
return 2.0/pi*arctan(exp(x))
1814
return log(tan(pi*q/2.0))
1816
return 0, pi*pi/4, 0, 2
1819
hypsecant = hypsecant_gen(name='hypsecant',longname="A hyperbolic secant",
1822
Hyperbolic secant distribution
1826
## Gauss Hypergeometric
1828
class gausshyper_gen(rv_continuous):
1829
def _argcheck(self, a, b, c, z):
1830
return (a > 0) & (b > 0) & (c==c) & (z==z)
1831
def _pdf(self, x, a, b, c, z):
1832
Cinv = gam(a)*gam(b)/gam(a+b)*special.hyp2f1(c,a,a+b,-z)
1833
return 1.0/Cinv * x**(a-1.0) * (1.0-x)**(b-1.0) / (1.0+z*x)**c
1834
def _munp(self, n, a, b, c, z):
1835
fac = special.beta(n+a,b) / special.beta(a,b)
1836
num = special.hyp2f1(c,a+n,a+b+n,-z)
1837
den = special.hyp2f1(c,a,a+b,-z)
1838
return fac*num / den
1839
gausshyper = gausshyper_gen(a=0.0, b=1.0, name='gausshyper',
1840
longname="A Gauss hypergeometric",
1844
Gauss hypergeometric distribution
1849
# special case of generalized gamma with c=-1
1852
class invgamma_gen(rv_continuous):
1853
def _pdf(self, x, a):
1854
return x**(-a-1) / special.gamma(a) * exp(-1.0/x)
1855
def _cdf(self, x, a):
1856
return 1.0-special.gammainc(a, 1.0/x)
1857
def _ppf(self, q, a):
1858
return 1.0/special.gammaincinv(a,1-q)
1859
def _munp(self, n, a):
1860
return special.gamma(a-n) / special.gamma(a)
1861
def _entropy(self, a):
1862
return a - (a+1.0)*special.psi(a) + special.gammaln(a)
1863
invgamma = invgamma_gen(a=0.0, name='invgamma',longname="An inverted gamma",
1864
shapes='a',extradoc="""
1866
Inverted gamma distribution
1871
## Inverse Normal Distribution
1872
# scale is gamma from DATAPLOT and B from Regress
1874
class invnorm_gen(rv_continuous):
1876
return rv._inst._Wald(mu,size=(self._size,))
1877
def _pdf(self, x, mu):
1878
return 1.0/sqrt(2*pi*x**3.0)*exp(-1.0/(2*x)*((x-mu)/mu)**2)
1879
def _cdf(self, x, mu):
1881
C1 = norm.cdf(fac*(x-mu)/mu)
1882
C1 += exp(2.0/mu)*norm.cdf(-fac*(x+mu)/mu)
1884
def _stats(self, mu):
1885
return mu, mu**3.0, 3*sqrt(mu), 15*mu
1886
invnorm = invnorm_gen(a=0.0, name='invnorm', longname="An inverse normal",
1887
shapes="mu",extradoc="""
1889
Inverse normal distribution
1895
class invweibull_gen(rv_continuous):
1896
def _pdf(self, x, c):
1900
def _cdf(self, x, c):
1903
def _ppf(self, q, c):
1904
return pow(-log(q),arr(-1.0/c))
1905
def _entropy(self, c):
1906
return 1+_EULER + _EULER / c - log(c)
1907
invweibull = invweibull_gen(a=0,name='invweibull',
1908
longname="An inverted Weibull",
1909
shapes='c',extradoc="""
1911
Inverted Weibull distribution
1917
class johnsonsb_gen(rv_continuous):
1918
def _argcheck(self, a, b):
1919
return (b > 0) & (a==a)
1920
def _pdf(self, x, a, b):
1921
trm = norm.pdf(a+b*log(x/(1.0-x)))
1922
return b*1.0/(x*(1-x))*trm
1923
def _cdf(self, x, a, b):
1924
return norm.cdf(a+b*log(x/(1.0-x)))
1925
def _ppf(self, q, a, b):
1926
return 1.0/(1+exp(-1.0/b*norm.ppf(q)-a))
1927
johnsonsb = johnsonsb_gen(a=0.0,b=1.0,name='johnsonb',
1928
longname="A Johnson SB",
1929
shapes="a,b",extradoc="""
1931
Johnson SB distribution
1936
class johnsonsu_gen(rv_continuous):
1937
def _argcheck(self, a, b):
1938
return (b > 0) & (a==a)
1939
def _pdf(self, x, a, b):
1941
trm = norm.pdf(a+b*log(x+sqrt(x2+1)))
1942
return b*1.0/sqrt(x2+1.0)*trm
1943
def _cdf(self, x, a, b):
1944
return norm.cdf(a+b*log(x+sqrt(x*x+1)))
1945
def _ppf(self, q, a, b):
1946
return sinh((norm.ppf(q)-a)/b)
1947
johnsonsu = johnsonsu_gen(name='johnsonsu',longname="A Johnson SU",
1948
shapes="a,b", extradoc="""
1950
Johnson SB distribution
1955
## Laplace Distribution
1957
class laplace_gen(rv_continuous):
1959
return 0.5*exp(-abs(x))
1961
return where(x > 0, 1.0-0.5*exp(-x), 0.5*exp(x))
1963
return where(q > 0.5, -log(2*(1-q)), log(2*q))
1968
laplace = laplace_gen(name='laplace', longname="A Laplace",
1971
Laplacian distribution
1976
## Levy Distribution
1978
class levy_gen(rv_continuous):
1980
return 1/sqrt(2*x)/x*exp(-1/(2*x))
1982
return 2*(1-norm._cdf(1/sqrt(x)))
1984
val = norm._ppf(1-q/2.0)
1985
return 1.0/(val*val)
1987
return inf, inf, nan, nan
1988
levy = levy_gen(a=0.0,name="levy", longname = "A Levy", extradoc="""
1994
## Left-skewed Levy Distribution
1996
class levy_l_gen(rv_continuous):
1999
return 1/sqrt(2*ax)/ax*exp(-1/(2*ax))
2002
return 2*norm._cdf(1/sqrt(ax))-1
2004
val = norm._ppf((q+1.0)/2)
2005
return -1.0/(val*val)
2007
return inf, inf, nan, nan
2008
levy_l = levy_l_gen(b=0.0,name="levy_l", longname = "A left-skewed Levy", extradoc="""
2010
Left-skewed Levy distribution
2014
## Levy-stable Distribution (only random variates)
2016
class levy_stable_gen(rv_continuous):
2017
def _rvs(self, alpha, beta):
2019
TH = uniform.rvs(loc=-pi/2.0,scale=pi,size=sz)
2020
W = expon.rvs(size=sz)
2022
return 2/pi*(pi/2+beta*TH)*tan(TH)-beta*log((pi/2*W*cos(TH))/(pi/2+beta*TH))
2027
return W/(cos(TH)/tan(aTH)+sin(TH))*((cos(aTH)+sin(aTH)*tan(TH))/W)**ialpha
2029
val0 = beta*tan(pi*alpha/2)
2030
th0 = arctan(val0)/alpha
2031
val3 = W/(cos(TH)/tan(alpha*(th0+TH))+sin(TH))
2032
res3 = val3*((cos(aTH)+sin(aTH)*tan(TH)-val0*(sin(aTH)-cos(aTH)*tan(TH)))/W)**ialpha
2035
def _argcheck(self, alpha, beta):
2040
return (alpha > 0) & (alpha <= 2) & (beta <= 1) & (beta >= -1)
2042
def _pdf(self, x, alpha, beta):
2043
raise NotImplementedError
2045
levy_stable = levy_stable_gen(name='levy_stable', longname="A Levy-stable",
2046
shapes="alpha, beta", extradoc="""
2048
Levy-stable distribution (only random variates available -- ignore other docs)
2053
## Logistic (special case of generalized logistic with c=1)
2056
class logistic_gen(rv_continuous):
2059
return ex / (1+ex)**2.0
2061
return 1.0/(1+exp(-x))
2063
return -log(1.0/q-1)
2065
return 0, pi*pi/3.0, 0, 6.0/5.0
2068
logistic = logistic_gen(name='logistic', longname="A logistic",
2071
Logistic distribution
2078
class loggamma_gen(rv_continuous):
2079
def _pdf(self, x, c):
2080
return exp(c*x-exp(x))/ special.gamma(c)
2081
def _cdf(self, x, c):
2082
return special.gammainc(c, exp(x))/ special.gamma(c)
2083
def _ppf(self, q, c):
2084
return log(special.gammaincinv(c,q*special.gamma(c)))
2085
loggamma = loggamma_gen(name='loggamma', longname="A log gamma",
2088
Log gamma distribution
2092
## Log-Laplace (Log Double Exponential)
2095
class loglaplace_gen(rv_continuous):
2096
def _pdf(self, x, c):
2098
c = where(x < 1, c, -c)
2100
def _cdf(self, x, c):
2101
return where(x < 1, 0.5*x**c, 1-0.5*x**(-c))
2102
def _ppf(self, q, c):
2103
return where(q < 0.5, (2.0*q)**(1.0/c), (2*(1.0-q))**(-1.0/c))
2104
def _entropy(self, c):
2105
return log(2.0/c) + 1.0
2106
loglaplace = loglaplace_gen(a=0.0, name='loglaplace',
2107
longname="A log-Laplace",shapes='c',
2110
Log-Laplace distribution (Log Double Exponential)
2114
## Lognormal (Cobb-Douglass)
2115
## std is a shape parameter and is the variance of the underlying
2117
## the mean of the underlying distribution is log(scale)
2119
class lognorm_gen(rv_continuous):
2121
return exp(s * norm.rvs(size=self._size))
2122
def _pdf(self, x, s):
2123
Px = exp(-log(x)**2 / (2*s**2))
2124
return Px / (s*x*sqrt(2*pi))
2125
def _cdf(self, x, s):
2126
return norm.cdf(log(x)/s)
2127
def _ppf(self, q, s):
2128
return exp(s*norm._ppf(q))
2129
def _stats(self, s):
2133
g1 = sqrt((p-1))*(2+p)
2134
g2 = scipy.polyval([1,2,3,0,-6.0],p)
2135
return mu, mu2, g1, g2
2136
def _entropy(self, s):
2137
return 0.5*(1+log(2*pi)+2*log(s))
2138
lognorm = lognorm_gen(a=0.0, name='lognorm',
2139
longname='A lognormal', shapes='s',
2142
Lognormal distribution
2146
# Gibrat's distribution is just lognormal with s=1
2148
class gilbrat_gen(lognorm_gen):
2150
return lognorm_gen._rvs(self, 1.0)
2152
return lognorm_gen._pdf(self, x, 1.0)
2154
return lognorm_gen._cdf(self, x, 1.0)
2156
return lognorm_gen._ppf(self, q, 1.0)
2158
return lognorm_gen._stats(self, 1.0)
2160
return 0.5*log(2*pi) + 0.5
2161
gilbrat = gilbrat_gen(a=0.0, name='gilbrat', longname='A Gilbrat',
2164
Gilbrat distribution
2170
# a special case of chi with df = 3, loc=0.0, and given scale = 1.0/sqrt(a)
2171
# where a is the parameter used in mathworld description
2173
class maxwell_gen(rv_continuous):
2175
return chi.rvs(3.0,size=self._size)
2177
return sqrt(2.0/pi)*x*x*exp(-x*x/2.0)
2179
return special.gammainc(1.5,x*x/2.0)
2181
return sqrt(2*special.gammaincinv(1.5,q))
2184
return 2*sqrt(2.0/pi), 3-8/pi, sqrt(2)*(32-10*pi)/val**1.5, \
2185
(-12*pi*pi + 160*pi - 384) / val**2.0
2187
return _EULER + 0.5*log(2*pi)-0.5
2188
maxwell = maxwell_gen(a=0.0, name='maxwell', longname="A Maxwell",
2191
Maxwell distribution
2195
# Mielke's Beta-Kappa
2197
class mielke_gen(rv_continuous):
2198
def _pdf(self, x, k, s):
2199
return k*x**(k-1.0) / (1.0+x**s)**(1.0+k*1.0/s)
2200
def _cdf(self, x, k, s):
2201
return x**k / (1.0+x**s)**(k*1.0/s)
2202
def _ppf(self, q, k, s):
2203
qsk = pow(q,s*1.0/k)
2204
return pow(qsk/(1.0-qsk),1.0/s)
2205
mielke = mielke_gen(a=0.0, name='mielke', longname="A Mielke's Beta-Kappa",
2206
shapes="k,s", extradoc="""
2208
Mielke's Beta-Kappa distribution
2214
class nakagami_gen(rv_continuous):
2215
def _pdf(self, x, nu):
2216
return 2*nu**nu/gam(nu)*(x**(2*nu-1.0))*exp(-nu*x*x)
2217
def _cdf(self, x, nu):
2218
return special.gammainc(nu,nu*x*x)
2219
def _ppf(self, q, nu):
2220
return sqrt(1.0/nu*special.gammaincinv(nu,q))
2221
def _stats(self, nu):
2222
mu = gam(nu+0.5)/gam(nu)/sqrt(nu)
2224
g1 = mu*(1-4*nu*mu2)/2.0/nu/mu2**1.5
2225
g2 = -6*mu**4*nu + (8*nu-2)*mu**2-2*nu + 1
2227
return mu, mu2, g1, g2
2228
nakagami = nakagami_gen(a=0.0, name="nakagami", longname="A Nakagami",
2229
shapes='nu', extradoc="""
2231
Nakagami distribution
2236
# Non-central chi-squared
2237
# nc is lambda of definition, df is nu
2239
class ncx2_gen(rv_continuous):
2240
def _rvs(self, df, nc):
2241
return rand.noncentral_chi2(df,nc,self._size)
2242
def _pdf(self, x, df, nc):
2244
Px = exp(-nc/2.0)*special.hyp0f1(a,nc*x/4.0)
2245
Px *= exp(-x/2.0)*x**(a-1) / arr(2**a * special.gamma(a))
2247
def _cdf(self, x, df, nc):
2248
return special.chndtr(x,df,nc)
2249
def _ppf(self, q, df, nc):
2250
return special.chndtrix(q,df,nc)
2251
def _stats(self, df, nc):
2253
return df + nc, 2*val, sqrt(8)*(val+nc)/val**1.5, \
2254
12.0*(val+2*nc)/val**2.0
2255
ncx2 = ncx2_gen(a=0.0, name='ncx2', longname="A non-central chi-squared",
2256
shapes="df,nc", extradoc="""
2258
Non-central chi-squared distribution
2264
class ncf_gen(rv_continuous):
2265
def _rvs(self, dfn, dfd, nc):
2266
return rand.noncentral_f(dfn,dfd,nc,self._size)
2267
def _pdf(self, x, dfn, dfd, nc):
2269
Px = exp(-nc/2+nc*n1*x/(2*(n2+n1*x)))
2270
Px *= n1**(n1/2) * n2**(n2/2) * x**(n1/2-1)
2271
Px *= (n2+n1*x)**(-(n1+n2)/2)
2272
Px *= special.gamma(n1/2)*special.gamma(1+n2/2)
2273
Px *= special.assoc_laguerre(-nc*n1*x/(2.0*(n2+n1*x)),n2/2,n1/2-1)
2274
Px /= special.beta(n1/2,n2/2)*special.gamma((n1+n2)/2.0)
2275
def _cdf(self, x, dfn, dfd, nc):
2276
return special.ncfdtr(dfn,dfd,nc,x)
2277
def _ppf(self, q, dfn, dfd, nc):
2278
return special.ncfdtri(dfn, dfd, nc, q)
2279
def _munp(self, n, dfn, dfd, nc):
2280
val = (dfn *1.0/dfd)**n
2281
val *= gam(n+0.5*dfn)*gam(0.5*dfd-n) / gam(dfd*0.5)
2282
val *= exp(-nc / 2.0)
2283
val *= special.hyp1f1(n+0.5*dfn, 0.5*dfn, 0.5*nc)
2285
def _stats(self, dfn, dfd, nc):
2286
mu = where(dfd <= 2, inf, dfd / (dfd-2.0)*(1+nc*1.0/dfn))
2287
mu2 = where(dfd <=4, inf, 2*(dfd*1.0/dfn)**2.0 * \
2288
((dfn+nc/2.0)**2.0 + (dfn+nc)*(dfd-2.0)) / \
2289
((dfd-2.0)**2.0 * (dfd-4.0)))
2290
return mu, mu2, None, None
2291
ncf = ncf_gen(a=0.0, name='ncf', longname="A non-central F distribution",
2292
shapes="dfn,dfd,nc", extradoc="""
2294
Non-central F distribution
2298
## Student t distribution
2300
class t_gen(rv_continuous):
2302
Y = f.rvs(df, df, size=self._size)
2304
return 0.5*sqrt(df)*(sY-1.0/sY)
2305
def _pdf(self, x, df):
2307
Px = exp(special.gammaln((r+1)/2)-special.gammaln(r/2))
2308
Px /= sqrt(r*pi)*(1+(x**2)/r)**((r+1)/2)
2310
def _cdf(self, x, df):
2311
return special.stdtr(df, x)
2312
def _ppf(self, q, df):
2313
return special.stdtrit(df, q)
2314
def _stats(self, df):
2315
mu2 = where(df > 2, df / (df-2.0), inf)
2316
g1 = where(df > 3, 0.0, nan)
2317
g2 = where(df > 4, 6.0/(df-4.0), nan)
2318
return 0, mu2, g1, g2
2319
t = t_gen(name='t',longname="Student's T",
2320
shapes="df", extradoc="""
2322
Student's T distribution
2326
## Non-central T distribution
2328
class nct_gen(rv_continuous):
2329
def _rvs(self, df, nc):
2330
return norm.rvs(loc=nc,size=self._size)*sqrt(df) / sqrt(chi2.rvs(df,size=self._size))
2331
def _pdf(self, x, df, nc):
2337
Px = n**(n/2) * special.gamma(n+1)
2338
Px /= arr(2.0**n*exp(nc*nc/2)*fac1**(n/2)*special.gamma(n/2))
2339
valF = ncx2 / (2*fac1)
2340
trm1 = sqrt(2)*nc*x*special.hyp1f1(n/2+1,1.5,valF)
2341
trm1 /= arr(fac1*special.gamma((n+1)/2))
2342
trm2 = special.hyp1f1((n+1)/2,0.5,valF)
2343
trm2 /= arr(sqrt(fac1)*special.gamma(n/2+1))
2346
def _cdf(self, x, df, nc):
2347
return special.nctdtr(df, nc, x)
2348
def _ppf(self, q, df, nc):
2349
return special.nctdtrit(df, nc, q)
2350
def _stats(self, df, nc, moments='mv'):
2351
mu, mu2, g1, g2 = None, None, None, None
2352
val1 = gam((df-1.0)/2.0)
2355
mu = nc*sqrt(df/2.0)*val1/val2
2357
var = (nc*nc+1.0)*df/(df-2.0)
2358
var -= nc*nc*df* val1**2 / 2.0 / val2**2
2361
g1n = 2*nc*sqrt(df)*val1*((nc*nc*(2*df-7)-3)*val2**2 \
2362
-nc*nc*(df-2)*(df-3)*val1**2)
2363
g1d = (df-3)*sqrt(2*df*(nc*nc+1)/(df-2) - \
2364
nc*nc*df*(val1/val2)**2) * val2 * \
2365
(nc*nc*(df-2)*val1**2 - \
2366
2*(nc*nc+1)*val2**2)
2369
g2n = 2*(-3*nc**4*(df-2)**2 *(df-3) *(df-4)*val1**4 + \
2370
2**(6-2*df) * nc*nc*(df-2)*(df-4)* \
2371
(nc*nc*(2*df-7)-3)*pi* gam(df+1)**2 - \
2372
4*(nc**4*(df-5)-6*nc*nc-3)*(df-3)*val2**4)
2373
g2d = (df-3)*(df-4)*(nc*nc*(df-2)*val1**2 - \
2374
2*(nc*nc+1)*val2)**2
2376
return mu, mu2, g1, g2
2377
nct = nct_gen(name="nct", longname="A Noncentral T",
2378
shapes="df,nc", extradoc="""
2380
Non-central Student T distribution
2386
class pareto_gen(rv_continuous):
2387
def _pdf(self, x, b):
2388
return b * x**(-b-1)
2389
def _cdf(self, x, b):
2391
def _ppf(self, q, b):
2392
return pow(1-q, -1.0/b)
2393
def _stats(self, b, moments='mv'):
2394
mu, mu2, g1, g2 = None, None, None, None
2397
bt = extract(mask,b)
2398
mu = valarray(shape(b),value=inf)
2399
insert(mu, mask, bt / (bt-1.0))
2402
bt = extract( mask,b)
2403
mu2 = valarray(shape(b), value=inf)
2404
insert(mu2, mask, bt / (bt-2.0) / (bt-1.0)**2)
2407
bt = extract( mask,b)
2408
g1 = valarray(shape(b), value=nan)
2409
vals = 2*(bt+1.0)*sqrt(b-2.0)/((b-3.0)*sqrt(b))
2410
insert(g1, mask, vals)
2413
bt = extract( mask,b)
2414
g2 = valarray(shape(b), value=nan)
2415
vals = 6.0*polyval([1.0,1.0,-6,-2],bt)/ \
2416
polyval([1.0,-7.0,12.0,0.0],bt)
2417
insert(g2, mask, vals)
2418
return mu, mu2, g1, g2
2419
def _entropy(self, c):
2420
return 1 + 1.0/c - log(c)
2421
pareto = pareto_gen(a=1.0, name="pareto", longname="A Pareto",
2422
shapes="b", extradoc="""
2428
# LOMAX (Pareto of the second kind.)
2429
# Special case of Pareto of the first kind (location=-1.0)
2431
class lomax_gen(rv_continuous):
2432
def _pdf(self, x, c):
2433
return c*1.0/(1.0+x)**(c+1.0)
2434
def _cdf(self, x, c):
2435
return 1.0-1.0/(1.0+x)**c
2436
def _ppf(self, q, c):
2437
return pow(1.0-q,-1.0/c)-1
2438
def _stats(self, c):
2439
mu, mu2, g1, g2 = pareto.stats(c, loc=-1.0, moments='mvsk')
2440
return mu, mu2, g1, g2
2441
def _entropy(self, c):
2442
return 1+1.0/c-log(c)
2443
lomax = lomax_gen(a=0.0, name="lomax",
2444
longname="A Lomax (Pareto of the second kind)",
2445
shapes="c", extradoc="""
2447
Lomax (Pareto of the second kind) distribution
2450
## Power-function distribution
2451
## Special case of beta dist. with d =1.0
2453
class powerlaw_gen(rv_continuous):
2454
def _pdf(self, x, a):
2456
def _cdf(self, x, a):
2458
def _ppf(self, q, a):
2459
return pow(q, 1.0/a)
2460
def _stats(self, a):
2461
return a/(a+1.0), a*(a+2.0)/(a+1.0)**2, \
2462
2*(1.0-a)*sqrt((a+2.0)/(a*(a+3.0))), \
2463
6*polyval([1,-1,-6,2],a)/(a*(a+3.0)*(a+4))
2464
def _entropy(self, a):
2465
return 1 - 1.0/a - log(a)
2466
powerlaw = powerlaw_gen(a=0.0, b=1.0, name="powerlaw",
2467
longname="A power-function",
2468
shapes="a", extradoc="""
2470
Power-function distribution
2476
class powerlognorm_gen(rv_continuous):
2477
def _pdf(self, x, c, s):
2478
return c/(x*s)*norm.pdf(log(x)/s)
2479
def _cdf(self, x, c, s):
2480
return 1.0 - pow(norm.cdf(-log(x)/s),c*1.0)
2481
def _ppf(self, q, c, s):
2482
return exp(-s*norm.ppf(pow(1.0-q,1.0/c)))
2483
powerlognorm = powerlognorm_gen(a=0.0, name="powerlognorm",
2484
longname="A power log-normal",
2485
shapes="c,s", extradoc="""
2487
Power log-normal distribution
2493
class powernorm_gen(norm_gen):
2494
def _pdf(self, x, c):
2495
return c*norm_gen._pdf(self, x)* \
2496
(norm_gen._cdf(self, -x)**(c-1.0))
2497
def _cdf(self, x, c):
2498
return 1.0-norm_gen._cdf(self, -x)**(c*1.0)
2499
def _ppf(self, q, c):
2500
return -norm_gen._ppf(self, pow(1.0-q,1.0/c))
2501
powernorm = powernorm_gen(name='powernorm', longname="A power normal",
2502
shapes="c", extradoc="""
2504
Power normal distribution
2508
# R-distribution ( a general-purpose distribution with a
2509
# variety of shapes.
2511
class rdist_gen(rv_continuous):
2512
def _pdf(self, x, c):
2513
return pow((1.0-x*x),c/2.0-1) / special.beta(0.5,c/2.0)
2514
def _cdf(self, x, c):
2515
return 0.5 + x/special.beta(0.5,c/2.0)* \
2516
special.hyp2f1(0.5,1.0-c/2.0,1.5,x*x)
2518
return (1-(n % 2))*special.beta((n+1.0)/2,c/2.0)
2519
rdist = rdist_gen(a=-1.0,b=1.0, name="rdist", longname="An R-distributed",
2520
shapes="c", extradoc="""
2526
# Rayleigh distribution (this is chi with df=2 and loc=0.0)
2527
# scale is the mode.
2529
class rayleigh_gen(rv_continuous):
2531
return chi.rvs(2,size=self._size)
2533
return r*exp(-r*r/2.0)
2535
return 1.0-exp(-r*r/2.0)
2537
return sqrt(-2*log(1-q))
2540
return pi/2, val/2, 2*(pi-3)*sqrt(pi)/val**1.5, \
2543
return _EULER/2.0 + 1 - 0.5*log(2)
2544
rayleigh = rayleigh_gen(a=0.0, name="rayleigh",
2545
longname="A Rayleigh",
2548
Rayleigh distribution
2552
# Reciprocal Distribution
2553
class reciprocal_gen(rv_continuous):
2554
def _argcheck(self, a, b):
2557
self.d = log(b*1.0 / a)
2558
return (a > 0) & (b > 0) & (b > a)
2559
def _pdf(self, x, a, b):
2560
# argcheck should be called before _pdf
2561
return 1.0/(x*self.d)
2562
def _cdf(self, x, a, b):
2563
return (log(x)-log(a)) / self.d
2564
def _ppf(self, q, a, b):
2565
return a*pow(b*1.0/a,q)
2566
def _munp(self, n, a, b):
2567
return 1.0/self.d / n * (pow(b*1.0,n) - pow(a*1.0,n))
2568
def _entropy(self,a,b):
2569
return 0.5*log(a*b)+log(log(b/a))
2570
reciprocal = reciprocal_gen(name="reciprocal",
2571
longname="A reciprocal",
2572
shapes="a,b", extradoc="""
2574
Reciprocal distribution
2580
class rice_gen(rv_continuous):
2581
def _pdf(self, x, b):
2582
return x*exp(-(x*x+b*b)/2.0)*special.i0(x*b)
2583
def _munp(self, n, b):
2587
return 2.0**(nd2)*exp(-b2)*special.gamma(n1) * \
2588
special.hyp1f1(n1,1,b2)
2589
rice = rice_gen(a=0.0, name="rice", longname="A Rice",
2590
shapes="b", extradoc="""
2596
# Reciprocal Inverse Gaussian
2598
class recipinvgauss_gen(rv_continuous):
2599
def _pdf(self, x, mu):
2600
return 1.0/sqrt(2*pi*x)*exp(-(1-mu*x)**2.0 / (2*x*mu**2.0))
2601
def _cdf(self, x, mu):
2605
return 1.0-norm.cdf(isqx*trm1)-exp(2.0/mu)*norm.cdf(-isqx*trm2)
2606
recipinvgauss = recipinvgauss_gen(a=0.0, name='recipinvgauss',
2607
longname="A reciprocal inverse Gaussian",
2608
shapes="mu", extradoc="""
2610
Reciprocal inverse Gaussian
2616
class semicircular_gen(rv_continuous):
2618
return 2.0/pi*sqrt(1-x*x)
2620
return 0.5+1.0/pi*(x*sqrt(1-x*x) + arcsin(x))
2622
return 0, 0.25, 0, -1.0
2624
return 0.64472988584940017414
2625
semicircular = semicircular_gen(a=-1.0,b=1.0, name="semicircular",
2626
longname="A semicircular",
2629
Semicircular distribution
2634
# up-sloping line from loc to (loc + c*scale) and then downsloping line from
2635
# loc + c*scale to loc + scale
2637
# _trstr = "Left must be <= mode which must be <= right with left < right"
2638
class triang_gen(rv_continuous):
2639
def _argcheck(self, c):
2640
return (c >= 0) & (c <= 1)
2641
def _pdf(self, x, c):
2642
return where(x < c, 2*x/c, 2*(1-x)/(1-c))
2643
def _cdf(self, x, c):
2644
return where(x < c, x*x/c, (x*x-2*x+c)/(c-1))
2645
def _ppf(self, q, c):
2646
return where(q < c, sqrt(c*q), 1-sqrt((1-c)*(1-q)))
2647
def _stats(self, c):
2648
return (c+1.0)/3.0, (1.0-c+c*c)/18, sqrt(2)*(2*c-1)*(c+1)*(c-2) / \
2649
(5*(1.0-c+c*c)**1.5), -3.0/5.0
2650
def _entropy(self,c):
2652
triang = triang_gen(a=0.0, b=1.0, name="triang", longname="A Triangular",
2653
shapes="c", extradoc="""
2655
Triangular distribution
2657
up-sloping line from loc to (loc + c*scale) and then downsloping
2658
for (loc + c*scale) to (loc+scale).
2660
- standard form is in the range [0,1] with c the mode.
2661
- location parameter shifts the start to loc
2662
- scale changes the width from 1 to scale
2666
# Truncated Exponential
2668
class truncexpon_gen(rv_continuous):
2669
def _argcheck(self, b):
2672
def _pdf(self, x, b):
2673
return exp(-x)/(1-exp(-b))
2674
def _cdf(self, x, b):
2675
return (1.0-exp(-x))/(1-exp(-b))
2676
def _ppf(self, q, b):
2677
return -log(1-q+q*exp(-b))
2678
def _munp(self, n, b):
2679
return gam(n+1)-special.gammainc(1+n,b)
2680
def _entropy(self, b):
2682
return log(eB-1)+(1+eB*(b-1.0))/(1.0-eB)
2683
truncexpon = truncexpon_gen(a=0.0, name='truncexpon',
2684
longname="A truncated exponential",
2685
shapes="b", extradoc="""
2687
Truncated exponential distribution
2693
class truncnorm_gen(rv_continuous):
2694
def _argcheck(self, a, b):
2697
self.nb = norm._cdf(b)
2698
self.na = norm._cdf(a)
2700
def _pdf(self, x, a, b):
2701
return norm._pdf(x) / (self.nb - self.na)
2702
def _cdf(self, x, a, b):
2703
return (norm._cdf(x) - self.na) / (self.nb - self.na)
2704
def _ppf(self, q, a, b):
2705
return norm._ppf(q*self.nb + self.na*(1.0-q))
2706
def _stats(self, a, b):
2707
nA, nB = self.na, self.nb
2709
pA, pB = norm._pdf(a), norm._pdf(b)
2711
mu2 = 1 + (a*pA - b*pB) / d - mu*mu
2712
return mu, mu2, None, None
2713
truncnorm = truncnorm_gen(name='truncnorm', longname="A truncated normal",
2714
shapes="a,b", extradoc="""
2716
Truncated Normal distribution.
2718
The standard form of this distribution is a standard normal truncated to the
2719
range [a,b] --- notice that a and b are defined over the domain
2720
of the standard normal. To convert clip values for a specific mean and
2721
standard deviation use a,b = (myclip_a-my_mean)/my_std, (myclip_b-my_mean)/my_std
2726
# A flexible distribution ranging from Cauchy (lam=-1)
2727
# to logistic (lam=0.0)
2728
# to approx Normal (lam=0.14)
2729
# to u-shape (lam = 0.5)
2730
# to Uniform from -1 to 1 (lam = 1)
2732
class tukeylambda_gen(rv_continuous):
2733
def _pdf(self, x, lam):
2734
Fx = arr(special.tklmbda(x,lam))
2735
Px = Fx**(lam-1.0) + (arr(1-Fx))**(lam-1.0)
2737
return where((lam > 0) & (abs(x) < 1.0/lam), Px, 0.0)
2738
def _cdf(self, x, lam):
2739
return special.tklmbda(x, lam)
2740
def _ppf(self, q, lam):
2742
vals1 = (q**lam - (1-q)**lam)/lam
2743
vals2 = log(q/(1-q))
2744
return where((lam == 0)&(q==q), vals2, vals1)
2745
def _stats(self, lam):
2746
mu2 = 2*gam(lam+1.5)-lam*pow(4,-lam)*sqrt(pi)*gam(lam)*(1-2*lam)
2747
mu2 /= lam*lam*(1+2*lam)*gam(1+1.5)
2748
mu4 = 3*gam(lam)*gam(lam+0.5)*pow(2,-2*lam) / lam**3 / gam(2*lam+1.5)
2749
mu4 += 2.0/lam**4 / (1+4*lam)
2750
mu4 -= 2*sqrt(3)*gam(lam)*pow(2,-6*lam)*pow(3,3*lam) * \
2751
gam(lam+1.0/3)*gam(lam+2.0/3) / (lam**3.0 * gam(2*lam+1.5) * \
2753
g2 = mu4 / mu2 / mu2 - 3.0
2755
return 0, mu2, 0, g2
2756
def _entropy(self, lam):
2758
return log(pow(p,lam-1)+pow(1-p,lam-1))
2759
return scipy.integrate.quad(integ,0,1)[0]
2760
tukeylambda = tukeylambda_gen(name='tukeylambda', longname="A Tukey-Lambda",
2761
shapes="lam", extradoc="""
2763
Tukey-Lambda distribution
2765
A flexible distribution ranging from Cauchy (lam=-1)
2766
to logistic (lam=0.0)
2767
to approx Normal (lam=0.14)
2768
to u-shape (lam = 0.5)
2769
to Uniform from -1 to 1 (lam = 1)
2774
# loc to loc + scale
2776
class uniform_gen(rv_continuous):
2778
return rand.uniform(0.0,1.0,self._size)
2786
return 0.5, 1.0/12, 0, -1.2
2789
uniform = uniform_gen(a=0.0,b=1.0, name='uniform', longname="A uniform",
2792
Uniform distribution
2794
constant between loc and loc+scale
2800
# if x is not in range or loc is not in range it assumes they are angles
2801
# and converts them to [-pi, pi] equivalents.
2803
eps = scipy_base.limits.double_epsilon
2805
class vonmises_gen(rv_continuous):
2807
return rv._inst._von_Mises(b,size=(self._size,))
2808
def _pdf(self, x, b):
2809
x = arr(angle(exp(1j*x)))
2810
Px = where(b < 100, exp(b*cos(x)) / (2*pi*special.i0(b)),
2811
norm.pdf(x, 0.0, sqrt(1.0/b)))
2813
def _cdf(self, x, b):
2814
x = arr(angle(exp(1j*x)))
2817
c_xsimple = atleast_1d((b==0)&(x==x))
2818
c_xiter = atleast_1d((b<100)&(b > 0)&(x==x))
2819
c_xnormal = atleast_1d((b>=100)&(x==x))
2820
c_bad = atleast_1d((b<=0) | (x != x))
2822
indxiter = nonzero(c_xiter)
2823
xiter = take(x, indxiter)
2825
vals = ones(len(c_xsimple),Num.Float)
2826
putmask(vals, c_bad, nan)
2827
putmask(vals, c_xsimple, x / 2.0/pi)
2829
st = where(isnan(st),0.0,st)
2830
putmask(vals, c_xnormal, norm.cdf(x, scale=st))
2832
biter = take(atleast_1d(b)*(x==x), indxiter)
2834
fac = special.i0(biter)
2837
for j in range(1,501):
2838
trm1 = special.iv(j,biter)/j/fac
2841
if all(trm1 < eps2):
2844
print "Warning: did not converge..."
2845
put(vals, indxiter, val)
2847
def _stats(self, b):
2848
return 0, None, 0, None
2849
vonmises = vonmises_gen(name='vonmises', longname="A Von Mises",
2850
shapes="b", extradoc="""
2852
Von Mises distribution
2854
if x is not in range or loc is not in range it assumes they are angles
2855
and converts them to [-pi, pi] equivalents.
2861
## Wald distribution (Inverse Normal with shape parameter mu=1.0)
2863
class wald_gen(invnorm_gen):
2865
return invnorm_gen._rvs(self, 1.0)
2867
return invnorm.pdf(x,1.0)
2869
return invnorm.cdf(x,1,0)
2871
return 1.0, 1.0, 3.0, 15.0
2872
wald = wald_gen(a=0.0, name="wald", longname="A Wald",
2884
class wrapcauchy_gen(rv_continuous):
2885
def _argcheck(self, c):
2886
return (c > 0) & (c < 1)
2887
def _pdf(self, x, c):
2888
return (1.0-c*c)/(2*pi*(1+c*c-2*c*cos(x)))
2889
def _cdf(self, x, c):
2891
val = (1.0+c)/(1.0-c)
2895
valp = extract(c1,val)
2897
valn = extract(c2,val)
2901
on = 1.0-1.0/pi*arctan(valn*yn)
2902
insert(output, c2, on)
2905
op = 1.0/pi*arctan(valp*yp)
2906
insert(output, c1, op)
2908
def _ppf(self, q, c):
2909
val = (1.0-c)/(1.0+c)
2910
rcq = 2*arctan(val*tan(pi*q))
2911
rcmq = 2*pi-2*arctan(val*tan(pi*(1-q)))
2912
return where(q < 1.0/2, rcq, rcmq)
2913
def _entropy(self, c):
2914
return log(2*pi*(1-c*c))
2915
wrapcauchy = wrapcauchy_gen(a=0.0,b=2*pi, name='wrapcauchy',
2916
longname="A wrapped Cauchy",
2917
shapes="c", extradoc="""
2919
Wrapped Cauchy distribution
2923
### DISCRETE DISTRIBUTIONS
2926
def entropy(pk,qk=None):
2927
"""S = entropy(pk,qk=None)
2929
calculate the entropy of a distribution given the p_k values
2930
S = -sum(pk * log(pk))
2932
If qk is not None, then compute a relative entropy
2933
S = -sum(pk * log(pk / qk))
2935
Routine will normalize pk and qk if they don't sum to 1
2938
pk = 1.0* pk / sum(pk)
2940
vec = where(pk == 0, 0.0, pk*log(pk))
2943
if len(qk) != len(pk):
2944
raise ValueError, "qk and pk must have same length."
2945
qk = 1.0*qk / sum(qk)
2946
# If qk is zero anywhere, then unless pk is zero at those places
2947
# too, the relative entropy is infinite.
2948
if any(take(pk,nonzero(qk==0.0))!=0.0):
2950
vec = where (pk == 0, 0.0, pk*log(pk / qk))
2954
## Handlers for generic case where xk and pk are given
2956
def _drv_pdf(self, xk, *args):
2962
def _drv_cdf(self, xk, *args):
2963
indx = argmax((self.xk>xk))-1
2964
return self.F[self.xk[indx]]
2966
def _drv_ppf(self, q, *args):
2967
indx = argmax((self.qvals>=q))
2968
return self.Finv[self.qvals[indx]]
2970
def _drv_nonzero(self, k, *args):
2973
def _drv_moment(self, n, *args):
2975
return sum(self.xk**n[NewAxis,...] * self.pk, axis=0)
2977
def _drv_moment_gen(self, t, *args):
2979
return sum(exp(self.xk * t[NewAxis,...]) * self.pk, axis=0)
2981
def _drv2_moment(self, n, *args):
2986
while (pos <= self.b) and ((pos >= (self.b + self.a)/2.0) and \
2987
(diff > self.moment_tol)):
2988
diff = pos**n * self._pdf(pos,*args)
2992
def _drv2_ppfsingle(self, q, *args): # Use basic bisection algorithm
2995
if isinf(b): # Be sure ending point is > q
2998
if b >= self.b: qb = 1.0; break
2999
qb = self._cdf(b,*args)
3000
if (qb < q): b += 10
3004
if isinf(a): # be sure starting point < q
3007
if a <= self.a: qb = 0.0; break
3008
qa = self._cdf(a,*args)
3009
if (qa > q): a -= 10
3012
qa = self._cdf(a, *args)
3022
qc = self._cdf(c, *args)
3032
def reverse_dict(dict):
3034
for key in dict.keys():
3035
newdict[dict[key]] = key
3038
def make_dict(keys, values):
3040
for key, value in zip(keys, values):
3044
# Must over-ride one of _pmf or _cdf or pass in
3045
# x_k, p(x_k) lists in initialization
3048
"""A generic discrete random variable.
3050
Discrete random variables are defined from a standard form.
3051
The standard form may require some other parameters to complete
3052
its specification. The distribution methods also take an optional location
3053
parameter using loc= keyword. The default is loc=0. The calling form
3054
of the methods follow:
3056
generic.rvs(<shape(s)>,loc=0)
3059
generic.pmf(x,<shape(s)>,loc=0)
3060
- probability mass function
3062
generic.cdf(x,<shape(s)>,loc=0)
3063
- cumulative density function
3065
generic.sf(x,<shape(s)>,loc=0)
3066
- survival function (1-cdf --- sometimes more accurate)
3068
generic.ppf(q,<shape(s)>,loc=0)
3069
- percent point function (inverse of cdf --- percentiles)
3071
generic.isf(q,<shape(s)>,loc=0)
3072
- inverse survival function (inverse of sf)
3074
generic.stats(<shape(s)>,loc=0,moments='mv')
3075
- mean('m'), variance('v'), skew('s'), and/or kurtosis('k')
3077
generic.entropy(<shape(s)>,loc=0)
3080
Alternatively, the object may be called (as a function) to fix
3081
the shape and location parameters returning a
3082
"frozen" discrete RV object:
3084
myrv = generic(<shape(s)>,loc=0)
3085
- frozen RV object with the same methods but holding the
3086
given shape and location fixed.
3088
You can construct an aribtrary discrete rv where P{X=xk} = pk
3089
by passing to the rv_discrete initialization method (through the values=
3090
keyword) a tuple of sequences (xk,pk) which describes only those values of
3091
X (xk) that occur with nonzero probability (pk).
3093
def __init__(self, a=0, b=inf, name=None, badvalue=None,
3094
moment_tol=1e-8,values=None,inc=1,longname=None,
3095
shapes=None, extradoc=None):
3096
if badvalue is None:
3098
self.badvalue = badvalue
3104
self.moment_tol = moment_tol
3106
self._cdfvec = sgf(self._cdfsingle,otypes='d')
3107
self.return_integers = 1
3108
self.vecentropy = vectorize(self._entropy)
3110
if values is not None:
3111
self.xk, self.pk = values
3112
self.return_integers = 0
3113
indx = argsort(ravel(self.xk))
3114
self.xk = take(ravel(self.xk),indx)
3115
self.pk = take(ravel(self.pk),indx)
3117
self.b = self.xk[-1]
3118
self.P = make_dict(self.xk, self.pk)
3119
self.qvals = scipy_base.cumsum(self.pk)
3120
self.F = make_dict(self.xk, self.qvals)
3121
self.Finv = reverse_dict(self.F)
3122
self._ppf = new.instancemethod(sgf(_drv_ppf,otypes='d'),
3124
self._pmf = new.instancemethod(sgf(_drv_pmf,otypes='d'),
3126
self._cdf = new.instancemethod(sgf(_drv_cdf,otypes='d'),
3128
self._nonzero = new.instancemethod(_drv_nonzero, self, rv_discrete)
3129
self.generic_moment = new.instancemethod(_drv_moment,
3131
self.moment_gen = new.instancemethod(_drv_moment_gen,
3135
self._vecppf = new.instancemethod(sgf(_drv2_ppfsingle,otypes='d'),
3137
self.generic_moment = new.instancemethod(sgf(_drv2_moment,
3140
cdf_signature = inspect.getargspec(self._cdf.im_func)
3141
numargs1 = len(cdf_signature[0]) - 2
3142
pmf_signature = inspect.getargspec(self._pmf.im_func)
3143
numargs2 = len(pmf_signature[0]) - 2
3144
self.numargs = max(numargs1, numargs2)
3146
if longname is None:
3147
if name[0] in ['aeiouAEIOU']: hstr = "An "
3149
longname = hstr + name
3150
if self.__doc__ is None:
3151
self.__doc__ = rv_discrete.__doc__
3152
if self.__doc__ is not None:
3153
self.__doc__ = self.__doc__.replace("A Generic",longname)
3154
if name is not None:
3155
self.__doc__ = self.__doc__.replace("generic",name)
3157
self.__doc__ = self.__doc__.replace("<shape(s)>,","")
3159
self.__doc__ = self.__doc__.replace("<shape(s)>",shapes)
3160
ind = self.__doc__.find("You can construct an arbitrary")
3161
self.__doc__ = self.__doc__[:ind].strip()
3162
if extradoc is not None:
3163
self.__doc__ = self.__doc__ + extradoc
3165
def _rvs(self, *args):
3166
return self._ppf(rand.sample(self._size),*args)
3168
def __fix_loc(self, args, loc):
3170
if N > self.numargs:
3171
if N == self.numargs + 1 and loc is None: # loc is given without keyword
3173
args = args[:self.numargs]
3178
def _nonzero(self, k, *args):
3181
def _argcheck(self, *args):
3187
def _pmf(self, k, *args):
3188
return self.cdf(k,*args) - self.cdf(k-1,*args)
3190
def _cdfsingle(self, k, *args):
3191
m = arange(int(self.a),k+1)
3192
return sum(self._pmf(m,*args))
3194
def _cdf(self, x, *args):
3196
return self._cdfvec(k,*args)
3198
def _sf(self, x, *args):
3199
return 1.0-self._cdf(x,*args)
3201
def _ppf(self, q, *args):
3202
return self._vecppf(q, *args)
3204
def _isf(self, q, *args):
3205
return self._ppf(1-q,*args)
3207
def _stats(self, *args):
3208
return None, None, None, None
3210
def _munp(self, n, *args):
3211
return self.generic_moment(n)
3214
def rvs(self, *args, **kwds):
3215
loc,size=map(kwds.get,['loc','size'])
3216
args, loc = self.__fix_loc(args, loc)
3217
cond = self._argcheck(*args)
3219
raise ValueError, "Domain error in arguments."
3224
self._size = product(size)
3225
if scipy.isscalar(size):
3229
vals = reshape(self._rvs(*args),size)
3230
if self.return_integers:
3232
if vals.typecode() not in scipy.typecodes['AllInteger']:
3233
vals = vals.astype(Num.Int)
3236
def pmf(self, k,*args, **kwds):
3237
"""Probability mass function at k of the given RV.
3241
The shape parameter(s) for the distribution (see docstring of the
3242
instance object for more information)
3246
loc - location parameter (default=0)
3248
loc = kwds.get('loc')
3249
args, loc = self.__fix_loc(args, loc)
3250
k,loc = map(arr,(k,loc))
3251
args = tuple(map(arr,args))
3253
cond0 = self._argcheck(*args)
3254
cond1 = (k >= self.a) & (k <= self.b) & self._nonzero(k,*args)
3255
cond = cond0 & cond1
3256
output = zeros(shape(cond),'d')
3257
insert(output,(1-cond0)*(cond1==cond1),self.badvalue)
3258
goodargs = argsreduce(cond, *((k,)+args))
3259
insert(output,cond,self._pmf(*goodargs))
3262
def cdf(self, k, *args, **kwds):
3263
"""Cumulative distribution function at k of the given RV
3267
The shape parameter(s) for the distribution (see docstring of the
3268
instance object for more information)
3272
loc - location parameter (default=0)
3274
loc = kwds.get('loc')
3275
args, loc = self.__fix_loc(args, loc)
3276
k,loc = map(arr,(k,loc))
3277
args = tuple(map(arr,args))
3279
cond0 = self._argcheck(*args)
3280
cond1 = (k >= self.a) & (k < self.b)
3281
cond2 = (k >= self.b)
3282
cond = cond0 & cond1
3283
output = zeros(shape(cond),'d')
3284
insert(output,(1-cond0)*(cond1==cond1),self.badvalue)
3285
insert(output,cond2*(cond0==cond0), 1.0)
3286
goodargs = argsreduce(cond, *((k,)+args))
3287
insert(output,cond,self._cdf(*goodargs))
3290
def sf(self,k,*args,**kwds):
3291
"""Survival function (1-cdf) at k of the given RV
3295
The shape parameter(s) for the distribution (see docstring of the
3296
instance object for more information)
3300
loc - location parameter (default=0)
3302
loc= kwds.get('loc')
3303
args, loc = self.__fix_loc(args, loc)
3304
k,loc = map(arr,(k,loc))
3305
args = tuple(map(arr,args))
3307
cond0 = self._argcheck(*args)
3308
cond1 = (k >= self.a) & (k <= self.b)
3309
cond2 = (k < self.a) & cond0
3310
cond = cond0 & cond1
3311
output = zeros(shape(cond),'d')
3312
insert(output,(1-cond0)*(cond1==cond1),self.badvalue)
3313
insert(output,cond2,1.0)
3314
goodargs = argsreduce(cond, *((k,)+args))
3315
insert(output,cond,self._sf(*goodargs))
3318
def ppf(self,q,*args,**kwds):
3319
"""Percent point function (inverse of cdf) at q of the given RV
3323
The shape parameter(s) for the distribution (see docstring of the
3324
instance object for more information)
3328
loc - location parameter (default=0)
3330
loc = kwds.get('loc')
3331
args, loc = self.__fix_loc(args, loc)
3332
q,loc = map(arr,(q,loc))
3333
args = tuple(map(arr,args))
3334
cond0 = self._argcheck(*args) & (loc == loc)
3335
cond1 = (q > 0) & (q < 1)
3336
cond2 = (q==1) & cond0
3337
cond = cond0 & cond1
3338
output = valarray(shape(cond),value=self.a-1)
3339
insert(output,(1-cond0)*(cond1==cond1), self.badvalue)
3340
insert(output,cond2,self.b)
3341
goodargs = argsreduce(cond, *((q,)+args+(loc,)))
3342
loc, goodargs = goodargs[-1], goodargs[:-1]
3343
insert(output,cond,self._ppf(*goodargs) + loc)
3346
def isf(self,q,*args,**kwds):
3347
"""Inverse survival function (1-sf) at q of the given RV
3351
The shape parameter(s) for the distribution (see docstring of the
3352
instance object for more information)
3356
loc - location parameter (default=0)
3359
loc = kwds.get('loc')
3360
args, loc = self.__fix_loc(args, loc)
3361
q,loc = map(arr,(q,loc))
3362
args = tuple(map(arr,args))
3363
cond0 = self._argcheck(*args) & (loc == loc)
3364
cond1 = (q > 0) & (q < 1)
3365
cond2 = (q==1) & cond0
3366
cond = cond0 & cond1
3367
output = valarray(shape(cond),value=self.b)
3368
insert(output,(1-cond0)*(cond1==cond1), self.badvalue)
3369
insert(output,cond2,self.a-1)
3370
goodargs = argsreduce(cond, *((q,)+args+(loc,)))
3371
loc, goodargs = goodargs[-1], goodargs[:-1]
3372
insert(output,cond,self._ppf(*goodargs) + loc)
3375
def stats(self, *args, **kwds):
3376
"""Some statistics of the given discrete RV
3380
The shape parameter(s) for the distribution (see docstring of the
3381
instance object for more information)
3385
loc - location parameter (default=0)
3386
moments - a string composed of letters ['mvsk'] specifying
3387
which moments to compute (default='mv')
3390
's' = (Fisher's) skew,
3391
'k' = (Fisher's) kurtosis.
3393
loc,moments=map(kwds.get,['loc','moments'])
3395
if N > self.numargs:
3396
if N == self.numargs + 1 and loc is None: # loc is given without keyword
3398
if N == self.numargs + 2 and moments is None: # loc, scale, and moments
3399
loc, moments = args[-2:]
3400
args = args[:self.numargs]
3401
if loc is None: loc = 0.0
3402
if moments is None: moments = 'mv'
3405
args = tuple(map(arr,args))
3406
cond = self._argcheck(*args) & (loc==loc)
3408
signature = inspect.getargspec(self._stats.im_func)
3409
if (signature[2] is not None) or ('moments' in signature[0]):
3410
mu, mu2, g1, g2 = self._stats(*args,**{'moments':moments})
3412
mu, mu2, g1, g2 = self._stats(*args)
3417
default = valarray(shape(cond), self.badvalue)
3420
# Use only entries that are valid in calculation
3421
goodargs = argsreduce(cond, *(args+(loc,)))
3422
loc, goodargs = goodargs[-1], goodargs[:-1]
3426
mu = self._munp(1.0,*goodargs)
3427
out0 = default.copy()
3428
insert(out0,cond,mu+loc)
3433
mu2p = self._munp(2.0,*goodargs)
3435
mu = self._munp(1.0,*goodargs)
3437
out0 = default.copy()
3438
insert(out0,cond,mu2)
3443
mu3p = self._munp(3.0,*goodargs)
3445
mu = self._munp(1.0,*goodargs)
3447
mu2p = self._munp(2.0,*goodargs)
3449
mu3 = mu3p - 3*mu*mu2 - mu**3
3451
out0 = default.copy()
3452
insert(out0,cond,g1)
3457
mu4p = self._munp(4.0,*goodargs)
3459
mu = self._munp(1.0,*goodargs)
3461
mu2p = self._munp(2.0,*goodargs)
3464
mu3p = self._munp(3.0,*goodargs)
3465
mu3 = mu3p - 3*mu*mu2 - mu**3
3466
mu4 = mu4p - 4*mu*mu3 - 6*mu*mu*mu2 - mu**4
3467
g2 = mu4 / mu2**2.0 - 3.0
3468
out0 = default.copy()
3469
insert(out0,cond,g2)
3472
if len(output) == 1:
3475
return tuple(output)
3477
def moment(self, n, *args, **kwds): # Non-central moments in standard form.
3479
raise ValueError, "Moment must be an integer."
3480
if (n < 0): raise ValueError, "Moment must be positive."
3481
if (n == 0): return 1.0
3482
if (n > 0) and (n < 5):
3483
signature = inspect.getargspec(self._stats.im_func)
3484
if (signature[2] is not None) or ('moments' in signature[0]):
3485
dict = {'moments':{1:'m',2:'v',3:'vs',4:'vk'}[n]}
3488
mu, mu2, g1, g2 = self._stats(*args,**dict)
3490
if mu is None: return self._munp(1,*args)
3493
if mu2 is None: return self._munp(2,*args)
3496
if g1 is None or mu2 is None: return self._munp(3,*args)
3497
else: return g1*(mu2**1.5)
3499
if g2 is None or mu2 is None: return self._munp(4,*args)
3500
else: return (g2+3.0)*(mu2**2.0)
3502
return self._munp(n,*args)
3504
def freeze(self, *args, **kwds):
3505
return rv_frozen(self, *args, **kwds)
3507
def _entropy(self, *args):
3508
if hasattr(self,'pk'):
3509
return entropy(self.pk)
3511
mu = int(self.stats(*args, **{'moments':'m'}))
3512
val = self.pmf(mu,*args)
3513
if (val==0.0): ent = 0.0
3514
else: ent = -val*log(val)
3517
while (abs(term) > eps):
3518
val = self.pmf(mu+k,*args)
3519
if val == 0.0: term = 0.0
3520
else: term = -val * log(val)
3521
val = self.pmf(mu-k,*args)
3522
if val != 0.0: term -= val*log(val)
3527
def entropy(self, *args, **kwds):
3528
loc= kwds.get('loc')
3529
args, loc = self.__fix_loc(args, loc)
3531
args = map(arr,args)
3532
cond0 = self._argcheck(*args) & (loc==loc)
3533
output = zeros(shape(cond0),'d')
3534
insert(output,(1-cond0),self.badvalue)
3535
goodargs = argsreduce(cond0, *args)
3536
insert(output,cond0,self.vecentropy(*goodargs))
3539
def __call__(self, *args, **kwds):
3540
return self.freeze(*args,**kwds)
3544
class binom_gen(rv_discrete):
3545
def _rvs(self, n, pr):
3546
return rand.binomial(n,pr,self._size)
3547
def _argcheck(self, n, pr):
3549
return (n>=0) & (pr >= 0) & (pr <= 1)
3550
def _cdf(self, x, n, pr):
3552
vals = special.bdtr(k,n,pr)
3554
def _sf(self, x, n, pr):
3556
return special.bdtrc(k,n,pr)
3557
def _ppf(self, q, n, pr):
3558
vals = ceil(special.bdtrik(q,n,pr))
3560
temp = special.bdtr(vals1,n,pr)
3561
return where(temp >= q, vals1, vals)
3562
def _stats(self, n, pr):
3566
g1 = (q-pr) / sqrt(n*pr*q)
3567
g2 = (1.0-6*pr*q)/(n*pr*q)
3568
return mu, var, g1, g2
3569
def _entropy(self, n, pr):
3571
vals = self._pmf(k,n,pr)
3572
lvals = where(vals==0,0.0,log(vals))
3573
return -sum(vals*lvals)
3574
binom = binom_gen(name='binom',shapes="n,pr",extradoc="""
3576
Binomial distribution
3578
Counts the number of successes in *n* independent
3579
trials when the probability of success each time is *pr*.
3582
# Bernoulli distribution
3584
class bernoulli_gen(binom_gen):
3586
return binom_gen._rvs(self, 1, pr)
3587
def _argcheck(self, pr):
3588
return (pr >=0 ) & (pr <= 1)
3589
def _cdf(self, x, pr):
3590
return binom_gen._cdf(self, x, 1, pr)
3591
def _sf(self, x, pr):
3592
return binom_gen._sf(self, x, 1, pr)
3593
def _ppf(self, q, pr):
3594
return binom_gen._ppf(self, q, 1, pr)
3595
def _stats(self, pr):
3596
return binom_gen._stats(self, 1, pr)
3597
def _entropy(self, pr):
3598
return -pr*log(pr)-(1-pr)*log(1-pr)
3599
bernoulli = bernoulli_gen(b=1,name='bernoulli',shapes="pr",extradoc="""
3601
Bernoulli distribution
3603
1 if binary experiment succeeds, 0 otherwise. Experiment
3604
succeeds with probabilty *pr*.
3609
class nbinom_gen(rv_discrete):
3610
def _rvs(self, n, pr):
3611
return rand.negative_binomial(n, pr, self._size)
3612
def _argcheck(self, n, pr):
3613
return (n >= 0) & (pr >= 0) & (pr <= 1)
3614
def _cdf(self, x, n, pr):
3616
return special.nbdtr(k,n,pr)
3617
def _sf(self, x, n, pr):
3619
return special.nbdtrc(k,n,pr)
3620
def _ppf(self, q, n, pr):
3621
vals = ceil(special.nbdtrik(q,n,pr))
3623
temp = special.nbdtr(vals1,n,pr)
3624
return where(temp >= q, vals1, vals)
3625
def _stats(self, n, pr):
3630
g1 = (Q+P)/sqrt(n*P*Q)
3631
g2 = (1.0 + 6*P*Q) / (n*P*Q)
3632
return mu, var, g1, g2
3633
nbinom = nbinom_gen(name='nbinom', longname="A negative binomial",
3634
shapes="n,pr", extradoc="""
3636
Negative binomial distribution
3640
## Geometric distribution
3642
class geom_gen(rv_discrete):
3644
return rv._inst._geom(pr,size=(self._size,))
3645
def _argcheck(self, pr):
3646
return (pr<=1) & (pr >= 0)
3647
def _pmf(self, k, pr):
3648
return (1-pr)**k * pr
3649
def _cdf(self, x, pr):
3651
return (1.0-(1.0-pr)**k)
3652
def _sf(self, x, pr):
3655
def _ppf(self, q, pr):
3656
vals = ceil(log(1.0-q)/log(1-pr))
3657
temp = 1.0-(1.0-pr)**(vals-1)
3658
return where((temp >= q) & (vals > 0), vals-1, vals)
3659
def _stats(self, pr):
3663
g1 = (2.0-pr) / sqrt(qr)
3664
g2 = scipy.polyval([1,-6,6],pr)/(1.0-pr)
3665
return mu, var, g1, g2
3666
geom = geom_gen(a=1,name='geom', longname="A geometric",
3667
shapes="pr", extradoc="""
3669
Geometric distribution
3673
## Hypergeometric distribution
3675
class hypergeom_gen(rv_discrete):
3676
def _rvs(self, M, n, N):
3677
return rv._inst._hypergeom(M,n,N,size=(self._size,))
3678
def _argcheck(self, M, n, N):
3679
cond = rv_discrete._argcheck(self,M,n,N)
3680
cond &= (n <= M) & (N <= M)
3684
def _pmf(self, k, M, n, N):
3688
return comb(good,k) * comb(bad,N-k) / comb(tot,N)
3689
def _stats(self, M, n, N):
3697
var = m*n*N*(tot-N)*1.0/(tot*tot*(tot-1))
3698
g1 = (m - n)*(tot-2*N) / (tot-2.0)*sqrt((tot-1.0)/(m*n*N*(tot-N)))
3699
m2, m3, m4, m5 = m**2, m**3, m**4, m**5
3700
n2, n3, n4, n5 = n**2, n**2, n**4, n**5
3701
g2 = m3 - m5 + n*(3*m2-6*m3+m4) + 3*m*n2 - 12*m2*n2 + 8*m3*n2 + n3 \
3702
- 6*m*n3 + 8*m2*n3 + m*n4 - n5 - 6*m3*N + 6*m4*N + 18*m2*n*N \
3703
- 6*m3*n*N + 18*m*n2*N - 24*m2*n2*N - 6*n3*N - 6*m*n3*N \
3704
+ 6*n4*N + N*N*(6*m2 - 6*m3 - 24*m*n + 12*m2*n + 6*n2 + \
3706
return mu, var, g1, g2
3707
def _entropy(self, M, n, N):
3708
k = r_[N-(M-n):min(n,N)+1]
3709
vals = self.pmf(k,M,n,N)
3710
lvals = where(vals==0.0,0.0,log(vals))
3711
return -sum(vals*lvals)
3712
hypergeom = hypergeom_gen(name='hypergeom',longname="A hypergeometric",
3713
shapes="M,n,N", extradoc="""
3715
Hypergeometric distribution
3717
Models drawing objects from a bin.
3718
M is total number of objects, n is total number of Type I objects.
3719
RV counts number of Type I objects in N drawn without replacement from
3724
## Logarithmic (Log-Series), (Series) distribution
3726
class logser_gen(rv_discrete):
3728
return rv._inst._logser(pr,size=(self._size,))
3729
def _argcheck(self, pr):
3730
return (pr > 0) & (pr < 1)
3731
def _pmf(self, k, pr):
3732
return -pr**k * 1.0 / k / log(1-pr)
3733
def _stats(self, pr):
3735
mu = pr / (pr - 1.0) / r
3736
mu2p = -pr / r / (pr-1.0)**2
3738
mu3p = -pr / r * (1.0+pr) / (1.0-pr)**3
3739
mu3 = mu3p - 3*mu*mu2p + 2*mu**3
3742
mu4p = -pr / r * (1.0/(pr-1)**2 - 6*pr/(pr-1)**3 + \
3743
6*pr*pr / (pr-1)**4)
3744
mu4 = mu4p - 4*mu3p*mu + 6*mu2p*mu*mu - 3*mu**4
3745
g2 = mu4 / var**2 - 3.0
3746
return mu, var, g1, g2
3747
logser = logser_gen(a=1,name='logser', longname='A logarithmic',
3748
shapes='pr', extradoc="""
3750
Logarithmic (Log-Series, Series) distribution
3754
## Poisson distribution
3756
class poisson_gen(rv_discrete):
3758
return rand.poisson(mu, self._size)
3759
def _pmf(self, k, mu):
3760
Pk = mu**k * exp(-mu) / arr(special.gamma(k+1))
3762
def _cdf(self, x, mu):
3764
return special.pdtr(k,mu)
3765
def _sf(self, x, mu):
3767
return special.pdtrc(k,mu)
3768
def _ppf(self, q, mu):
3769
vals = ceil(special.pdtrik(q,mu))
3770
temp = special.pdtr(vals-1,mu)
3771
return where((temp >= q), vals1, vals)
3772
def _stats(self, mu):
3774
g1 = 1.0/arr(sqrt(mu))
3776
return mu, var, g1, g2
3777
poisson = poisson_gen(name="poisson", longname='A Poisson',
3778
shapes="mu", extradoc="""
3780
Poisson distribution
3784
## (Planck) Discrete Exponential
3786
class planck_gen(rv_discrete):
3787
def _argcheck(self, lambda_):
3796
return 0 # lambda_ = 0
3797
def _pmf(self, k, lambda_):
3798
fact = (1-exp(-lamba))
3799
return fact*exp(-lambda_(k))
3800
def _cdf(self, x, lambda_):
3802
return 1-exp(-lambda_*(k+1))
3803
def _ppf(self, q, lambda_):
3804
val = ceil(-1.0/lambda_ * log(1-q)-1)
3806
def _stats(self, lambda_):
3807
mu = 1/(exp(lambda_)-1)
3808
var = exp(-lambda_)/(1-exp(-lambda_))**2
3809
g1 = 2*cosh(lambda_/2.0)
3810
g2 = 4+2*cosh(lambda_)
3811
return mu, var, g1, g2
3812
def _entropy(self, lambda_):
3815
return l*exp(-l)/C - log(C)
3816
planck = planck_gen(name='planck',longname='A discrete exponential ',
3820
Planck (Discrete Exponential)
3822
pmf is p(k; b) = (1-exp(-b))*exp(-b*k) for kb >= 0
3826
class boltzmann_gen(rv_discrete):
3827
def _pmf(self, k, lambda_, N):
3828
fact = (1-exp(-lambda_))/(1-exp(-lambda_*N))
3829
return fact*exp(-lambda_(k))
3830
def _cdf(self, x, lambda_, N):
3832
return (1-exp(-lambda_*(k+1)))/(1-exp(-lambda_*N))
3833
def _ppf(self, q, lambda_, N):
3834
qnew = q*(1-exp(-lambda_*N))
3835
val = ceil(-1.0/lambda_ * log(1-qnew)-1)
3837
def _stats(self, lambda_, N):
3839
zN = exp(-lambda_*N)
3840
mu = z/(1.0-z)-N*zN/(1-zN)
3841
var = z/(1.0-z)**2 - N*N*zN/(1-zN)**2
3843
trm2 = (z*trm**2 - N*N*zN)
3844
g1 = z*(1+z)*trm**3 - N**3*zN*(1+zN)
3845
g1 = g1 / trm2**(1.5)
3846
g2 = z*(1+4*z+z*z)*trm**4 - N**4 * zN*(1+4*zN+zN*zN)
3847
g2 = g2 / trm2 / trm2
3848
return mu, var, g1, g2
3850
boltzmann = boltzmann_gen(name='boltzmann',longname='A truncated discrete exponential ',
3854
Boltzmann (Truncated Discrete Exponential)
3856
pmf is p(k; b, N) = (1-exp(-b))*exp(-b*k)/(1-exp(-b*N)) for k=0,..,N-1
3865
class randint_gen(rv_discrete):
3866
def _argcheck(self, min, max):
3870
def _pmf(self, k, min, max):
3871
fact = 1.0 / (max - min)
3873
def _cdf(self, x, min, max):
3875
return (k-min+1)*1.0/(max-min)
3876
def _ppf(self, q, min, max):
3877
val = ceil(q*(max-min)+min)
3879
def _stats(self, min, max):
3880
m2, m1 = arr(max), arr(min)
3881
mu = (m2 + m1 - 1.0) / 2
3883
var = (d-1)*(d+1.0)/12.0
3885
g2 = -6.0/5.0*(d*d+1.0)/(d-1.0)*(d+1.0)
3886
return mu, var, g1, g2
3887
def rvs(self, min, max=None, size=None):
3888
"""An array of *size* random integers >= min and < max.
3890
If max is None, then range is >=0 and < min
3895
U = random(size=size)
3896
val = floor((max-min)*U + min)
3897
return arr(val).astype(Num.Int)
3898
def _entropy(self, min, max):
3900
randint = randint_gen(name='randint',longname='A discrete uniform '\
3901
'(random integer)', shapes="min,max",
3906
Random integers >=min and <max.
3913
class zipf_gen(rv_discrete):
3915
return rv._inst._Zipf(a, size=(self._size,))
3916
def _argcheck(self, a):
3918
def _pmf(self, k, a):
3919
Pk = 1.0 / arr(special.zeta(a,1) * k**a)
3921
def _munp(self, n, a):
3922
return special.zeta(a-n,1) / special.zeta(a,1)
3923
def _stats(self, a):
3925
fac = arr(special.zeta(a,1))
3926
mu = special.zeta(a-1.0,1)/fac
3927
mu2p = special.zeta(a-2.0,1)/fac
3929
mu3p = special.zeta(a-3.0,1)/fac
3930
mu3 = mu3p - 3*mu*mu2p + 2*mu**3
3931
g1 = mu3 / arr(var**1.5)
3933
mu4p = special.zeta(a-4.0,1)/fac
3935
mu4 = mu4p - 4*mu3p*mu + 6*mu2p*mu*mu - 3*mu**4
3936
g2 = mu4 / arr(var**2) - 3.0
3937
return mu, var, g1, g2
3938
zipf = zipf_gen(a=1,name='zipf', longname='A Zipf',
3939
shapes="a", extradoc="""
3946
# Discrete Laplacian
3948
class dlaplace_gen(rv_discrete):
3949
def _pmf(self, k, a):
3950
return tanh(a/2.0)*exp(-a*abs(k))
3951
def _cdf(self, x, a):
3955
return where(ind, 1.0-exp(-a*k)/const, exp(a*(k+1))/const)
3956
def _ppf(self, q, a):
3957
const = 1.0/(1+exp(-a))
3960
return ceil(1.0/a*where(ind, log(q*cons2)-1, -log((1-q)*cons2)))
3961
def _stats(self, a):
3966
mu2 = 2* (e2a + ea) / (1-ea)**3.0
3967
mu4 = 2* (e4a + 11*e3a + 11*e2a + ea) / (1-ea)**5.0
3968
return 0.0, mu2, 0.0, mu4 / mu2**2.0 - 3
3969
def _entropy(self, a):
3970
return a / sinh(a) - log(tanh(a/2.0))
3971
dlaplace = dlaplace_gen(a=-inf,
3972
name='dlaplace', longname='A discrete Laplacian',
3973
shapes="a", extradoc="""
3975
Discrete Laplacian distribution.
added
removed
Lib/stats/distributions.py
