Changeset 137

Show
Ignore:
Timestamp:
04/05/08 13:51:58 (8 months ago)
Author:
edsuom
Message:

New base class 'Weave' for convenience with Weave inline operations; working on MCMC for normal random walk

Files:

Legend:

Unmodified
Added
Removed
Modified
Copied
Moved

  • projects/AsynCluster/trunk/svpmc/sample.py

    r136 r137  
    2121 
    2222import scipy as s 
    23 from scipy import random, weave 
     23from scipy import random 
    2424from scipy.stats import distributions as dists 
    2525 
    26  
    27 RESAMPLE_CODE = r""" 
    28 int ja = Na[0] - 1; 
    29 int jb = Nb[0] - 1; 
    30 int ka, kb; 
    31 while (ja >= 0 && jb >= 0) { 
    32   ka = A1(ja); 
    33   kb = B1(jb); 
    34   X2(ka,1) = kb; 
    35   X2(kb,0) = X2(kb,0) + X2(ka,0) - 1.0; 
    36   if (X2(kb,0) < 1) { 
    37     A1(ja) = kb; 
    38     jb--; 
    39   } else { 
    40     ja--; 
    41   } 
    42 } 
    43 """ 
     26from util import Weaver 
    4427 
    4528 
    46 def resample(W, N=None, logWeights=False): 
     29class Resampler(Weaver): 
    4730    """ 
    48     Returns an index array that resamples I{N} values from some external 
    49     array. Each element of the array will have a probability of being sampled 
    50     (with replacement) that is proportional to its corresponding weight in the 
    51     1-D array I{W}. 
    52      
     31    Call an instance of me with a 1-D array of weights to get an index array 
     32    that resamples values from some external array. Each element of the array 
     33    will have a probability of being sampled (with replacement) that is 
     34    proportional to its corresponding weight in the array. 
     35 
    5336    Uses Walker's alias algorithm as set forth in L. Devroye, Non-Uniform 
    5437    Random Variate Generation, p. 109, 
    55     http://cg.scs.carleton.ca/~luc/rnbookindex.html. 
     38    U{http://cg.scs.carleton.ca/~luc/rnbookindex.html}. 
     39    """ 
    5640 
    57     If N is not specified, it defaults to the length of W, i.e., importance 
    58     resampling. 
     41    def __init__(self, logWeights=False): 
     42        self.logWeights = logWeights 
    5943 
    60     If I{logWeights} is set C{True}, the weights are transformed from log 
    61     to linear. 
    62     """ 
    63     if logWeights: 
    64         W = s.exp(W - W.max()) 
    65     W = W.astype(float) / sum(W) 
    66     I0 = (s.isfinite(W) * s.greater(W, 1E-7)).nonzero()[0] 
    67     K = len(I0) 
    68     if N is None: 
    69         N = K 
    70     if K == 0: 
    71         return [] 
    72     if K == 1: 
    73         return [0]*N             
    74     x = s.zeros((K,2)) 
    75     x[:,0] = K*W[I0] 
    76     a = s.less(x[:,0], 1.0).nonzero()[0] 
    77     b = s.setdiff1d(s.arange(0, K), a) 
    78     weave.inline( 
    79         RESAMPLE_CODE, 
    80         ['x', 'a', 'b'], 
    81         extra_compile_args=['-w']) 
    82     V = s.rand(N) 
    83     RI = s.random.randint(0, K, N) 
    84     I1 = s.greater(V, x[RI,0]).choose(RI, x[RI,1].astype(int)) 
    85     return I0[I1] 
     44    def __call__(self, W, N=None): 
     45        """ 
     46        If N is not specified, it defaults to the length of W, i.e., importance 
     47        resampling. 
     48 
     49        If my I{logWeights} attribute is set C{True}, the weights are 
     50        transformed from log to linear. 
     51        """ 
     52        if self.logWeights: 
     53            W = s.exp(W - W.max()) 
     54        W = W.astype(float) / sum(W) 
     55        I0 = (s.isfinite(W) * s.greater(W, 1E-7)).nonzero()[0] 
     56        K = len(I0) 
     57        if N is None: 
     58            N = K 
     59        if K == 0: 
     60            return [] 
     61        if K == 1: 
     62            return [0]*N             
     63        x = s.zeros((K,2)) 
     64        x[:,0] = K*W[I0] 
     65        a = s.less(x[:,0], 1.0).nonzero()[0] 
     66        x = self.c('x', 'a', 'b', x=x, a=a, b=s.setdiff1d(s.arange(0, K), a)) 
     67        V = s.rand(N) 
     68        RI = s.random.randint(0, K, N) 
     69        I1 = s.greater(V, x[RI,0]).choose(RI, x[RI,1].astype(int)) 
     70        return I0[I1] 
    8671 
    8772 
    88 class NormalWalk(object): 
     73class NormalWalk(Weaver): 
    8974    """ 
    90     I run an array of values through a random walk within a zero-mean, unit 
    91     variance normal distribution. Instantiate me with the dimensions of the 
    92     array as one or more arguments. 
     75    I run a 2-D array of values through a random walk within a zero-mean, unit 
     76    variance normal distribution. 
     77     
     78    Call an instance of me with an array I{w} of fractional wiggle values in 
     79    the range 0.0 to +1.0, with one array element for each of my random 
     80    walkers. The elements of my value array will be perturbed accordingly and 
     81    updated array returned. 
    9382    """ 
    94     def __init__(self, *dims): 
    95         self._dims = dims 
    96         self.distObj = dists.norm() 
    97         self.X = s.rand(*dims) 
    98  
    99     def cinv(self, p): 
    100         """ 
    101         The inverse of the cumulative distribution function. 
    102  
    103         Returns the random variate value(s) occurring with probability less 
    104         than or equal to the specified probabilit(ies) from 0.0 to +1.0. 
    105         """ 
    106         # TODO: Write C-weave code to do this lightning-fast with the following 
    107         # approximation. Because of the symmetry of the normal distribution, we 
    108         # need only consider 0 < p < 0.5. If you have p > 0.5, then apply the 
    109         # algorithm below to q = 1-p, and then negate the value for X obtained. 
    110         # 
    111         # t = sqrt[ ln(1/p^2) ] 
    112         # 
    113         #             c_0 + c_1*t + c_2*t^2 
    114         # X = t -  ------------------------------ 
    115         #          1 + d_1*t + d_2*t^2 + d_3*t^3 
    116         # 
    117         # c_0 = 2.515517 
    118         # c_1 = 0.802853 
    119         # c_2 = 0.010328 
    120         #  
    121         # d_1 = 1.432788 
    122         # d_2 = 0.189269 
    123         # d_3 = 0.001308 
    124         # 
    125         # See Abramowitz and Stegun; Press, et al. 
    126         return self.distObj.ppf(p) 
     83    def __init__(self, rows, cols): 
     84        self.x = s.rand(rows, cols) 
    12785     
    128     def step(self, W): 
    129         """ 
    130         Call this method with an array I{W} of fractional wiggle values in the 
    131         range 0.0 to +1.0, with one array element for each of my random 
    132         walkers. The elements of my value array will be perturbed accordingly. 
    133  
    134         Two like-dimensioned arrays are returned, one with the target values 
    135         after the jumps and another containing the probability densities of the 
    136         jumps. 
    137         """ 
    138         if s.shape(W) != self._dims: 
     86    def __call__(self, w): 
     87        if s.shape(w) != self.w.shape: 
    13988            raise ValueError( 
    14089                "Wiggle array must have same dimensions as walker array") 
    141         # TODO: Make work with 2-D array 
    142         P = s.empty_like(W) 
    143         for k, dx in enumerate(0.5*W): 
    144             distObj = dists.truncnorm(self.X[k]-dx, self.X[k]+dx) 
    145             self.X[k] = distObj.rvs(1)[0] 
    146             P[k] = distObj.pdf(self.X[k]) 
    147         return self.X, P 
    148      
     90        return self.c('x', 'w', w=w) 
     91         
    14992             
    15093         

  • projects/AsynCluster/trunk/svpmc/test/test_sample.py

    r136 r137  
    2828 
    2929class Test_Resample(util.TestCase): 
     30    def setUp(self): 
     31        self.resample = sample.Resampler() 
     32     
    3033    def test_oneElement(self): 
    31         I = sample.resample(s.array([1])) 
     34        I = self.resample(s.array([1])) 
    3235        self.failUnlessEqual(I[0], 0) 
    3336 
    3437    def test_uniform(self): 
    3538        W = s.array([1,1,1]) 
    36         I = sample.resample(W, 10000) 
     39        I = self.resample(W, 10000) 
    3740        counts = [s.equal(I, k).sum() for k in xrange(len(W))] 
    3841        self.failUnless(stats.chisquare(counts)[1] > 0.05) 
    …… 
    4043    def test_basic(self): 
    4144        W = s.array([1,2,3]) 
    42         I = sample.resample(W, 10000) 
     45        I = self.resample(W, 10000) 
    4346        counts = [s.equal(I, k).sum() for k in xrange(len(W))] 
    4447        counts = s.array(counts) * s.array([2, 1, 2./3]) 
    …… 
    5356        samples, AR = [], [] 
    5457        for k in xrange(10): 
    55             I = sample.resample(W) 
     58            I = self.resample(W) 
    5659            samples.append(X[I]) 
    5760            AR.append(len(s.unique(I))) 
    …… 
    6669class Test_NormalWalk(util.TestCase): 
    6770    def test_univariate(self): 
    68         walker = sample.NormalWalk(1) 
    69         W = s.array([0.2]) 
    70         X = s.array([walker.step(W)[0][0] for k in xrange(50000)]) 
     71        walker = sample.NormalWalk(1, 1) 
     72        w = s.array([[0.2]]) 
     73        x = s.array([walker(w)[0][0] for k in xrange(50000)]) 
    7174        fig = self.fig 
    7275        sp = fig.add_subplot(211) 
    73         sp.plot(X) 
     76        sp.plot(x) 
    7477        sp = fig.add_subplot(212) 
    75         sp.hist(X, bins=20) 
     78        sp.hist(x, bins=20) 
    7679 
    7780