Changeset 134
- Timestamp:
- 04/03/08 22:14:41 (8 months ago)
- Files:
-
- projects/AsynCluster/trunk/svpmc (moved) (moved from projects/AsynCluster/trunk/pypmc)
- projects/AsynCluster/trunk/svpmc/model.py (copied) (copied from projects/AsynCluster/trunk/pypmc/model.py)
- projects/AsynCluster/trunk/svpmc/pmc.py (copied) (copied from projects/AsynCluster/trunk/pypmc/pmc.py)
- projects/AsynCluster/trunk/svpmc/sample.py (copied) (copied from projects/AsynCluster/trunk/pypmc/sample.py) (2 diffs)
- projects/AsynCluster/trunk/svpmc/test (copied) (copied from projects/AsynCluster/trunk/pypmc/test)
- projects/AsynCluster/trunk/svpmc/test/project.py (modified) (1 diff)
- projects/AsynCluster/trunk/svpmc/test/test_sample.py (modified) (1 diff)
- projects/AsynCluster/trunk/svpmc/test/test_tseries.py (modified) (1 diff)
- projects/AsynCluster/trunk/svpmc/tseries.py (copied) (copied from projects/AsynCluster/trunk/pypmc/tseries.py)
- projects/AsynCluster/trunk/svpmc/__init__.py (copied) (copied from projects/AsynCluster/trunk/pypmc/__init__.py)
Legend:
- Unmodified
- Added
- Removed
- Modified
- Copied
- Moved
projects/AsynCluster/trunk/svpmc/sample.py
r132 r134 22 22 import scipy as s 23 23 from scipy import random, weave 24 from scipy.stats.distributions as dists 24 25 25 26 … … 83 84 I1 = s.greater(V, x[RI,0]).choose(RI, x[RI,1].astype(int)) 84 85 return I0[I1] 86 87 88 class NormalWalk(object): 89 """ 90 What is a good approximation to the inverse of the cumulative normal 91 distribution function? 92 93 Mathematically we can write this as: Find X such that Q(X) = p for 94 any 0 < p < 1. (Note: this computes an upper tail probability.) 95 96 Again, it is not possible to write this as a closed form expression, 97 so we resort to approximations. Because of the symmetry of the normal 98 distribution, we need only consider 0 < p < 0.5. If you have p > 0.5, 99 then apply the algorithm below to q = 1-p, and then negate the value 100 for X obtained. (This approximation is also from Abramowitz and 101 Stegun.) 102 103 t = sqrt[ ln(1/p^2) ] 104 105 c_0 + c_1*t + c_2*t^2 106 X = t - ------------------------------ 107 1 + d_1*t + d_2*t^2 + d_3*t^3 108 109 c_0 = 2.515517 110 c_1 = 0.802853 111 c_2 = 0.010328 112 113 d_1 = 1.432788 114 d_2 = 0.189269 115 d_3 = 0.001308 116 117 See Abramowitz and Stegun; Press, et al. 118 """ 119 def cinv(self, p): 120 """ 121 The inverse of the cumulative distribution function. 122 123 Returns the random variate value(s) occurring with probability less 124 than or equal to the specified probabilit(ies) from 0.0 - 1.0. 125 """ 126 # TODO: Write C-weave code to do this lightning-fast with the above 127 # approximation. 128 if not hasattr(self, 'distObj'): 129 self.distObj = dists.norm() 130 return self.distObj.ppf(p) 131 132 def kernel(self, x, sdev): 133 """ 134 For each current random variate value in I{x}, returns the transition 135 kernel... 136 """ projects/AsynCluster/trunk/svpmc/test/project.py
r133 r134 21 21 22 22 import os.path 23 from pypmc.tseries import TimeSeries23 from svpmc.tseries import TimeSeries 24 24 25 25 def ts(fileName): projects/AsynCluster/trunk/svpmc/test/test_sample.py
r132 r134 17 17 18 18 """ 19 Unit tests for pypmc.sample19 Unit tests for svpmc.sample 20 20 """ 21 21 projects/AsynCluster/trunk/svpmc/test/test_tseries.py
r133 r134 17 17 18 18 """ 19 Unit tests for pypmc.tseries19 Unit tests for svpmc.tseries 20 20 """ 21 21
