Changeset 151

Show
Ignore:
Timestamp:
04/19/08 20:26:56 (8 months ago)
Author:
edsuom
Message:

Tested FlexArray?; last work on Model.likelihood before switching it to Hybrid Gibbs approach

Files:

Legend:

Unmodified
Added
Removed
Modified
Copied
Moved

  • projects/AsynCluster/tags/20080419/svpmc/model.py

    r150 r151  
    7070 
    7171    def likelihood(paramContainer, wiggle): 
    72         v, L = M.likelihood(paramContainer, wiggle) 
    73         return Packer(v, L)() 
     72        L, v = M.likelihood(paramContainer, wiggle) 
     73        return Packer(L, v)() 
    7474 
    7575    ########################################################################### 
    …… 
    112112    def likelihood(self, paramContainer, wiggle, remote=False, local=False): 
    113113        """ 
    114         Returns the log-likelihood of the parameters in the supplied 
     114        Computes the log-likelihood of the parameters in the supplied 
    115115        I{paramContainer} for my model. 
    116116 
    117         Updates in-place the container's sole latent parameter, an array of 
    118         volatility shocks after some MCMC involved with the likelihood 
    119         compution. 
     117        Updates the container in-place with the log-likelihood and an array of 
     118        volatility shocks on which the likelihood computation is based, after 
     119        undergoing some nested MCMC. 
     120 
     121        Returns a deferred that fires (with C{None}) when the likelihood 
     122        computation is done and the paramContainer has been updated. 
    120123        """ 
    121124        def gotPackedLikelihood(result): 
    …… 
    123126                # An error from the remote likelihood method call is treated as 
    124127                # infinitely low likelihood 
    125                 return -s.inf 
     128                paramContainer.pLogLikelihood = -s.inf 
    126129            up = pack.Unpacker(result) 
    127             paramContainer.v = up() 
    128             return up() 
     130            paramContainer.pLogLikelihood = up() 
     131            paramContainer.latent = up() 
    129132 
    130133        def gotLikelihood(result): 
    131             paramContainer.v = result[0] 
    132             return result[1] 
     134            paramContainer.pLogLikelihood = result[0] 
     135            paramContainer.latent = result[1] 
    133136         
    134137        if (remote or self.remoteMode) and not local: 
    …… 
    210213    def _get_var(self): 
    211214        if 'var' not in self.cache: 
    212             self.cache['var'] = VAR(self.tsList) 
     215            self.cache['var'] = VAR(self.y) 
    213216        return self.cache['var'] 
    214217    var = property(_get_var) 
    …… 
    254257            if hs != xs: 
    255258                raise ValueError("Need one log-volatility per innovation") 
    256             x *= s.exp(-0.5*h) 
     259            x = s.exp(-0.5*h) * x.copy() 
    257260        if 'decorrelate' not in self.cache: 
    258261            self.cache['decorrelate'] = linalg.inv( 
    …… 
    268271        The result will be the total log-likelihood of the innovations, given 
    269272        my volatility/innovation shock correlations C{g} and the volatility 
    270         shocks currently in my random walker C{nw}. A constant term in the 
    271         log-likelihood is disregarded, C{sqrt(2*pi)}. 
    272         """ 
    273         mu = self.g * self.nw() 
     273        shocks currently in my random walker C{nw}. 
     274        """ 
     275        mu = +self.g * self.nw() 
    274276        sigma2 = 1 - self.g**2 
    275         logProbs = -(y - mu)**2 / (2*sigma2) - 0.5 * s.log(sigma2) 
     277        logProbs = -(y - mu)**2 / (2*sigma2) -  0.5 * s.log(2*s.pi*sigma2) 
    276278        return logProbs.sum() 
    277279 
    278     def logUniform(self, N): 
     280    def logUniform(self): 
    279281        """ 
    280282        Returns a log-uniform variate taken from a large array that is 
    …… 
    284286        """ 
    285287        logU = getattr(self, '_logU_array', []) 
    286         if self._logU_index >= len(logU): 
     288        if self._logU_index >= len(logU)-1: 
    287289            self._logU_index = -1 
    288290            # The efficient array operation 
    …… 
    296298    def likelihood(self, paramContainer, wiggle): 
    297299        """ 
    298         Call this method with a parameter object that defines (1) a set of 
    299         values for my parameters and (2) an array of volatility shocks to be 
    300         used as a starting point for another simulation draw of volatilities. 
     300        Call this method with a parameter object that defines: 
     301 
     302            1. a set of values for my parameters, and 
     303 
     304            2. a latent parameter consisting of an array of volatility shocks 
     305               to be used as a starting point for another simulation draw of 
     306               volatilities, and 
     307 
     308            3. a 1-D array of log-likelihoods for each observation time given 
     309               the volatility shocks. 
    301310 
    302311        A new 2-D array of log-volatilities will be drawn from a relatively 
    …… 
    309318 
    310319        The value of that log-likelihood at the end of the log-volatility 
    311         simulation is returned in a tuple, preceded by the volatility shocks 
     320        simulation is returned in a tuple, followed by the volatility shocks 
    312321        underlying the log-volatilities drawn at that point. 
    313322 
    …… 
    318327        what, making any call to this method with it a waste of computing time. 
    319328        """ 
    320         if paramContainer.v: 
    321             self.nw(paramContainer.v) 
     329        if len(paramContainer.latent): 
     330            self.nw(paramContainer.latent) 
    322331        self.setParamSequence(paramContainer.paramSequence()) 
    323332        # Obtain innovations as residuals of the reversed VAR(1) 
    …… 
    330339            # Compute the total log-likelihood conditional on the proposed 
    331340            # log-volatilities 
     341            #Lp = self.likelihoodOfInnovations(s.exp(-0.5*h)*x) 
    332342            Lp = self.likelihoodOfInnovations(self.decorrelate(x, h)) 
    333343            # Metropolis-Hastings accept/reject of the proposal 
    334344            if Lp > L or (Lp-L) > self.logUniform(): 
     345                if j > 0: 
     346                    print "\n%d -> %d" % (L, Lp) 
    335347                L = Lp 
    336348                vShocks = self.nw() 
    337         return vShocks, L 
     349            else: 
     350                print ".", 
     351        print "\n" 
     352        return L, vShocks 

  • projects/AsynCluster/tags/20080419/svpmc/params.py

    r150 r151  
    2727 
    2828 
     29SHAPE_MISMATCH = ValueError( 
     30    "shape mismatch: objects cannot be broadcast to a single shape") 
     31 
     32 
    2933PARAM_NAMES = [ 
    3034    # Approximately p*(3+2.5) parameters total, e.g., 125 for p=5 
    …… 
    5054 
    5155 
     56class FlexArray(object): 
     57    """ 
     58    I am an array of arbitrary objects, accessible in much the same manner as the 
     59    elements of a SciPy array object. Initially, all my objects are C{None}. 
     60 
     61    Call an instance of me with a hash array and you'll get a new instance with 
     62    just the subset of objects referenced by the array. 
     63 
     64    Call the instance with a single object ID and you'll get that object, just 
     65    like you'd get a numeric scalar by addressing a single element of a SciPy 
     66    array. 
     67 
     68    B{Caveat}: Unlike the behavior of SciPy arrays, Setting elements of a 
     69      subset or raveled version of an instance of me does I{not} set the 
     70      corresponding element of the original instance. 
     71     
     72    """ 
     73    def __init__(self, *shape): 
     74        self.shape = shape 
     75        self.N = 1 
     76        for dim in shape: 
     77            self.N *= dim 
     78        # Object list 
     79        self.O = [None] * self.N 
     80        # Index array 
     81        self.I = s.arange(self.N).reshape(*shape) 
     82 
     83    def __call__(self, I): 
     84        if isinstance(I, int): 
     85            return self.O[I] 
     86        newVersion = self.__class__(*I.shape) 
     87        for j, k in enumerate(I.ravel()): 
     88            newVersion.O[j] = self.O[k] 
     89        return newVersion 
     90 
     91    def __len__(self): 
     92        return self.N 
     93 
     94    def __iter__(self): 
     95        self.k = -1 
     96        return self 
     97 
     98    def next(self): 
     99        self.k += 1 
     100        if self.k < self.N: 
     101            return self(self.I[self.k]) 
     102        raise StopIteration 
     103     
     104    def __getitem__(self, key): 
     105        return self(self.I[key]) 
     106 
     107    def __setitem__(self, key, value): 
     108        I = self.I[key] 
     109        if isinstance(I, int): 
     110            self.O[I] = value 
     111        elif I.shape == s.shape(value): 
     112            if not isinstance(value, (list, tuple)): 
     113                value = value.ravel() 
     114            for j, k in enumerate(I.ravel()): 
     115                self.O[k] = value[j] 
     116        else: 
     117            raise SHAPE_MISMATCH 
     118 
     119    def ravel(self): 
     120        return self(self.I.ravel()) 
     121 
     122    def reshape(self, *shape): 
     123        return self(self.I.reshape(*shape)) 
     124 
     125    def __add__(self, other): 
     126        if other.shape != self.shape: 
     127            raise SHAPE_MISMATCH 
     128        newVersion = self(self.I) 
     129        for k, thisObj in enumerate(newVersion.O): 
     130            thisObj += other.O[k] 
     131        return newVersion 
     132 
     133    def call(self, methodName, *args, **kw): 
     134        """ 
     135        For each of my objects, call the method specified by I{methodName} with 
     136        any further args and keywords supplied. 
     137 
     138        Returns the result (must be a numeric scalar) of each call in an array 
     139        of the same shape as my object array. 
     140        """ 
     141        x = s.empty(self.shape) 
     142        for k in self.I.ravel(): 
     143            thisValue = getattr(self.O[k], methodName)(*args, **kw) 
     144            if not isinstance(thisValue, (int, long, float)): 
     145                raise ValueError( 
     146                    "Called method must return a numeric scalar value") 
     147            x.ravel()[k] = thisValue 
     148        return x 
     149 
    52150 
    53151class PriorContainer(object): 
    …… 
    61159        self.spec = fh.read() 
    62160        fh.close() 
    63      
    64  
    65161     
    66162 
    …… 
    140236 
    141237 
    142 class ParamArray(object): 
    143     """ 
    144     My items can be accessed in the same manner as those of a SciPy array 
    145     object. 
    146     """ 
    147     def __init__(self, iterable): 
    148         self.pcList = list(iterable) 
    149  
    150     def _get_indices(self): 
    151         return s.arange(len(self.pcList)) 
    152     I = property(_get_indices) 
    153      
    154     def __len__(self): 
    155         return len(self.pcList) 
    156      
    157     def __getitem__(self, key): 
    158         I = self.I[key] 
    159         if isinstance(I, int): 
    160             return self.pcList[I] 
    161         return self.__class__([self.pcList[x] for x in I]) 
    162  
    163     def __getslice__(self, i, j): 
    164         return self.__class__(self.pcList[i:j]) 
    165  
    166     def append(self, pc): 
    167         self.pcList.append(pc) 
    168  
    169          
    170  
    171238class ParameterContainer(Parameterized): 
    172239    """ 
    173240    """ 
    174     keyAttrs = {'v':[]} 
     241    keyAttrs = {'latent':[], 'pLogLikelihood':None, 'pLogProposal':None} 
    175242    paramNames = PARAM_NAMES 
    176243 
    177244    def __add__(self, other): 
    178         newVersion = self.__class__(v=self.v) 
     245        newVersion = self.__class__(latent=self.latent) 
    179246        newVersion.setParamSequence([ 
    180247            getattr(self, name) + getattr(other, name) 

  • projects/AsynCluster/tags/20080419/svpmc/pmc.py

    r150 r151  
    3333from asyncluster.master import jobs 
    3434 
    35 import model 
     35import model, params 
    3636from sample import resample 
    3737 
    …… 
    194194        """ 
    195195        k = 0 
    196         XP, W = params.ParamArray, [] 
     196        XP, W = [], [] 
    197197        wfd = defer.waitForDeferred( 
    198198            self.queue.call(self.proposals, X, v, series='proposals')) 
    …… 
    214214                W.append(L[If] + XPJ[1][Ic][If] - XPJ[2][Ic][If]) 
    215215        if XP: 
    216             yield s.row_stack(XP), s.concatenate(W) 
     216            arrayXP = params.FlexArray(len(XP)) 
     217            for k, thisXP in enumerate(XP): 
     218                arrayXP[k] = thisXP 
     219            yield arrayXP, s.concatenate(W) 
    217220        else: 
    218221            yield [], [] 

  • projects/AsynCluster/tags/20080419/svpmc/sample.py

    r150 r151  
    8787        if x is None: 
    8888            return self.x.copy() 
    89         if s.shape(x) != self.x.shape: 
     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: 
    9096            raise ValueError( 
    9197                "You can only update me with an equivalent-dimensioned array") 
    …… 
    134140        vector supplied will be re-used. 
    135141        """ 
    136         if correlations: 
    137             # Generate the correlation matrix and save in case not supplied 
    138             # next time 
    139             self.r = linalg.cholesky( 
    140                 self.covarMatrix(correlations), lower=True) 
     142        if correlations is not None: 
    141143            self.y = s.empty_like(self.x) 
     144            if correlations: 
     145                # Generate the correlation matrix and save in case not supplied 
     146                # next time 
     147                self.r = linalg.cholesky( 
     148                    self.covarMatrix(correlations), lower=True) 
     149            else: 
     150                self.r = None 
     151             
    142152        n = self.x.shape[0] 
    143         if s.shape(getattr(self, 'r', None)) != (n, n): 
     153        if n == 1: 
     154            return self.x.copy() 
     155        if not hasattr(self, 'r') or s.shape(self.r) != (n, n): 
    144156            raise ValueError( 
    145157                "Must have a [%d x %d] cross-correlation matrix defined" \ 

  • projects/AsynCluster/tags/20080419/svpmc/test/test_model.py

    r149 r151  
    2222import scipy as s 
    2323from scipy import stats, random 
     24from scipy.signal import signaltools 
    2425 
    2526from twisted.internet import reactor, defer 
    …… 
    2829from twisted_goodies.pybywire import pack, params 
    2930 
    30 import model 
     31import model, tseries 
    3132import util 
    3233 
    …… 
    150151                "'Uncorrelated' hypothesis fails with p=%f" % corrInfo[1]) 
    151152 
     153    def _check_corr(self, x, y, plot=True): 
     154        corrInfo = stats.spearmanr(x, y) 
     155        failed = abs(corrInfo[0]) < 0.5 and corrInfo[1] > 0.05 
     156        regressInfo = stats.linregress(x, y) 
     157        if (plot and VERBOSE) or failed: 
     158            m, b = regressInfo[:2] 
     159            sp = self.fig.add_subplot(111) 
     160            sp.plot(x, y, '.', x, m*x+b, '-') 
     161            sp.grid(True) 
     162            sp.set_title("r=%f, p=%f" % corrInfo) 
     163        if failed: 
     164            self.fail( 
     165                "Correlation is %f; " % corrInfo[0] +\ 
     166                "'Uncorrelated' hypothesis not disproved, p=%f" % corrInfo[1]) 
     167     
    152168    def _check_xcorr(self, dPeak): 
    153169        n = self.model.n 
    …… 
    354370 
    355371class Test_Model_other(Multivariate_BaseCase): 
     372    N = 100 
     373    corr = 0.0 
     374    params = { 
     375        'a': [0], 
     376        'b': [[0]], 
     377        'c': [], 
     378        'd': [0.01], 
     379        'e': [[0.8]], 
     380        'f': [], 
     381        'g': [corr], 
     382        } 
     383    for key, value in params.iteritems(): 
     384        params[key] = s.array(value) 
     385 
     386    class Mock_ParameterContainer(object): 
     387        latent = [] 
     388        def paramSequence(self): 
     389            return [self.params[x] for x in "abcdefg"] 
     390 
     391    @classmethod 
     392    def pcFactory(cls): 
     393        pc = cls.Mock_ParameterContainer() 
     394        pc.params = cls.params 
     395        return pc 
     396     
    356397    def setUp(self): 
    357         self.model = model.Model() 
    358  
     398        # Parameters 
     399        # Simulate a single time series of data per the parameters 
     400        self.iv = random.multivariate_normal( 
     401            [0, 0], [[1.0, self.corr], [self.corr, 1.0]], self.N).transpose() 
     402        h = signaltools.lfilter( 
     403            [1], [1.0, self.params['e'][0][0]], self.iv[1,:]) 
     404        self.x = s.exp(0.5*h) * self.iv[1,:] 
     405        # Instantiate a model to be tested and set its parameters 
     406        self.model = model.Model(tsList=[tseries.TimeSeries( 
     407            'test', data=s.column_stack([s.arange(self.N), self.x]))]) 
     408        for name, value in self.params.iteritems(): 
     409            setattr(self.model, name, value) 
     410     
    359411    def test_likelihoodOfInnovations(self): 
    360         self.fail("TODO") 
     412        self.model.nw(self.iv[1,:]) 
     413        pLog = self.model.likelihoodOfInnovations(self.iv[0,:]) 
     414        self.failUnlessAlmostEqual( 
     415            4.0 * s.exp(pLog/self.N) * s.sqrt(1-self.corr**2), 1.0, 0) 
     416        #   <+>   <-+-------------->   <-+---------------------------> 
     417        #    |      |                    |  
     418        #    |      |                    +- Divide by scaling term of 
     419        #    |      |                       joint density 
     420        #    |      | 
     421        #    |      +- Linear density, averaged 
     422        #    | 
     423        #    +--- Unscale by 1/2 total probability for each, squared 
    361424 
    362425    def test_logUniform(self): 
    363         self.fail("TODO") 
     426        N = 35000 
     427        import timeit 
     428        tObserved = timeit.Timer( 
     429            "M.logUniform()", 
     430            "import model\nM = model.Model()").timeit(N) 
     431        tBaseline = timeit.Timer( 
     432            "s.log(s.rand())", 
     433            "import scipy as s").timeit(N) 
     434        percent = 100 * tObserved/tBaseline 
     435        if VERBOSE: 
     436            print "The efficient method took %2d%% as long as the stupid way" \ 
     437                  % percent 
     438        self.failUnless(percent < 40) 
    364439         
    365440    def test_likelihood(self): 
    366         self.fail("TODO") 
     441        N_plot = 100 
     442        wiggle = 0.001 
     443        self.model.M = 10000 
     444        paramContainer = self.pcFactory() 
     445        # DEBUG 
     446        paramContainer.latent = 0.5*self.iv[1,:] + 0.5*s.randn(self.N) 
     447        L, vShocks = self.model.likelihood(paramContainer, wiggle) 
     448        #import pdb; pdb.set_trace() 
     449        if VERBOSE: 
     450            fig = self.fig 
     451            sp = fig.add_subplot(211) 
     452            sp.plot(self.iv[1,:N_plot]) 
     453            sp = fig.add_subplot(212) 
     454            sp.plot(vShocks[0][:N_plot]) 
     455            self._check_corr(self.iv[1,:], vShocks[0]) 
     456         

  • projects/AsynCluster/tags/20080419/svpmc/test/test_sample.py

    r148 r151  
    8989            if k % M == 0: 
    9090                yield x 
     91 
     92    def test_call(self): 
     93        walker = sample.NormalWalk(1, 100) 
     94        walker(s.ones(100)) 
     95        walker(s.ones((1,100))) 
     96        walker(s.ones((100,1))) 
    9197     
    9298    def test_univariate_bigJumps(self): 

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

    r150 r151  
    7070 
    7171    def likelihood(paramContainer, wiggle): 
    72         v, L = M.likelihood(paramContainer, wiggle) 
    73         return Packer(v, L)() 
     72        L, v = M.likelihood(paramContainer, wiggle) 
     73        return Packer(L, v)() 
    7474 
    7575    ########################################################################### 
    …… 
    112112    def likelihood(self, paramContainer, wiggle, remote=False, local=False): 
    113113        """ 
    114         Returns the log-likelihood of the parameters in the supplied 
     114        Computes the log-likelihood of the parameters in the supplied 
    115115        I{paramContainer} for my model. 
    116116 
    117         Updates in-place the container's sole latent parameter, an array of 
    118         volatility shocks after some MCMC involved with the likelihood 
    119         compution. 
     117        Updates the container in-place with the log-likelihood and an array of 
     118        volatility shocks on which the likelihood computation is based, after 
     119        undergoing some nested MCMC. 
     120 
     121        Returns a deferred that fires (with C{None}) when the likelihood 
     122        computation is done and the paramContainer has been updated. 
    120123        """ 
    121124        def gotPackedLikelihood(result): 
    …… 
    123126                # An error from the remote likelihood method call is treated as 
    124127                # infinitely low likelihood 
    125                 return -s.inf 
     128                paramContainer.pLogLikelihood = -s.inf 
    126129            up = pack.Unpacker(result) 
    127             paramContainer.v = up() 
    128             return up() 
     130            paramContainer.pLogLikelihood = up() 
     131            paramContainer.latent = up() 
    129132 
    130133        def gotLikelihood(result): 
    131             paramContainer.v = result[0] 
    132             return result[1] 
     134            paramContainer.pLogLikelihood = result[0] 
     135            paramContainer.latent = result[1] 
    133136         
    134137        if (remote or self.remoteMode) and not local: 
    …… 
    210213    def _get_var(self): 
    211214        if 'var' not in self.cache: 
    212             self.cache['var'] = VAR(self.tsList) 
     215            self.cache['var'] = VAR(self.y) 
    213216        return self.cache['var'] 
    214217    var = property(_get_var) 
    …… 
    254257            if hs != xs: 
    255258                raise ValueError("Need one log-volatility per innovation") 
    256             x *= s.exp(-0.5*h) 
     259            x = s.exp(-0.5*h) * x.copy() 
    257260        if 'decorrelate' not in self.cache: 
    258261            self.cache['decorrelate'] = linalg.inv( 
    …… 
    268271        The result will be the total log-likelihood of the innovations, given 
    269272        my volatility/innovation shock correlations C{g} and the volatility 
    270         shocks currently in my random walker C{nw}. A constant term in the 
    271         log-likelihood is disregarded, C{sqrt(2*pi)}. 
    272         """ 
    273         mu = self.g * self.nw() 
     273        shocks currently in my random walker C{nw}. 
     274        """ 
     275        mu = +self.g * self.nw() 
    274276        sigma2 = 1 - self.g**2 
    275         logProbs = -(y - mu)**2 / (2*sigma2) - 0.5 * s.log(sigma2) 
     277        logProbs = -(y - mu)**2 / (2*sigma2) -  0.5 * s.log(2*s.pi*sigma2) 
    276278        return logProbs.sum() 
    277279 
    278     def logUniform(self, N): 
     280    def logUniform(self): 
    279281        """ 
    280282        Returns a log-uniform variate taken from a large array that is 
    …… 
    284286        """ 
    285287        logU = getattr(self, '_logU_array', []) 
    286         if self._logU_index >= len(logU): 
     288        if self._logU_index >= len(logU)-1: 
    287289            self._logU_index = -1 
    288290            # The efficient array operation 
    …… 
    296298    def likelihood(self, paramContainer, wiggle): 
    297299        """ 
    298         Call this method with a parameter object that defines (1) a set of 
    299         values for my parameters and (2) an array of volatility shocks to be 
    300         used as a starting point for another simulation draw of volatilities. 
     300        Call this method with a parameter object that defines: 
     301 
     302            1. a set of values for my parameters, and 
     303 
     304            2. a latent parameter consisting of an array of volatility shocks 
     305               to be used as a starting point for another simulation draw of 
     306               volatilities, and 
     307 
     308            3. a 1-D array of log-likelihoods for each observation time given 
     309               the volatility shocks. 
    301310 
    302311        A new 2-D array of log-volatilities will be drawn from a relatively 
    …… 
    309318 
    310319        The value of that log-likelihood at the end of the log-volatility 
    311         simulation is returned in a tuple, preceded by the volatility shocks 
     320        simulation is returned in a tuple, followed by the volatility shocks 
    312321        underlying the log-volatilities drawn at that point. 
    313322 
    …… 
    318327        what, making any call to this method with it a waste of computing time. 
    319328        """ 
    320         if paramContainer.v: 
    321             self.nw(paramContainer.v) 
     329        if len(paramContainer.latent): 
     330            self.nw(paramContainer.latent) 
    322331        self.setParamSequence(paramContainer.paramSequence()) 
    323332        # Obtain innovations as residuals of the reversed VAR(1) 
    …… 
    330339            # Compute the total log-likelihood conditional on the proposed 
    331340            # log-volatilities 
     341            #Lp = self.likelihoodOfInnovations(s.exp(-0.5*h)*x) 
    332342            Lp = self.likelihoodOfInnovations(self.decorrelate(x, h)) 
    333343            # Metropolis-Hastings accept/reject of the proposal 
    334344            if Lp > L or (Lp-L) > self.logUniform(): 
     345                if j > 0: 
     346                    print "\n%d -> %d" % (L, Lp) 
    335347                L = Lp 
    336348                vShocks = self.nw() 
    337         return vShocks, L 
     349            else: 
     350                print ".", 
     351        print "\n" 
     352        return L, vShocks 

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

    r150 r151  
    2727 
    2828 
     29SHAPE_MISMATCH = ValueError( 
     30    "shape mismatch: objects cannot be broadcast to a single shape") 
     31 
     32 
    2933PARAM_NAMES = [ 
    3034    # Approximately p*(3+2.5) parameters total, e.g., 125 for p=5 
    …… 
    5054 
    5155 
     56class FlexArray(object): 
     57    """ 
     58    I am an array of arbitrary objects, accessible in much the same manner as the 
     59    elements of a SciPy array object. Initially, all my objects are C{None}. 
     60 
     61    Call an instance of me with a hash array and you'll get a new instance with 
     62    just the subset of objects referenced by the array. 
     63 
     64    Call the instance with a single object ID and you'll get that object, just 
     65    like you'd get a numeric scalar by addressing a single element of a SciPy 
     66    array. 
     67 
     68    B{Caveat}: Unlike the behavior of SciPy arrays, Setting elements of a 
     69      subset or raveled version of an instance of me does I{not} set the 
     70      corresponding element of the original instance. 
     71     
     72    """ 
     73    def __init__(self, *shape): 
     74        self.shape = shape 
     75        self.N = 1 
     76        for dim in shape: 
     77            self.N *= dim 
     78        # Object list 
     79        self.O = [None] * self.N 
     80        # Index array 
     81        self.I = s.arange(self.N).reshape(*shape) 
     82 
     83    def __call__(self, I): 
     84        if isinstance(I, int): 
     85            return self.O[I] 
     86        newVersion = self.__class__(*I.shape) 
     87        for j, k in enumerate(I.ravel()): 
     88            newVersion.O[j] = self.O[k] 
     89        return newVersion 
     90 
     91    def __len__(self): 
     92        return self.N 
     93 
     94    def __iter__(self): 
     95        self.k = -1 
     96        return self 
     97 
     98    def next(self): 
     99        self.k += 1 
     100        if self.k < self.N: 
     101            return self(self.I[self.k]) 
     102        raise StopIteration 
     103     
     104    def __getitem__(self, key): 
     105        return self(self.I[key]) 
     106 
     107    def __setitem__(self, key, value): 
     108        I = self.I[key] 
     109        if isinstance(I, int): 
     110            self.O[I] = value 
     111        elif I.shape == s.shape(value): 
     112            if not isinstance(value, (list, tuple)): 
     113                value = value.ravel() 
     114            for j, k in enumerate(I.ravel()): 
     115                self.O[k] = value[j] 
     116        else: 
     117            raise SHAPE_MISMATCH 
     118 
     119    def ravel(self): 
     120        return self(self.I.ravel()) 
     121 
     122    def reshape(self, *shape): 
     123        return self(self.I.reshape(*shape)) 
     124 
     125    def __add__(self, other): 
     126        if other.shape != self.shape: 
     127            raise SHAPE_MISMATCH 
     128        newVersion = self(self.I) 
     129        for k, thisObj in enumerate(newVersion.O): 
     130            thisObj += other.O[k] 
     131        return newVersion 
     132 
     133    def call(self, methodName, *args, **kw): 
     134        """ 
     135        For each of my objects, call the method specified by I{methodName} with 
     136        any further args and keywords supplied. 
     137 
     138        Returns the result (must be a numeric scalar) of each call in an array 
     139        of the same shape as my object array. 
     140        """ 
     141        x = s.empty(self.shape) 
     142        for k in self.I.ravel(): 
     143            thisValue = getattr(self.O[k], methodName)(*args, **kw) 
     144            if not isinstance(thisValue, (int, long, float)): 
     145                raise ValueError( 
     146                    "Called method must return a numeric scalar value") 
     147            x.ravel()[k] = thisValue 
     148        return x 
     149 
    52150 
    53151class PriorContainer(object): 
    …… 
    61159        self.spec = fh.read() 
    62160        fh.close() 
    63      
    64  
    65161     
    66162 
    …… 
    140236 
    141237 
    142 class ParamArray(object): 
    143     """ 
    144     My items can be accessed in the same manner as those of a SciPy array 
    145     object. 
    146     """ 
    147     def __init__(self, iterable): 
    148         self.pcList = list(iterable) 
    149  
    150     def _get_indices(self): 
    151         return s.arange(len(self.pcList)) 
    152     I = property(_get_indices) 
    153      
    154     def __len__(self): 
    155         return len(self.pcList) 
    156      
    157     def __getitem__(self, key): 
    158         I = self.I[key] 
    159         if isinstance(I, int): 
    160             return self.pcList[I] 
    161         return self.__class__([self.pcList[x] for x in I]) 
    162  
    163     def __getslice__(self, i, j): 
    164         return self.__class__(self.pcList[i:j]) 
    165  
    166     def append(self, pc): 
    167         self.pcList.append(pc) 
    168  
    169          
    170  
    171238class ParameterContainer(Parameterized): 
    172239    """ 
    173240    """ 
    174     keyAttrs = {'v':[]} 
     241    keyAttrs = {'latent':[], 'pLogLikelihood':None, 'pLogProposal':None} 
    175242    paramNames = PARAM_NAMES 
    176243 
    177244    def __add__(self, other): 
    178         newVersion = self.__class__(v=self.v) 
     245        newVersion = self.__class__(latent=self.latent) 
    179246        newVersion.setParamSequence([ 
    180247            getattr(self, name) + getattr(other, name) 

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

    r150 r151  
    3333from asyncluster.master import jobs 
    3434 
    35 import model 
     35import model, params 
    3636from sample import resample 
    3737 
    …… 
    194194        """ 
    195195        k = 0 
    196         XP, W = params.ParamArray, [] 
     196        XP, W = [], [] 
    197197        wfd = defer.waitForDeferred( 
    198198            self.queue.call(self.proposals, X, v, series='proposals')) 
    …… 
    214214                W.append(L[If] + XPJ[1][Ic][If] - XPJ[2][Ic][If]) 
    215215        if XP: 
    216             yield s.row_stack(XP), s.concatenate(W) 
     216            arrayXP = params.FlexArray(len(XP)) 
     217            for k, thisXP in enumerate(XP): 
     218                arrayXP[k] = thisXP 
     219            yield arrayXP, s.concatenate(W) 
    217220        else: 
    218221            yield [], [] 

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

    r150 r151  
    8787        if x is None: 
    8888            return self.x.copy() 
    89         if s.shape(x) != self.x.shape: 
     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: 
    9096            raise ValueError( 
    9197                "You can only update me with an equivalent-dimensioned array") 
    …… 
    134140        vector supplied will be re-used. 
    135141        """ 
    136         if correlations: 
    137             # Generate the correlation matrix and save in case not supplied 
    138             # next time 
    139             self.r = linalg.cholesky( 
    140                 self.covarMatrix(correlations), lower=True) 
     142        if correlations is not None: 
    141143            self.y = s.empty_like(self.x) 
     144            if correlations: 
     145                # Generate the correlation matrix and save in case not supplied 
     146                # next time 
     147                self.r = linalg.cholesky( 
     148                    self.covarMatrix(correlations), lower=True) 
     149            else: 
     150                self.r = None 
     151             
    142152        n = self.x.shape[0] 
    143         if s.shape(getattr(self, 'r', None)) != (n, n): 
     153        if n == 1: 
     154            return self.x.copy() 
     155        if not hasattr(self, 'r') or s.shape(self.r) != (n, n): 
    144156            raise ValueError( 
    145157                "Must have a [%d x %d] cross-correlation matrix defined" \ 

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

    r149 r151  
    2222import scipy as s 
    2323from scipy import stats, random 
     24from scipy.signal import signaltools 
    2425 
    2526from twisted.internet import reactor, defer 
    …… 
    2829from twisted_goodies.pybywire import pack, params 
    2930 
    30 import model 
     31import model, tseries 
    3132import util 
    3233 
    …… 
    150151                "'Uncorrelated' hypothesis fails with p=%f" % corrInfo[1]) 
    151152 
     153    def _check_corr(self, x, y, plot=True): 
     154        corrInfo = stats.spearmanr(x, y) 
     155        failed = abs(corrInfo[0]) < 0.5 and corrInfo[1] > 0.05 
     156        regressInfo = stats.linregress(x, y) 
     157        if (plot and VERBOSE) or failed: 
     158            m, b = regressInfo[:2] 
     159            sp = self.fig.add_subplot(111) 
     160            sp.plot(x, y, '.', x, m*x+b, '-') 
     161            sp.grid(True) 
     162          &nbs