"""Modelling random variates.
A random variate's model needs to be able to generate `unpredictably'
(chaotically or pseudo-randomly) from the vairate's distribution; it also needs
to be able to perform the integration needed to determine moments (e.g. the
mean) and the probability of the variate falling in any given range of its
permitted values.
$Id: variate.py,v 1.4 2007/06/03 16:42:54 eddy Exp $
"""
from integrate import Integrator
from study.snake.lazy import Lazy
class Variate (Integrator, Lazy):
__upinit = Integrator.__init__
# Integrator provides between(), beyond() and before().
def __init__(self, func, width=None):
self.__upinit(func, width)
self.__moments = []
def __moment(self, i, s):
g = self.measure(lambda x, j=i: x**j)
return g.before(s) + g.between(s,s) + g.beyond(s)
def moments(self, n):
"""Returns expected values of various powers of the variate.
Single argument, n, will be the length of the tuple returned: its
entry[i] is the expected value of self.sample()**(1+i). Thus the first
entry is the mean, the second is the `mean square' and so on. """
seed = self.sample()
total = self.before(seed) + self.beyond(seed) # self.__moment(0, s)
return tuple(map(lambda i, s=seed, t=total, m=self.__moment: m(i,s)/t,
range(1, 1+n)))
def _lazy_get_mean_(self, ignored):
return self.moments(1)[0]
def _lazy_get_variance_(self, ignored):
self.mean, two = self.moments(2)
return two - self.mean**2
def sample(self):
"""Returns a sample value from this distribution.
Each call returns a fresh value independently of previous calls.
This method should be over-ridden in derived classes.
See RatioGenerator for a generic algorithm for implementing samples.
"""
raise NotImplementedError
from random import random # 0 <= uniform < 1
class Uniform (Variate):
__upinit = Variate.__init__
def __init__(self, lo=1, hi=None):
if hi is None: lo, hi = 0, lo # i.e. lo has default 0, really hi was supplied ...
self.min, self.max, self.width = lo, hi, hi - lo
self.__height = 1./self.width
self.__upinit(self.__p, self.width)
def __p(self, x):
if self.min < x < self.max: return self.__height
return 0 * self.__height # may have units of measurement ...
def sample(self):
return random() * self.width + self.min
class RatioGenerator:
"""Illustration of a technique.
If source() generates samples from a distribution g, ratio = (: f(y)/g(y)
←y :), and (|ratio|source) is bounded below by 0 and above by bound,
calls to the object returned by RatioGenerator(source, ratio, bound) will be
samples from the distribution f. """
def __init__(self, source, ratio, bound):
self.__source, self.__ratio, self.__bound = source, ratio, bound
def __call__(self):
"""The ratio test algorithm."""
while 1:
ans = self.__source()
if self.__ratio(ans) > random() * bound:
return ans