← Back to branch summary

~ubuntu-branches/ubuntu/karmic/python-scipy/karmic

« back to all changes in this revision

Viewing changes to scipy/sandbox/montecarlo/montecarlo.py

  • Committer: Bazaar Package Importer
  • Author(s): Ondrej Certik
  • Date: 2008-06-16 22:58:01 UTC
  • mfrom: (2.1.24 intrepid)
  • Revision ID: james.westby@ubuntu.com-20080616225801-irdhrpcwiocfbcmt
Tags: 0.6.0-12
http://bugs.debian.org/489149
* The description updated to match the current SciPy (Closes: #489149).
* Standards-Version bumped to 3.8.0 (no action needed)
* Build-Depends: netcdf-dev changed to libnetcdf-dev

Show diffs side-by-side

added added

removed removed

Lines of Context:
scipy/sandbox/montecarlo/montecarlo.py
 
1
# montecarlo.py: Routines for Monte Carlo simulation
 
2
 
 
3
# Copyright: Ed Schofield, 2005-2006
 
4
# License: BSD-style (see LICENSE.txt at root of scipy tree)
 
5
 
 
6
__author__ = "Ed Schofield"
 
7
 
 
8
from __future__ import division
 
9
import numpy
 
10
import scipy
 
11
from scipy.sandbox.montecarlo._intsampler import _intsampler
 
12
 
 
13
 
 
14
class genericsampler(object):
 
15
    """A base class for other samplers.
 
16
    """
 
17
    def __init__(self):
 
18
        self.format = self.__class__.__name__[:3]
 
19
        if self.__class__.__name__ == "genericsampler":
 
20
            raise TypeError, "this class cannot be instantiated directly"
 
21
    
 
22
    def seed(self, myseed):
 
23
        """Initializes the sampler with a given seed.  If the seed is 0, uses a
 
24
        seed from /dev/urandom or the system time.
 
25
        """
 
26
        self.sampler.seed(myseed)
 
27
 
 
28
 
 
29
class intsampler(genericsampler):
 
30
    """A class that samples objects from a given discrete distribution.
 
31
    The distribution is defined on an integer-valued sample space and
 
32
    specified with a PMF as a list or array like this:
 
33
 
 
34
    >>> table = [10, 15, 20]
 
35
 
 
36
    representing this pmf:
 
37
        x       0       1       2
 
38
        p(x)    10/45   15/45   20/45
 
39
 
 
40
    The output will be something like:
 
41
    >>> sampler = intsampler(table)
 
42
 
 
43
    >> sampler.sample(10)
 
44
    array([c, b, b, b, b, b, c, b, b, b], dtype=object)
 
45
 
 
46
    The class is a thin wrapper around the '_intsampler' class, which
 
47
    uses the 5 table compressed lookup sampler of Marsaglia (2004).  It
 
48
    is very fast, and requires near-constant time in the size of the
 
49
    sample space.  It provides a friendlier interface than _intsampler,
 
50
    allowing use of an integer-valued PMF (unnormalized) and able to
 
51
    return pmf values of the sample points in addition to the sample
 
52
    itself.
 
53
    """
 
54
    def __init__(self, mytable):
 
55
        self.probs = numpy.array(mytable, float)
 
56
        s = self.probs.sum()
 
57
        if s > 0:
 
58
            self.probs /= s
 
59
        else:
 
60
            raise ValueError, "sum of table frequencies must be > 0"
 
61
 
 
62
        self.sampler =  _intsampler(self.probs)
 
63
 
 
64
 
 
65
    def sample(self, size, return_probs=0):
 
66
        """Generates a sample of the given size from the specified
 
67
        discrete distribution, optionally returning the probabilities
 
68
        under the distribution.
 
69
 
 
70
        The optional argument return_probs, 0 by default, has the
 
71
        following meaning:
 
72
            0: don't return pmf values at each sample point
 
73
            1: return pmf values at each sample point
 
74
            2: return log pmf values at each sample point
 
75
        """
 
76
        sample = self.sampler.sample(size)
 
77
        if return_probs == 0:
 
78
            return sample
 
79
        elif return_probs > 0:
 
80
            # More fancy indexing:
 
81
            sampleprobs = self.probs[sample]
 
82
            if return_probs == 1:
 
83
                return (sample, sampleprobs)
 
84
            elif return_probs == 2:
 
85
                return (sample, numpy.log(sampleprobs))
 
86
 
 
87
 
 
88
 
 
89
class dictsampler(genericsampler):
 
90
    """A class that samples objects from a given discrete distribution.
 
91
    The distribution is specified as a dictionary representing a PMF
 
92
    like this:
 
93
    >>> table = {}
 
94
    >>> table['a'] = 10
 
95
    >>> table['b'] = 15
 
96
    >>> table['c'] = 20
 
97
 
 
98
    The output will be something like:
 
99
    >>> sampler = dictsampler(table)
 
100
 
 
101
    >> sampler.sample(10)
 
102
    array([c, b, b, b, b, b, c, b, b, b], dtype=object)
 
103
 
 
104
    giving a sample from this distribution:
 
105
    x       'a'       'b'       'c'
 
106
    p(x)   10/180   150/180   20/180
 
107
 
 
108
    The function uses the constant-time 'intsampler' class, and should be
 
109
    very fast.
 
110
    """
 
111
    def __init__(self, mydict):
 
112
        # We can't use this:
 
113
        #   self.labels = numpy.array(mydict.keys(), object)
 
114
        # since numpy's construction of object arrays is dodgy.  Instead,
 
115
        # create an empty object array and fill it:
 
116
        self.labels = numpy.empty(len(mydict), dtype=object)
 
117
        for i, label in enumerate(mydict):
 
118
            self.labels[i] = label
 
119
        self.probs = numpy.array(mydict.values(), float)
 
120
        s = self.probs.sum()
 
121
        if s > 0:
 
122
            self.probs /= s
 
123
        else:
 
124
            raise ValueError, "sum of table frequencies must be > 0"
 
125
 
 
126
        self.sampler =  intsampler(self.probs)
 
127
    
 
128
    
 
129
    def sample(self, size, return_probs=0):
 
130
        """Generates a sample of the given size from the specified
 
131
        discrete distribution, optionally returning the probabilities
 
132
        under the distribution.
 
133
 
 
134
        The optional argument return_probs, 0 by default, has the
 
135
        following meaning:
 
136
            0: don't return pmf values at each sample point
 
137
            1: return pmf values at each sample point
 
138
            2: return log pmf values at each sample point
 
139
        """
 
140
        sampleindices = self.sampler.sample(size)
 
141
        # Fancy indexing with the object array of labels
 
142
        sample = self.labels[sampleindices]
 
143
        if return_probs == 0:
 
144
            return sample
 
145
        elif return_probs > 0:
 
146
            # More fancy indexing:
 
147
            sampleprobs = self.probs[sampleindices]
 
148
            if return_probs == 1:
 
149
                return (sample, sampleprobs)
 
150
            elif return_probs == 2:
 
151
                return (sample, numpy.log(sampleprobs))
 
152
 
 
153
 
 
154
def _test():
 
155
    import doctest
 
156
    doctest.testmod()
 
157
 
 
158
if __name__ == "__main__":
 
159
    _test()