Changeset 168

Show
Ignore:
Timestamp:
04/30/08 20:13:51 (7 months ago)
Author:
edsuom
Message:

Got model unit testing working again; more complete unit testing of params; general cleanup

Files:

Legend:

Unmodified
Added
Removed
Modified
Copied
Moved

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

    r167 r168  
    2020""" 
    2121 
    22 import os.path, re 
    2322import scipy as s 
    2423from scipy import linalg 
    25 from scipy.stats import distributions 
    2624 
    2725from twisted.internet import defer, reactor 
    …… 
    2927from asyncluster.master import jobs 
    3028from twisted_goodies.pybywire.params import Parameterized 
    31 from twisted_goodies.pybywire.pack import Packer, Unpacker 
    32  
    33 import tseries, sample, weave, params 
     29from twisted_goodies.pybywire import pack 
     30 
     31import weave 
    3432 
    3533 
    …… 
    136134                # Unpack string result from remote call 
    137135                result = list(pack.Unpacker(result)) 
    138             for k, name in enumerate(('Lx', 'z', 'h')): 
     136            for k, name in enumerate(('Lx', 'h', 'z')): 
    139137                setattr(paramContainer, name, result[k]) 
    140138            return paramContainer 
    141          
     139 
    142140        if (remote or self.remoteMode) and not local: 
    143141            d = self.client.run( 
    …… 
    393391            1. The value of the log-likelihood. 
    394392 
    395             2. The IID normal variates underlying the log-volatilities. 
    396  
    397             3. The last sample's multivariate log-volatility value. 
     393            2. The last sample's multivariate log-volatility value. 
     394 
     395            3. The IID normal variates underlying the log-volatilities. 
     396 
    398397         
    399398        The returned likelihood does not consider the prior probability of the 

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

    r166 r168  
    290290        deviation width of a normal distribution.) 
    291291        """ 
    292         return self.rdsObj.rvs(1)[0] 
     292        return wiggle * self.rdsObj.rvs(1)[0] 
    293293 
    294294    def pJump(self, x, wiggle): 
    …… 
    318318 
    319319    """ 
     320    attempts = 20 
     321    typeChecking = True 
     322     
    320323    def __init__(self, p, n): 
    321324        self.p = p 
    …… 
    359362                "Specify a positive wiggle value, '%s' obj invalid" \ 
    360363                % str(wiggle)) 
    361         if not isinstance(paramContainer, ParameterContainer): 
     364        if self.typeChecking and \ 
     365               not isinstance(paramContainer, ParameterContainer): 
    362366            raise TypeError( 
    363367                "You must supply a ParameterContainer as the first "+\ 
    …… 
    376380                raise ValueError("Failed to generate a valid proposal") 
    377381            Lp += s.log(pPriors).sum() 
    378             Lj += s.log(priorFlexArray.call('pJump', jumps)).sum() 
     382            Lj += s.log(priorFlexArray.call('pJump', jumps, wiggle)).sum() 
    379383            setattr(newVersion, name, paramArray) 
    380384        newVersion.Lp = Lp 

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

    r166 r168  
    1717 
    1818""" 
    19 Population Monte Carlo sampling. 
     19Population Monte Carlo sampling 
    2020""" 
    2121 
    …… 
    2929import params 
    3030from sample import resample 
     31 
     32 
     33class Allocator(object): 
     34    """ 
     35    Construct me with a population size I{N} and a sequence I{V} of jump 
     36    deviations. 
     37    """ 
     38    def __init__(self, N, V): 
     39        self.N = N 
     40        self.V = V 
     41        self.P = len(V) 
     42        self.rMin = max([2, int(round(0.01 * N))]) 
     43        self.rMax = N - self.rMin*(self.P-1) 
     44        self.updateAllocations() 
     45     
     46    def subsetIndex(self, k): 
     47        """ 
     48        Returns a subset index for population members that correspond to the 
     49        jump deviation for the supplied index I{k}. 
     50        """ 
     51        I = sampleWOR(self.Is, self.R[k]) 
     52        self.Is = s.setdiff1d(self.Is, I) 
     53        return I 
     54     
     55    def assembler(self, resultList): 
     56        """ 
     57        Call this method with a reference to an empty list that will hold 
     58        assembled results. 
     59 
     60        For each of my jump deviations, the method will yield the deviation 
     61        value, a 1-D array of indices for the subset of my population allocated 
     62        to that deviation value, and a fresh deferred already fired with 
     63        C{None}. You must add a callback to the deferred for each iteration 
     64        that returns one or more 1-D arrays of the same length as the yielded 
     65        index array. 
     66         
     67        Each returned array for each subset will be assembled back into a 
     68        single population-size array and placed into the supplied list, in 
     69        order. 
     70        """ 
     71        def gotResults(results, I): 
     72            if not resultList: 
     73                for result in results: 
     74                    if isinstance(result, (int, float)): 
     75                        resultList.append(s.empty(self.N)) 
     76                    else: 
     77                        resultList.append(params.FlexArray(self.N)) 
     78            for k, result in enumerate(results): 
     79                resultList[k][I] = result 
     80         
     81        if resultList != []: 
     82            raise ValueError( 
     83                "You must supply an empty list to hold your assembled results") 
     84        self.II = [None] * self.P 
     85        for k in s.flipud(s.argsort(self.R)): 
     86            I = self.II[k] = self.subsetIndex(k) 
     87            d = defer.succeed(None) 
     88            yield self.V[k], I, d 
     89            d.addCallback(gotResults, I) 
     90 
     91    def updateAllocations(self, I=None): 
     92        """ 
     93        Rescales my allocations based on the supplied 1-D array of resampling 
     94        Highly successful allocations will have a large portion of the next 
     95        allocation even if they last operated on only a few members of the 
     96        population. 
     97        """ 
     98        if I is None: 
     99            R = self.N * s.ones(self.P) / self.P 
     100        elif hasattr(self, 'II'): 
     101            R = s.array( 
     102                [sum([x in self.II[k] for x in I]) for k in xrange(self.P)]) 
     103            R = R.astype(float) / R.sum() 
     104        else: 
     105            raise AttributeError( 
     106                "You can only update with an index array after completion "+\ 
     107                "of an assembler run") 
     108 
     109        # Turn the scaled array into an array of subsample sizes, with a 
     110        # minimum size of two apiece. 
     111        R = s.clip(s.round_(self.N*R), self.rMin, self.rMax).astype(int) 
     112        # Twiddle the biggest one as needed to keep sum = Number of members 
     113        R[s.argmax(R)] +=  - sum(R) 
     114        # Replace the old list and subset index 
     115        self.R = R 
     116        self.Is = s.arange(self.N) 
    31117 
    32118 
    …… 
    157243            X = XP[I] 
    158244            allocator.updateAllocations(I) 
    159         self.recorder.done() 
    160  
    161  
    162 class Allocator(object): 
    163     """ 
    164     Construct me with a population size I{N} and a sequence I{V} of jump 
    165     deviations. 
    166     """ 
    167     def __init__(self, N, V): 
    168         self.N = N 
    169         self.V = V 
    170         self.P = len(V) 
    171         self.rMin = max([2, int(round(0.01 * N))]) 
    172         self.rMax = N - self.rMin*(self.P-1) 
    173         self.updateAllocations() 
    174      
    175     def subsetIndex(self, k): 
    176         """ 
    177         Returns a subset index for population members that correspond to the 
    178         jump deviation for the supplied index I{k}. 
    179         """ 
    180         I = sampleWOR(self.Is, self.R[k]) 
    181         self.Is = s.setdiff1d(self.Is, I) 
    182         return I 
    183      
    184     def assembler(self, resultList): 
    185         """ 
    186         Call this method with a reference to an empty list that will hold 
    187         assembled results. 
    188  
    189         For each of my jump deviations, the method will yield the deviation 
    190         value, a 1-D array of indices for the subset of my population allocated 
    191         to that deviation value, and a fresh deferred already fired with 
    192         C{None}. You must add a callback to the deferred for each iteration 
    193         that returns one or more 1-D arrays of the same length as the yielded 
    194         index array. 
    195          
    196         Each returned array for each subset will be assembled back into a 
    197         single population-size array and placed into the supplied list, in 
    198         order. 
    199         """ 
    200         def gotResults(results, I): 
    201             if not resultList: 
    202                 for result in results: 
    203                     if isinstance(result, (int, float)): 
    204                         resultList.append(s.empty(self.N)) 
    205                     else: 
    206                         resultList.append(params.FlexArray(self.N)) 
    207             for k, result in enumerate(results): 
    208                 resultList[k][I] = result 
    209          
    210         if resultList != []: 
    211             raise ValueError( 
    212                 "You must supply an empty list to hold your assembled results") 
    213         self.II = [None] * self.P 
    214         for k in s.flipud(s.argsort(self.R)): 
    215             I = self.II[k] = self.subsetIndex(k) 
    216             d = defer.succeed(None) 
    217             yield self.V[k], I, d 
    218             d.addCallback(gotResults, I) 
    219  
    220     def updateAllocations(self, I=None): 
    221         """ 
    222         Rescales my allocations based on the supplied 1-D array of resampling 
    223         Highly successful allocations will have a large portion of the next 
    224         allocation even if they last operated on only a few members of the 
    225         population. 
    226         """ 
    227         if I is None: 
    228             R = self.N * s.ones(self.P) / self.P 
    229         elif hasattr(self, 'II'): 
    230             R = s.array( 
    231                 [sum([x in self.II[k] for x in I]) for k in xrange(self.P)]) 
    232             R = R.astype(float) / R.sum() 
    233         else: 
    234             raise AttributeError( 
    235                 "You can only update with an index array after completion "+\ 
    236                 "of an assembler run") 
    237  
    238         # Turn the scaled array into an array of subsample sizes, with a 
    239         # minimum size of two apiece. 
    240         R = s.clip(s.round_(self.N*R), self.rMin, self.rMax).astype(int) 
    241         # Twiddle the biggest one as needed to keep sum = Number of members 
    242         R[s.argmax(R)] +=  - sum(R) 
    243         # Replace the old list and subset index 
    244         self.R = R 
    245         self.Is = s.arange(self.N) 
    246  
     245        self.pm.done() 
     246 
     247 

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

    r167 r168  
    240240        var[self.iteration, :] = parameters.Lx.astype('f') 
    241241        # Done writing this iteration of the CDF record 
     242        self.cdf.sync() 
    242243        self.iteration += 1 
    243244 
    …… 
    247248        """ 
    248249        self.cdf.close() 
     250 

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

    r145 r168  
    3939  } 
    4040} 
    41  
    42  
    43  
    44 // NormalWalk.step 
    45 // 
    46 // Supplied variables 
    47 // ---------------------------------------------------------------------------- 
    48 // x        2-D array of random walker values 
    49 // p        2-D array of log-probabilities for each value in x 
    50 // m        Oversampling ratio 
    51 // wiggle   Wiggle (0-1) for scaling the uniform random walk proposals 
    52  
    53 int i, j, k; 
    54 double xp, dp, p_xp; 
    55  
    56 // For each row... 
    57 for(i=0; i<Nx[0]; i++) { 
    58   // For each column... 
    59   for(j=0; j<Nx[1]; j++) { 
    60     // For each oversampling iteration... 
    61     for(k=0; k<m; k++) { 
    62       // Generate a uniform, symmetric, random walk proposal 
    63       xp = X2(i,j) + wiggle * (drand48() - 0.5); 
    64       // Do the ratio in log space 
    65       p_xp = -0.5 * pow(xp, 2); 
    66       dp = p_xp - P2(i,j); 
    67       // The uniform random variate thus needs to be in logspace as well, but 
    68       // it's not always even needed. Thus the log operation is done just once, 
    69       // and only some of the time. 
    70       if(dp > 0  || log(drand48()) < dp) { 
    71         // Update current value and its log probability when proposal accepted 
    72         X2(i,j) = xp; 
    73         P2(i,j) = p_xp; 
    74       } 
    75     } 
    76   } 
    77 } 
    78  
    79  
    80  
    81 // NormalWalk.correlate 
    82 // 
    83 // Supplied variables 
    84 // ---------------------------------------------------------------------------- 
    85 // y    2-D array of correlated random variates (to be overwritten) 
    86 // x    2-D array of independent random variates 
    87 // r    2-D array of cross-correlations between values in each column of x 
    88  
    89 int i, j, k; 
    90 double sum; 
    91  
    92 // For each column of independent variates to be correlated... 
    93 for(i=0; i<Nx[1]; i++) { 
    94   // For each of the values in the column... 
    95   for(j=0; j<Nx[0]; j++) { 
    96     // The matrix product of the cross-correlations and the independent 
    97     // variates... 
    98     sum = 0; 
    99     for(k=0; k<Nx[0]; k++) 
    100       sum += R2(j,k) * X2(k,i); 
    101     // ...is the correlated output 
    102     Y2(j,i) = sum; 
    103   } 
    104 } 

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

    r164 r168  
    2121 
    2222import scipy as s 
    23 from scipy import random, linalg 
     23from scipy import random 
    2424 
    2525from weave import Weaver 
    …… 
    6868        I1 = s.greater(V, x[RI,0]).choose(RI, x[RI,1].astype(int)) 
    6969        return I0[I1] 
    70  
    71  
    72 class NormalWalk(Weaver): 
    73     """ 
    74     I run a 2-D array of values through a random walk within a zero-mean, unit 
    75     variance normal distribution. 
    76  
    77     You can call me to access my current array of values, with no argument to 
    78     get a copy of the array or with a conforming-dimension array of new values. 
    79     """ 
    80     m = 10 
    81      
    82     def __init__(self, rows, cols): 
    83         self.x = s.randn(rows, cols) 
    84         self.p = s.zeros_like(self.x) 
    85  
    86     def __call__(self, x=None): 
    87         if x is None: 
    88             return self.x.copy() 
    89         s1, s0 = s.shape(x), self.x.shape 
    90         if len(s1) == 1 or (len(s1) == 2 and s1[1] == 1): 
    91             # Resize doesn't seem to work in all cases, so just build another 
    92             # single-row array 
    93             x = s.row_stack([x]) 
    94             s1 = (1, s1[0]) 
    95         if s1 != s0: 
    96             raise ValueError( 
    97                 "You can only update me with an equivalent-dimensioned array") 
    98         self.x = x.copy() 
    99      
    100     def step(self, wiggle): 
    101         """ 
    102         Call with a fractional wiggle value in the range 0.0 to +1.0. The 
    103         elements of my value array will be perturbed via a random walk MCMC 
    104         step with a symmetric, uniform proposal over +/- half the wiggle value. 
    105  
    106         A reference to the updated array is returned. 
    107         """ 
    108         self.inline('x', 'p', 'm', 'wiggle') 
    109         return self.x 
    110  
    111     def covarMatrix(self, correlations): 
    112         """ 
    113         Returns a covariance matrix from the vector or sequence of 
    114         I{correlations} supplied in the form used by L{correlate}. 
    115  
    116         Assumes that the independent components have unit variance, as is the 
    117         case with my random walkers C{x}. 
    118         """ 
    119         i = 0 
    120         n = self.x.shape[0] 
    121         cv = s.ones((n, n)) 
    122         for j in xrange(n - 1): 
    123             for k in xrange(j+1, n): 
    124                 cv[j, k] = cv[k, j] = correlations[i] 
    125                 i += 1 
    126         return cv 
    127      
    128     def correlate(self, correlations=None): 
    129         """ 
    130         Returns a version of my current values correlated along the first axis, 
    131         given the supplied vector or sequence of concurrent I{correlations}, 
    132         the elements of which are in the following form:: 
    133  
    134             [r01, r02,..., r0n, r12,..., r1n,..] 
    135  
    136         The values in each column of the result are correlated. The values from 
    137         one column to the next remain independent. 
    138  
    139         The correlations vector can be omitted after the first call to this 
    140         method of a given instance of me. In such case, the last correlation 
    141         vector supplied will be re-used. 
    142         """ 
    143         if correlations is not None: 
    144             self.y = s.empty_like(self.x) 
    145             if correlations: 
    146                 # Generate the correlation matrix and save in case not supplied 
    147                 # next time 
    148                 self.r = linalg.cholesky( 
    149                     self.covarMatrix(correlations), lower=True) 
    150             else: 
    151                 self.r = None 
    152              
    153         n = self.x.shape[0] 
    154         if n == 1: 
    155             return self.x.copy() 
    156         if not hasattr(self, 'r') or s.shape(self.r) != (n, n): 
    157             raise ValueError( 
    158                 "Must have a [%d x %d] cross-correlation matrix defined" \ 
    159                 % (n, n)) 
    160         y = s.empty_like(self.x) 
    161         self.inline('y', 'x', 'r') 
    162         return y 
    163          
    164          
    165          

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

    r167 r168  
    3333 
    3434 
    35 VERBOSE = True 
     35VERBOSE = util.Mock.verbose() 
    3636 
    3737 
    …… 
    4545 
    4646    def likelihood(self, paramContainer, sigma): 
    47         return likelihoodFunc(paramContainer, self.drift) 
     47        h = s.zeros(self.p) 
     48        z = s.zeros(self.n) 
     49        return [likelihoodFunc(paramContainer, self.drift), h, z] 
    4850 
    4951 
    …… 
    8082        self.JobClient = model.jobs.JobClient 
    8183        model.jobs.JobClient = Mock_Client 
    82         self.mgr = model.ModelManager(self.model) 
     84        self.mgr = model.ModelManager( 
     85            util.PARAM_NAMES, s.zeros((self.model.p, self.model.n))) 
     86        self.mgr.modelObj = self.model 
    8387 
    8488    def tearDown(self): 
    …… 
    8892        def gotResult(result): 
    8993            self.failUnlessEqual( 
    90                 result[0], likelihoodFunc(paramContainer, 1.5)) 
     94                result.Lx, likelihoodFunc(paramContainer, 1.5)) 
    9195            self.nextCall(Mock_Model.calls) 
    9296            self.nextCall('likelihood') 
    …… 
    9498        Mock_Model.clearCalls() 
    9599        paramContainer = util.Mock_ParameterContainer() 
    96         d = self.mgr.call('likelihood', paramContainer, 1.0) 
     100        d = self.mgr.likelihood(paramContainer, 1.0) 
    97101        d.addCallback(gotResult) 
    98102        return d 
    …… 
    100104    def test_setRemoteMode_generic(self): 
    101105        def gotResult(result): 
    102             self.failUnlessElementsEqual(result, s.exp(-X**2)) 
     106            self.failUnlessEqual( 
     107                result.Lx, likelihoodFunc(paramContainer, 1.5)) 
    103108            self.nextCall(Mock_Client.calls) 
    104109            self.nextCall( 
    …… 
    107112                'update', ('newModel', self.model)) 
    108113            self.nextCall( 
    109                 'run', ('callModel', 'likelihood', pack.Packer(X)())) 
     114                'run', ('likelihood', paramContainer, 1.0)) 
    110115            for call in Mock_Client.calls: 
    111116                if call[0].endswith('.update') and call[1][0] == 'parallow': 
    …… 
    118123             
    119124        Mock_Client.clearCalls() 
    120         X = s.array([0.0, 0.1]) 
     125        paramContainer = util.Mock_ParameterContainer() 
    121126        d = self.mgr.setRemoteMode(None) 
    122         d.addCallback(lambda _: self.mgr.likelihood(X)) 
     127        d.addCallback(lambda _: self.mgr.likelihood(paramContainer, 1.0)) 
    123128        d.addCallback(gotResult) 
    124129        return d 
    …… 
    142147        Instantiate a model to be tested and set its parameters 
    143148        """ 
    144         kw['paramNames'] = ['a', 'b', 'cr', 'd', 'e', 'fs', 'fr', 'g'] 
     149        kw['paramNames'] = util.PARAM_NAMES 
    145150        for name, value in self.params.iteritems(): 
    146151            if name != 'y': 
    …… 
    325330        self.failUnless(percent < 40) 
    326331 
    327     def _check_normality(self, x): 
    328         p = stats.normaltest(x)[1] 
    329         if VERBOSE: 
    330             fig = self.fig 
    331             sp = fig.add_subplot(211) 
    332             sp.plot(x) 
    333             sp = fig.add_subplot(212) 
    334             sp.hist(x, bins=20) 
    335         self.failUnless( 
    336             p > 0.05, "Deviation from normality with significance p=%f" % p) 
    337  
    338332    def test_zProposal(self): 
    339333        N = 1000 
    …… 
    346340            z1[k] = z[0,0] 
    347341            z2[k] = z[1,1] 
    348         self._check_normality(z1[0.7*N:]) 
    349         self._check_normality(z2[0.7*N:]) 
     342        self.failUnlessNormal(z1[0.7*N:]) 
     343        self.failUnlessNormal(z2[0.7*N:]) 
    350344        self._check_uncorr(z1, z2) 
    351345 

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

    r167 r168  
    2020""" 
    2121 
     22import sys 
    2223import scipy as s 
    2324 
    …… 
    207208 
    208209class Test_PriorContainer(util.TestCase): 
    209     p = 3 
    210     shapes = [(p,), (p,p), (0.5*(p**2+p),)  ] 
    211     scales = [ 0.1,  0.2,   0.3             ] 
     210    MPC = util.Mock_ParameterContainer 
    212211     
    213212    def setUp(self): 
    214         self.priorContainer = params.PriorContainer(3, 100) 
    215         self.priorContainer.paramNames = ['a', 'b', 'c'] 
     213        self.priorContainer = params.PriorContainer(self.MPC.p, self.MPC.N) 
     214        self.priorContainer.typeChecking = False 
     215        self.priorContainer.paramNames = self.MPC.paramNames 
    216216        kw = {'dname':'norm', 'loc':0} 
    217         for i , name in enumerate(self.priorContainer.paramNames): 
    218             shape = self.shapes[i] 
    219             kw['scale'] = self.scales[i] 
     217        for i, name in enumerate(self.priorContainer.paramNames): 
     218            shape = self.MPC.paramShapes[i] 
     219            kw['scale'] = 0.1 * (i+1) 
    220220            priorArray = self.priorContainer.priorArray(name, *shape) 
    221221            for j in xrange(int(shape[0])): 
    …… 
    229229        for k, name in enumerate(self.priorContainer.paramNames): 
    230230            priorArray = getattr(self.priorContainer, name) 
    231             self.failUnlessEqual(priorArray.shape, self.shapes[k]) 
    232  
    233     def test_new(self): 
    234         paramContainer = self.priorContainer.new() 
     231            self.failUnlessEqual(priorArray.shape, self.MPC.paramShapes[k]) 
     232            thisScale = 0.1 * (k+1) 
     233            self.failUnlessEqual(priorArray.scale.max(), thisScale) 
     234            self.failUnlessEqual(priorArray.scale.min(), thisScale) 
     235 
     236    def _check_basic(self, paramContainer): 
    235237        self.failUnless(isinstance(paramContainer, params.ParameterContainer)) 
    236238        self.failUnless(isinstance(paramContainer.Lp, float)) 
    237239        for k, name in enumerate(self.priorContainer.paramNames): 
    238240            paramArray = getattr(paramContainer, name) 
    239             self.failUnlessEqual(paramArray.shape, self.shapes[k]) 
    240  
    241     def test_proposal(self): 
    242         self.fail("TODO") 
     241            self.failUnlessEqual(paramArray.shape, self.MPC.paramShapes[k]) 
     242 
     243    def _check_dist(self, paramArray, scales=None, means=False): 
     244        for k, name in enumerate(self.priorContainer.paramNames): 
     245            if scales is None: 
     246                thisScale = 0.1 * (k+1) 
     247            else: 
     248                thisScale = scales[k] 
     249            arrays = getattr(paramArray, name)  
     250            if means: 
     251                x = arrays.call('mean') 
     252            else: 
     253                x = arrays.call('item', 0) 
     254            self.failUnlessNormal(x) 
     255            self.failUnlessAlmostEqual(x.mean()/thisScale, 0.0, 0) 
     256            stdExpected = x.std()/thisScale 
     257            if means: 
     258                stdExpected /= s.sqrt(len(arrays[0].ravel())) 
     259            self.failUnlessAlmostEqual(stdExpected, 1.0, 0) 
     260         
     261    def test_new_basic(self): 
     262        return self._check_basic(self.priorContainer.new()) 
     263 
     264    def test_new_dist(self): 
     265        N = 200 
     266        paramArray = params.FlexArray(N) 
     267        for k in xrange(N): 
     268            paramArray[k] = self.priorContainer.new() 
     269        return self._check_dist(paramArray) 
     270 
     271    def test_proposal_basic(self): 
     272        old = util.Mock_ParameterContainer(0.0) 
     273        new = self.priorContainer.proposal(old, 1.0) 
     274        return self._check_basic(new) 
    243275     
     276    def test_proposal_dist(self): 
     277        N = 5000 
     278        sigma = 0.12 
     279        paramArray = params.FlexArray(N) 
     280        print "\nRunning MCMC: " 
     281        X = util.Mock_ParameterContainer(0.0) 
     282        #X = self.priorContainer.new() 
     283        for k in xrange(N): 
     284            XP = self.priorContainer.proposal(X, sigma) 
     285            if k == 0 or XP.Lp > X.Lp or s.log(s.rand()) < (XP.Lp - X.Lp): 
     286                print "+", 
     287                X = XP 
     288            else: 
     289                print ".", 
     290            sys.stdout.flush() 
     291            paramArray[k] = X 
     292        print "\n" 
     293        scale = 0.1 * sigma * s.arange( 
     294            1, len(self.priorContainer.paramNames)+1) 
     295        return self._check_dist(paramArray, scale) 

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

    r159 r168  
    1 # PyFolio: Portfolio Performance Simulation 
     1# svpmc: Stochastic Volatility Inference via Population Monte Carlo 
    22# 
    3 # Copyright (C) 2007 by Edwin A. Suominen, http://www.eepatents.com 
     3# Copyright (C) 2007-2008 by Edwin A. Suominen, http://www.eepatents.com 
    44# 
    55# This program is free software; you can redistribute it and/or modify it under 
    …… 
    1717 
    1818""" 
    19 Unit tests for pyfolio.mcmc 
     19Unit tests for svpmc.pmc 
    2020""" 
    2121 
    …… 
    2727from asynqueue import ThreadQueue 
    2828 
    29 import model, mcmc 
     29import model, pmc 
    3030from sample import resample 
    3131import util 
    3232 
    3333 
    34 VERBOSE = False 
     34VERBOSE = util.Mock.verbose() 
    3535 
    3636 
    …… 
    4141 
    4242 
    43 class Mock_ModelManager(util.Mock): 
     43class Mock_ProjectManager(util.Mock): 
    4444    def __init__(self, modelObj): 
    4545        self.modelObj = modelObj 
    …… 
    6969 
    7070class BaseTC(util.TestCase): 
    71     class IterationConsumer(object): 
    72         implements(interfaces.IConsumer) 
    73         def __init__(self, mgr, modulus): 
    74             self.mgr, self.modulus = mgr, modulus 
    75         def registerProducer(self, producer, streaming): 
    76             self.X, self.stats = [], [] 
    77             print "\nSample means, standard deviations, and sequential steps" 
    78         def unregisterProducer(self): 
    79             pass 
    80         def write(self, X): 
    81             def a2s(V): 
    82                 joined = ", ".join(["%4.3f" % x for x in V]) 
    83                 return "[%s]" % joined 
    84  
    85             count = len(self.X) 
    86             if count % self.modulus == 0: 
    87                 mean, std = X.mean(0), X.std(0) 
    88                 self.mgr.caller.resetCalls() 
    89                 if VERBOSE: 
    90                     print "%4d: %15s, %15s" % (count, a2s(mean), a2s(std)) 
    91                 self.stats.append((mean, std)) 
    92             self.X.append(X.copy()) 
    93  
    94     def _get_consumer(self): 
    95         if getattr(self, '_consumer', None) is None: 
    96             self._consumer = self.IterationConsumer(self.mgr, 5) 
    97         return self._consumer 
    98     def _set_consumer(self, value): 
    99         self._consumer = value 
    100     consumer = property(_get_consumer, _set_consumer) 
    101      
    10271    def setUp(self): 
    10372        self.mgr = model.ModelManager(projectManager) 

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

    r167 r168  
    115115        self.cdfPath = os.path.join(projectDir, "svpmc.nc") 
    116116        self.specFile = os.path.join(projectDir, "project-spec.txt") 
    117         self.mgr = project.ProjectManager(self.specFile, m=3) 
     117        self.mgr = project.ProjectManager(self.specFile, m=5) 
    118118 
    119119    def _makeParameters(self, *scales): 
    …… 
    125125    def test_writeParams(self): 
    126126        # Write them 
    127         parameters = self._makeParameters(7, 11, 13) 
     127        parameters = self._makeParameters(3, 7, 11, 13, 17) 
    128128        self.mgr.writeParams(parameters) 
    129129        self.mgr.cdf.close() 

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

    r167 r168  
    2525import os, copy, time, imp, tarfile 
    2626 
    27 import scipy as s  
     27import scipy as s 
     28from scipy import stats 
    2829from twisted.internet import defer, reactor 
    2930from twisted.spread import pb 
    …… 
    4142    x, os.path.join(os.path.dirname(__file__), "%s-us.dat" % x)).v2r() 
    4243    for x in ('dm', 'jy')] 
     44 
     45PARAM_NAMES = ['a', 'b', 'cr', 'd', 'e', 'fs', 'fr', 'g'] 
    4346 
    4447 
    …… 
    163166 
    164167class Mock_ParameterContainer(Mock): 
     168    N = 937 
    165169    p, xcorrs = 3, 6 
    166     paramNames =  ['a', 'b', 'cr', 'd', 'e', 'fs', 'fr', 'g'] 
     170    paramNames =  PARAM_NAMES 
    167171    paramShapes = [(p,), (p,p), (xcorrs,), (p,), (p,p), (p,), (xcorrs,), (p,)] 
     172 
     173    Lp = -s.inf 
     174    z = s.zeros(N) 
    168175 
    169176    def __init__(self, scale=1.0): 
    …… 
    173180        self.h = self.p * scale * s.ones(self.p) 
    174181        self.Lx = 0.0 
     182 
     183    def paramSequence(self): 
     184        result = [] 
     185        for k, name in enumerate(self.paramNames): 
     186            value = getattr(self, name, None) 
     187            if value is None: 
     188                value = s.zeros(paramShapes[k]) 
     189            result.append(value) 
     190        return result 
    175191 
    176192 
    …… 
    251267            error = sum(error) 
    252268        self.failUnlessAlmostEqual(error, 0.0, places) 
     269 
     270    def failUnlessNormal(self, x): 
     271        p = stats.normaltest(x)[1] 
     272        if Mock.verbose(): 
     273            fig = self.fig 
     274            sp = fig.add_subplot(211) 
     275            sp.plot(x) 
     276            sp = fig.add_subplot(212) 
     277            sp.hist(x, bins=20) 
     278        self.failUnless( 
     279            p > 0.05, "Deviation from normality with significance p=%f" % p) 
    253280 
    254281    @defer.deferredGenerator