Changeset 132

Show
Ignore:
Timestamp:
03/31/08 19:21:24 (8 months ago)
Author:
edsuom
Message:

Working on refactoring pypmc stuff in from separate project repo; tseries and sample testing OK

Files:

Legend:

Unmodified
Added
Removed
Modified
Copied
Moved

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

    r131 r132  
    2626from twisted.internet import defer, reactor 
    2727 
    28 #from asynqueue import ThreadQueue 
    2928from asyncluster.master import jobs 
    3029from twisted_goodies.pybywire.params import Parameterized 
    3130from twisted_goodies.pybywire.pack import Packer, Unpacker 
    3231 
    33 from dist import Distribution 
    34  
    35  
    36 class ThreadQueue(object): 
     32 
     33class NullQueue(object): 
     34    """ 
     35    I act like an AsynQueue except that I run everything in the main 
     36    thread. I'm mostly useful for debugging. 
     37    """ 
    3738    def __init__(self, *args): 
    3839        pass 
    …… 
    4546 
    4647 
    47 class ModelManager(object): 
    48     """ 
    49     I manage an ensemble of L{Model} instances, with a parameter vector that is 
    50     the joint set of all of the models' parameters. 
    51  
    52     I am constructed with a sequence of parameter names, an equal-length 
    53     sequence of L{Prior} objects for the parameters, and a sequence I{models} 
    54     of 2-tuples for the models. 
    55  
    56     Each tuple in the models sequence contains: 
    57  
    58       1. A reference to one of the models, and 
    59  
    60       2. A sequence of expressions for the normalization and distribution 
    61          parameters as supplied to that model in each parameter vector. 
    62      
    63     """ 
    64     def __init__(self, projectManager): 
    65         ID = id(self) 
    66         self.mgr = projectManager 
    67         for name in ('paramNames', 'priors', 'expressions'): 
    68             setattr(self, name, getattr(projectManager, name)) 
    69         self.N_params = len(self.paramNames) 
    70         self.modelMap, self.peMap = {}, {} 
    71         for modelObj, paramExprs in projectManager.models: 
    72             if modelObj.N_params() != len(paramExprs): 
    73                 exprString = ", ".join(["'%s'" % x for x in paramExprs]) 
    74                 msg = "Supplied expressions (%s) " % exprString +\ 
    75                       "do not match %d model parameters" % modelObj.N_params() 
    76                 raise ValueError(msg) 
    77             self.modelMap[ID] = modelObj 
    78             self.peMap[ID] = paramExprs 
    79             ID += 1 
    80         self.caller = ModelCaller(self.modelMap) 
    81  
    82     def getID(self, modelObj): 
    83         """ 
    84         Returns the integer ID for the supplied I{modelObj}. 
    85         """ 
    86         for ID, thisModelObj in self.modelMap.iteritems(): 
    87             if thisModelObj == modelObj: 
    88                 return ID 
    89      
    90     def setLocalMode(self): 
    91         """ 
    92         See L{ModelCaller.setLocalMode}. 
    93         """ 
    94         return self.caller.setLocalMode() 
    95  
    96     def setRemoteMode(self, socket): 
    97         """ 
    98         See L{ModelCaller.setRemoteMode}. 
    99         """ 
    100         return self.caller.setRemoteMode(socket) 
    101  
    102     def modelParams(self, X): 
    103         """ 
    104         Based on the parameter vectors supplied in the 2d array I{X}, or a 
    105         single parameter vector supplied as I{X}, returns a dict of parameter 
    106         vectors evaluated for each of my models from their parameter 
    107         expressions in I{peMap}. The vectors are keyed by the ID of the model 
    108         to which they apply. 
    109         """ 
    110         if len(X.shape) == 1: 
    111             X = X.reshape(1, X.shape[0]) 
    112         XM = {} 
    113         namespace = {'U':s.ones(len(X))} 
    114         for k, paramName in enumerate(self.paramNames): 
    115             namespace[paramName] = X[:,k] 
    116         for name, expression in self.expressions.iteritems(): 
    117             namespace[name] = eval(expression, namespace) 
    118         for ID, paramExprs in self.peMap.iteritems(): 
    119             evalCode = "[%s]" % ",".join(paramExprs) 
    120             XM[ID] = s.column_stack(eval(evalCode, namespace)) 
    121         return XM 
    122  
    123     def rPriors(self, N): 
    124         """ 
    125         Returns I{N} parameter vectors as random draws from the prior 
    126         distributions of the parameters. 
    127         """ 
    128         Y = s.empty((N, self.N_params), s.float32) 
    129         for k, p in enumerate(self.priors): 
    130             Y[:,k] = p.rvs(N) 
    131         return Y 
    132  
    133     def pPriors(self, X): 
    134         """ 
    135         Returns a vector of prior probability densities for each parameter in 
    136         X, or rows of X. 
    137         """ 
    138         if len(X.shape) == 1: 
    139             X = X.reshape(1, X.shape[0]) 
    140         N = X.shape[0] 
    141         Y = s.zeros((N, self.N_params), s.float32) 
    142         for k, p in enumerate(self.priors): 
    143             Y[:,k] = p.pdf(X[:,k]) 
    144         return Y 
    145  
    146     def rProposal(self, N, wiggle): 
    147         """ 
    148         Returns an array of I{N} rows of parameter offsets drawn from the prior 
    149         distributions scaled by a I{wiggle} value between 0 and 1. 
    150         """ 
    151         Y = s.empty((N, self.N_params), s.float32) 
    152         for k, p in enumerate(self.priors): 
    153             Y[:,k] = p.rds(N, wiggle) 
    154         return Y 
    155  
    156     def pProposal(self, X, wiggle): 
    157         """ 
    158         Returns a vector of probability densities for each parameter offset in 
    159         X, or rows of X, under the assumption that they were generated from my 
    160         L{rProposal} method with the specified I{wiggle}. 
    161         """ 
    162         if len(X.shape) == 1: 
    163             X = X.reshape(1, X.shape[0]) 
    164         N = X.shape[0] 
    165         Y = s.zeros((N, self.N_params), s.float32) 
    166         for k, p in enumerate(self.priors): 
    167             Y[:,k] = p.pds(X[:,k], wiggle) 
    168         return Y 
    169  
    170     @defer.deferredGenerator 
    171     def rDistribution(self, X, N): 
    172         """ 
    173         Given the parameter vector I{X}, returns a dict of time series vectors 
    174         of I{N} random variates from the probability distribution for each 
    175         model, keyed by model ID. 
    176         """ 
    177         results = {} 
    178         Distribution.rvCache.clear() 
    179         XM = self.modelParams(X) 
    180         for ID, Xm in XM.iteritems(): 
    181             wfd = defer.waitForDeferred( 
    182                 self.caller.call(ID, 'setParamVector', Xm[0])) 
    183             yield wfd; wfd.getResult() 
    184             wfd = defer.waitForDeferred( 
    185                 self.caller.call(ID, 'rDistribution', N)) 
    186             yield wfd 
    187             results[ID] = wfd.getResult() 
    188         yield results 
    189  
    190     def pDistribution(self, X, N, M=1): 
    191         """ 
    192         Given the parameter vector I{X}, returns a dict of 2-tuples that each 
    193         contain a 1d array of I{N} equispaced points over the range of the 
    194         normalized data for each model, and a 1d array of probability densities 
    195         for the data points. 
    196  
    197         The returned vectors are keyed by the ID of the model to which they 
    198         apply. 
    199  
    200         You can compute the probability densities over a multiple of each 
    201         model's data ranges by setting the keyword I{M} to something greater 
    202         than 1. 
    203         """ 
    204         def called(result, ID): 
    205             results[ID] = result 
    206          
    207         dList, results = [], {} 
    208         XM = self.modelParams(X) 
    209         for ID, Xm in XM.iteritems(): 
    210             d = self.caller.call(ID, 'pDistribution', Xm[0], N, M) 
    211             d.addCallback(called, ID) 
    212             dList.append(d) 
    213         return defer.DeferredList(dList).addCallback(lambda _: results) 
    214      
    215     def likelihoods(self, X): 
    216         """ 
    217         Returns the total log likelihoods of the global parameter vector(s) in 
    218         the supplied array I{X}. 
    219  
    220         Runs via thread or remotely, with a deferred result either way. Tries 
    221         twice for each parameter vector, returning infinitely negative 
    222         likelihood for repeated failures. 
    223         """ 
    224         def called(result, paramVectors, k): 
    225             if not isinstance(result, float): 
    226                 d = self.caller.likelihood(IDs, paramVectors) 
    227                 d.addBoth(triedAgain, k) 
    228                 return d 
    229             L[k] = result 
    230              
    231         def triedAgain(result, k): 
    232             if not isinstance(result, float): 
    233                 result = -s.inf 
    234             L[k] = result 
    235  
    236         XM = self.modelParams(X) 
    237         IDs = XM.keys() 
    238         dList = [] 
    239         N = XM.values()[0].shape[0] 
    240         L = s.empty(N) 
    241         for k in xrange(N): 
    242             paramVectors = [XM[ID][k] for ID in IDs] 
    243             d = self.caller.likelihood(IDs, paramVectors) 
    244             d.addBoth(called, paramVectors, k) 
    245             dList.append(d) 
    246         return defer.DeferredList(dList).addCallback(lambda _: L) 
    247  
    248  
    24948class ModelCaller(object): 
    25049    """ 
    251     I handle deferred calls to models in an ensemble, in either local or remote 
    252     mode. Calls in remote mode are run through a separate thread. 
    253  
    254     Instantiate me with a dict of the models, keyed by unique integer IDs. 
     50    I handle deferred calls to a L{Model} object in either local or remote 
     51    mode. Calls in local mode are run through an instance of L{NullQueue} in 
     52    the main thread. 
    25553    """ 
    25654    remoteMode = False 
    …… 
    26260    from twisted_goodies.pybywire import pack 
    26361     
    264     G = {} 
    265  
    266     def newModel(ID, model): 
    267         G[ID] = model 
    268  
    269     def callModel(ID, methodName, *possiblyPackedArgs): 
    270         method = getattr(G[ID], methodName) 
    271         return packwrap(method)(*possiblyPackedArgs) 
    272  
    273     def likelihood(IDs, paramVectors): 
    274         L = 0 
    275         for k, X in enumerate(pack.Unpacker(paramVectors)): 
    276             ID = IDs[k] 
    277             L += G[ID].likelihood(X) 
    278             if not isfinite(L): 
    279                 break 
     62    M = None 
     63 
     64    def newModel(modelObj): 
     65        global M 
     66        M = modelObj 
     67 
     68    def callModel(methodName, *possiblyPackedArgs): 
     69        method = getattr(M, methodName) 
     70        return pack.packwrap(method)(*possiblyPackedArgs) 
     71 
     72    def likelihood(packedVector): 
     73        paramVector = pack.Unpacker(packedVector)() 
     74        L = M.likelihood(paramVector) 
    28075        return str(float(L)) 
    28176 
    28277    ########################################################################### 
    28378    """ 
    284      
    285     def __init__(self, modelMap): 
    286         self.modelMap = modelMap 
    287         self.threadQueue = ThreadQueue(1) 
     79    def __init__(self, modelObj): 
     80        self.modelObj = modelObj 
     81        self.queue = NullQueue(1) 
    28882        reactor.addSystemEventTrigger( 
    289             'before', 'shutdown',  self.threadQueue.shutdown) 
     83            'before', 'shutdown',  self.queue.shutdown) 
    29084 
    29185    def _oops(self, failure): 
    …… 
    30498        remotely-run operations can commence. 
    30599        """ 
    306         def dispatchModels(null): 
    307             dList = [ 
    308                 self.client.update('newModel', ID, modelObj) 
    309                 for ID, modelObj in self.modelMap.iteritems()] 
    310             return defer.DeferredList(dList) 
    311  
    312100        if hasattr(self, 'client'): 
    313101            return defer.succeed(None) 
    …… 
    318106        d.addCallback( 
    319107            lambda _: self.client.registerClasses(*Parameterized.registry)) 
    320         d.addCallback(dispatchModels) 
     108        d.addCallback( 
     109            lambda _: self.client.update('newModel', modelObj)) 
    321110        d.addCallback(lambda _: setattr(self, 'remoteMode', True)) 
    322111        d.addErrback(self._oops) 
    …… 
    330119        return result 
    331120 
    332     def call(self, ID, methodName, *args, **kw): 
    333         """ 
    334         For the model specified by integer I{ID}, calls the method identified 
    335         by I{methodName} with any supplied args. 
     121    def call(self, methodName, *args, **kw): 
     122        """ 
     123        Calls my model's method identified by I{methodName} with any supplied 
     124        args. 
    336125 
    337126        If I am in local mode, runs the call in a dedicated thread of the 
    …… 
    352141                args = (Packer(*args)(),) 
    353142            d = self.client.run( 
    354                 'callModel', ID, methodName, *args, **{'timeout':self.timeout}) 
     143                'callModel', methodName, *args, **{'timeout':self.timeout}) 
    355144            d.addCallback(self._gotPossiblyPackedResult) 
    356145        else: 
    357             d = self.threadQueue.call( 
    358                 getattr(self.modelMap[ID], methodName), *args) 
     146            d = self.queue.call(getattr(self.modelObj, methodName), *args) 
    359147        d.addErrback(self._oops) 
    360148        return d 
    361149 
    362     def likelihood(self, IDs, paramVectors, remote=False, local=False): 
    363         """ 
    364         Returns the total log likelihood of the parameter vectors for the 
    365         models referenced by ID. The first argument is a sequence of the model 
    366         I{IDs} and the second is a sequence of the models' parameter vectors, 
    367         in order. 
     150    def likelihood(self, paramVector, remote=False, local=False): 
     151        """ 
     152        Returns the log likelihood of the supplied parameter vector for my 
     153        model. 
    368154        """ 
    369155        def gotLikelihood(value): 
    …… 
    375161         
    376162        if (remote or self.remoteMode) and not local: 
    377             packedVectors = Packer(*paramVectors) 
     163            packedVector = Packer(*paramVector) 
    378164            d = self.client.run( 
    379                 'likelihood', IDs, packedVectors(), timeout=self.timeout) 
     165                'likelihood', packedVector(), timeout=self.timeout) 
    380166            d.addCallback(gotLikelihood) 
    381167        else: 
    382             dList = [self.threadQueue.call( 
    383                 self.modelMap[ID].likelihood, paramVectors[k]) 
    384                 for k, ID in enumerate(IDs)] 
    385             d = defer.gatherResults(dList).addCallback(lambda x: sum(x)) 
     168            d = self.queue.call(self.modelObj.likelihood, paramVector) 
    386169        return d.addErrback(self._oops) 
    387170 
    …… 
    451234        """ 
    452235        Returns the probability densities of the supplied zero-mean deviates, 
    453         under the assumption that they were generated from my L{rds} method 
    454         with the specified I{wiggle}. 
    455         """ 
    456         # Taking out divide-by-wiggle helps avoid low-likelihood tagalongs 
    457         # but makes for very lumpy posteriors. Which is the correct way? 
     236        under the assumption that they were generated from my Gaussian-PDF 
     237        L{rds} method with the specified I{wiggle}. 
     238 
     239        For my Gaussian-PDF deviate generator, the probability density of a 
     240        given deviate is inversely proportional to the scaling to its standard 
     241        deviation imparted by I{wiggle}. Thus, the unit-normal PDF is divided 
     242        by the I{wiggle} in the result. 
     243        """ 
    458244        return self.rdsObj.pdf(X/wiggle) / wiggle 
    459245 
    …… 
    461247class Model(Parameterized): 
    462248    """ 
    463     Stochastic Model for fitting a random variable, or a linear combination 
    464     thereof, to a set of observations. 
    465  
    466     The probability distributions of the random variable(s) is embodied in the 
    467     my I{distribution} attribute, an instance of L{dist.Distribution} or, for a 
     249    Stochastic Model for fitting a C{Kx1} vector of random variables to a 
     250    C{KxN} matrix of observations. 
     251 
     252    The probability distributions of the random variables is embodied in the my 
     253    I{distribution} attribute, an instance of L{dist.Distribution} or, for a 
    468254    combination of random variables, an instance of L{dist.Combo_Distribution}. 
    469255 

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

    r131 r132  
    1 # PyFolio: Portfolio Performance Simulation 
     1# PyPMC: Population Monte Carlo implemented with Python and Twisted 
    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 
    …… 
    3030 
    3131import asynqueue 
     32from twisted_goodies.pybywire.pack import Packer, Unpacker 
    3233from asyncluster.master import jobs 
    33 from twisted_goodies.pybywire.pack import Packer, Unpacker 
    3434 
    3535import model 
    …… 
    9191 
    9292 
    93 class MC_Base(object): 
     93class PMC(object): 
    9494    """ 
    95     I provide a common foundation for Monte Carlo simulators. 
     95    Population MCMC with per Cappe et al. 
    9696    """ 
     97    wiggle = 1.0 
     98    chunkSize = 5000 
     99    V = [1.0, 0.5, 0.1, 0.05, 0.001] 
     100 
    97101    def __init__(self, modelManager, **kw): 
    98102        self.mgr = modelManager 
    …… 
    131135        """ 
    132136        self.iterationProducer.registerConsumer(consumer) 
    133  
    134  
    135 class MCMC(MC_Base): 
    136     """ 
    137     I am a conventional Metropolis-Hastings MCMC sampler, and serve as a base 
    138     class for the DE-MCMC sampler. 
    139     """ 
    140     wiggle = 0.05 
    141     reportInterval = 2 
    142  
    143     def getXL(self, N): 
    144         """ 
    145         Generates a (N_pop, N_params+1) population array I{XL}. Each row 
    146         contains the one set of distribution parameters concatenated with the 
    147         scalar value of the log-likelihood of those parameters, given the data. 
    148  
    149         Returns a deferred that fires with the array. 
    150         """ 
    151         def gotLikelihood(L, m, X): 
    152             XL[m,:-1] = X 
    153             XL[m,-1] = L[0] 
    154  
    155         def gotAllLikelihoods(null): 
    156             XL[:,-1] += s.log(self.mgr.pPriors(XL[:,:-1])).sum(1) 
    157             return XL 
    158  
    159         # TODO: Use array processing in model manager's likelihoods method 
    160         XL = s.empty((N, self.N_params+1), s.float32) 
    161         X = self.mgr.rPriors(N) 
    162         dList = [] 
    163         for m, Xm in enumerate(X): 
    164             d = self.mgr.likelihoods(Xm) 
    165             d.addCallback(gotLikelihood, m, Xm) 
    166             dList.append(d) 
    167         d = defer.DeferredList(dList) 
    168         d.addCallback(gotAllLikelihoods) 
    169         return d 
    170  
    171     def E(self, wiggle): 
    172         """ 
    173         Returns a random parameter offset vector drawn from the model's 
    174         proposal distribution. 
    175  
    176         Draws the vectors from the model in batches for improved efficiency. 
    177         """ 
    178         if not hasattr(self, '_OV') or self._kOV == 499: 
    179             self._kOV = -1 
    180             self._OV = self.mgr.rProposal(500, wiggle) 
    181         self._kOV += 1 
    182         return self._OV[self._kOV,:] 
    183  
    184     @defer.deferredGenerator 
    185     def run(self, N_chains, **kw): 
    186         """ 
    187         Does an MCMC or DE-MCMC run with I{N_chains} parallel chains or 
    188         population members. 
    189  
    190         Returns a deferred that fires when the run is done. No output value is 
    191         provided via the deferred, however. 
    192  
    193         Use a provider of L{interfaces.IConsumer} to get the MCMC samples as 
    194         they are produced. Each attached consumer is updated after every 
    195         iteration with my array of parameter vectors (one vector per row) for 
    196         each iteration. 
    197  
    198         @param N_iter: The number of iterations to produce after burn-in. 
    199  
    200         @param N_bi: The number of burn-in iterations to discard before 
    201           producing any. 
    202          
    203         """ 
    204         N_iter, N_bi = self.setupRun(N_chains, **kw) 
    205         wfd = defer.waitForDeferred(self.getXL(N_chains)) 
    206         yield wfd 
    207         self.XL = wfd.getResult() 
    208         i = 0 
    209         acceptances = [] 
    210         while i < N_iter + N_bi: 
    211             d = defer.gatherResults(self.iteration()) 
    212             t0 = time.time() 
    213             wfd = defer.waitForDeferred(d) 
    214             yield wfd; 
    215             acceptances.extend(wfd.getResult()) 
    216             if i >= N_bi: 
    217                 self.iterationProducer.gotNext(self.XL[:,:-1]) 
    218             if not i % self.reportInterval: 
    219                 acceptanceRate = 100.0 * sum(acceptances) / len(acceptances) 
    220                 acceptances = [] 
    221                 msg = "%4d, %4.2f sec, AR=%3d%%, LPmax=%d: %s" % ( 
    222                     i, time.time()-t0, acceptanceRate, self.XL[:,-1].max(), 
    223                     ", ".join(["%8.3g" % x for x in self.XL[:,:-1].mean(0)])) 
    224                 print msg 
    225             i += 1 
    226         self.iterationProducer.stopProducing() 
    227  
    228     def iteration(self): 
    229         """ 
    230         Does a conventional MCMC iteration. 
    231         """ 
    232         return [ 
    233             self.metropolisHastings(m) 
    234             for m in xrange(self.N_chains)] 
    235  
    236     def metropolisHastings(self, m, Xd=None): 
    237         """ 
    238         Proposes a new parameter vector for population member or chain I{m} 
    239         that is offset from the last one by a draw from the model's proposal 
    240         distribution plus any offset vector I{Xd} that is supplied. 
    241  
    242         The parameter is accepted with probability proportional to the 
    243         Metropolis-Hastings ratio of its likelihood versus the likelihood of 
    244         the existing vector. 
    245  
    246         If accepted, the proposal vector overwrites the existing one in my XL 
    247         matrix. 
    248  
    249         Returns a deferred that fires with C{True} if the proposal was 
    250         accepted, C{False} if not. 
    251         """ 
    252         def gotL(Lp): 
    253             Lp = Lp[0] 
    254             # The difference between log probabilities is the same as the ratio 
    255             # of the probabilities themselves. 
    256             if s.isfinite(Lp): 
    257                 Lp += s.log(P).sum() 
    258                 if (Lp - self.XL[m,-1]) > s.log(s.rand()): 
    259                     # The proposal was accepted, so write the new parameter 
    260                     # vector and its log likelihood to the XL matrix. 
    261                     self.XL[m,:-1] = Xp 
    262                     self.XL[m,-1] = Lp 
    263                     return True 
    264             # The proposal was rejected, so I'll leave the current copy alone. 
    265             return False 
    266  
    267         if Xd is None: 
    268             # Standard MCMC, where the wiggle scales the jump 
    269             Xp = self.XL[m,:-1] + self.E(self.wiggle) 
    270         else: 
    271             # DE-MCMC, where the wiggle scales the offset vector (=gamma) 
    272             Xp = self.XL[m,:-1] + self.wiggle*Xd + self.E(0.0005) 
    273         P = self.mgr.pPriors(Xp) 
    274         if s.greater(P, 0).all(): 
    275             d = self.mgr.likelihoods(Xp) 
    276             # A failure in the likelihood call is treated like a rejected 
    277             # proposal 
    278             d.addCallbacks(gotL, lambda _: False) 
    279         else: 
    280             # A zero-probability prior guarantees a rejection, so don't bother 
    281             # computing the likelihood 
    282             d = defer.succeed(False) 
    283         return d 
    284  
    285  
    286 class DE_MCMC(MCMC): 
    287     """ 
    288     Population MCMC with an asynchronous, parallelized version of the 
    289     Differential Evolution solution of Cajo J. F. Ter Braak, Stat Comput (2006) 
    290     16:239-249.  The samples of each chain at any given iteration are members 
    291     of a population at one evolutionary generation. 
    292  
    293     Each proposal is based on a scaled difference of two parameter vectors from 
    294     the population, other than the vector being modified, and a random variate 
    295     vector with the same length as the parameter vector. The random variates 
    296     are scaled and then added to the scaled difference vector, in accordance 
    297     with Ter Braak (2). 
    298     """ 
    299     reportInterval = 1 
    300      
    301     def _get_wiggle(self): 
    302         """ 
    303         Implements gamma=1 every tenth iteration, per Ter Braak, p. 245. 
    304         """ 
    305         if hasattr(self, '_wiggle'): 
    306             return self._wiggle 
    307         self._gammaCallCount = getattr(self, '_gammaCallCount', 0) + 1 
    308         if self._gammaCallCount % 10: 
    309             return 2.38 / s.sqrt(2*self.N_params) 
    310         return 1.0 
    311     def _set_wiggle(self, value): 
    312         self._wiggle = value 
    313     wiggle = property(_get_wiggle, _set_wiggle) 
    314  
    315     def getXds(self): 
    316         """ 
    317         Returns an array of offset vectors, one for the evolution of each 
    318         population member in the next iteration, with no scaling. 
    319         """ 
    320         Xds = s.empty((self.XL.shape[0], self.XL.shape[1]-1), s.float32) 
    321         for m0 in xrange(self.N_chains): 
    322             m12 = [] 
    323             while len(m12) < 2: 
    324                 r = random.randint(0, self.N_chains) 
    325                 if r not in m12: 
    326                     m12.append(r) 
    327             Xds[m0,:] = self.XL[m12[1],:-1] - self.XL[m12[0],:-1] 
    328         return Xds 
    329  
    330     def iteration(self): 
    331         """ 
    332         Does a DE-MCMC iteration. 
    333         """ 
    334         Xds = self.getXds() 
    335         return [ 
    336             self.metropolisHastings(m0, Xd) 
    337             for m0, Xd in enumerate(Xds)] 
    338          
    339  
    340 class PMC(MC_Base): 
    341     """ 
    342     Population MCMC with per Cappe et al. 
    343     """ 
    344     wiggle = 1.0 
    345     chunkSize = 5000 
    346     V = [1.0, 0.5, 0.1, 0.05, 0.001] 
    347137 
    348138    def _get_queue(self): 

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

    r131 r132  
    1 # PyFolio: Portfolio Performance Simulation 
     1# PyPMC: Population Monte Carlo implemented with Python and Twisted 
    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 
    …… 
    8383    I1 = s.greater(V, x[RI,0]).choose(RI, x[RI,1].astype(int)) 
    8484    return I0[I1] 
    85  
    86  
    87 class TrapSampler(object): 
    88     """ 
    89     """ 
    90     def __init__(self, X, Y): 
    91         self.X, self.Y = X, Y 
    92         self.W = s.diff(X) * (Y[:-1] + Y[1:]) 
    93         self.W /= sum(self.W) 
    94  
    95     def trapSample(self, x0, y0, x1, y1): 
    96         """ 
    97         The probability density function (pdf) for variates within a selected 
    98         trapezoidal portion of the piece-wise approximated distribution curve 
    99         is:: 
    100  
    101                 (          (y1 - y0)      ) 
    102           y = C ( (x - x0) --------- + y0 ) 
    103                 (          (x1 -x0)       ) 
    104  
    105         where C is a normalizing constant. 
    106  
    107         The cumulative density function (cdf) for the trapezoidal pdf is 
    108         based on the definite integral of the pdf:: 
    109          
    110                            2         2 
    111                 (2 x x0 - x ) y1 + (x  - 2 x x1) y0 
    112           y = - ----------------------------------- 
    113                        2                   2 
    114                     (x1  - 2 x0 x1) y1 + x1  y0 
    115  
    116         Solving for x yields the inverse cdf. With inverse transform sampling, 
    117         y is set to a uniform random variate U(0,1). 
    118         """ 
    119         N = len(x0) 
    120         y = s.rand(N) 
    121         x = x0 + (x1-x0)*y 
    122         I = s.not_equal(y0, y1).nonzero() 
    123         xI  = -y0[I] + s.sqrt(y0[I]**2 + y[I]*(y1[I]**2 - y0[I]**2)) 
    124         xI *= (x1[I] - x0[I])/(y1[I] - y0[I]) 
    125         x[I] = xI + x0[I] 
    126         return x 
    127      
    128     def __call__(self, N): 
    129         """ 
    130         Call me to get I{N} random variates. 
    131         """ 
    132         I = resample(self.W, N) 
    133         return self.trapSample(self.X[I], self.Y[I], self.X[I+1], self.Y[I+1]) 
    134