"""The gaussian distribution.
The primitive distribution is given by function gauss(). For (one-dimensional)
gaussian distributions, use the class Normal(mean, stddev), derived from
variate.Variate, q.v., so as to provide probabilities as integrals of gauss(),
with suitable scalings.
Some day I may add support for multi-dimensional gaussian distributions. I'll
use Gaussian(mean, variance) as its interface; mean will be a vector quantity,
variance will be a tensor quantity of the same rank as mean*mean. For the
present, a function Gaussian(mean, variance) is defined, working only for
one-dimensional data, implemented using Normal.
See also gamma.py for Gamma distributions.
Note that a Gaussian distribution always has non-zero probability of delivering
negative values. Many variates are incapable of delivering negative values (for
example, the heights of human beings) so clearly cannot be normally distributed
(though they may be Gamma distributed). However, provided the variate is also
incapable of being zero, we can take the logarithm of the variate (divided, if
necessary, by some suitable unit) to obtain a variate which may be negative or
positive, so may be well modelled by a normal distribution. From this we can
then infer a distribution for the original variate, the `lognormal'
distribution, which does guarantee positive values, yet has roughly the form of
the normal distribution.
From HAKMEM, http://www.inwap.com/pdp10/hbaker/hakmem/random.html:
ITEM 26 (? via Salamin):
A mathematically exact method of generating a Gaussian distribution
from a uniform distribution: let x be uniform on [0,1] and y uniform
on [0, 2 pi], x and y independent. Calculate r = sqrt(-log x). Then
r cos y and r sin y are two independent Gaussian distributed random
numbers.
ITEM 27 (Salamin):
PROBLEM: Generate random unit vectors in N-space uniform on the unit
sphere. SOLUTION: Generate N Gaussian random numbers and normalize
to unit length.
$Id: gauss.py,v 1.2 2005/01/04 23:49:44 eddy Exp $
"""
import math
def gauss(x, e=math.exp, n=(2*math.pi)**.5): return e(-x*x/2)/n
del math
from variate import Variate
class Normal (Variate):
__upinit = Variate.__init__
def __init__(self, mean, stddev):
"""Initialises a (one-dimensional) Normal distribution.
Takes two arguments; mean and stddev, the mean and standard deviation
of the normal distribution. """
scalar = 1
try:
mean.width, mean.best
if callable(mean.evaluate): scalar = None
except AttributeError: pass
try:
stddev.width, stddev.best
if callable(stddev.evaluate): scalar = None
except AttributeError: pass
if scalar:
# ordinary python numeric types, or functional equivalents
def func(val, m=mean, s=stddev, g=gauss):
return g((val-m)/s)/s
else:
# Quantity(), supporting sophistication :^)
def func(val, m=mean, s=stddev, g=gauss):
return ((val-m)/s).evaluate(g)/s
self.__upinit(func, stddev)
self.mean, self.sigma = mean, stddev
import random
def sample(self, g=random.gauss):
return g(self.mean, self.sigma)
del random
def Gaussian(mean, variance):
return Normal(mean, variance**.5)
class logNormal (Variate):
"""Distribution of a variate whose logarithm is normally distributed.
If the variate is X, then there is some constant m for which log(X/m) is
normally distributed with expected value 0. [To find such an m; chose any
value X could take, say x, and find the mean, say q, of log(X/x); then
log(X/x) -q is normally distributed with mean 0; and log(X/x) -q is just
log(X/x.exp(q)) so use m = x.exp(q).]
"""
__upinit = Variate.__init__
def __init__(self, mean, vary):
pass
del Variate
_rcs_log = """
$Log: gauss.py,v $
Revision 1.2 2005/01/04 23:49:44 eddy
Added notes from HAKMEM.
Revision 1.1 2002/10/08 21:35:02 eddy
Initial revision
"""