Changeset 140

Show
Ignore:
Timestamp:
04/07/08 22:25:54 (9 months ago)
Author:
edsuom
Message:

Got cool & convenient Weaver mixin written and tested; refactored svpmc.sample to use it; working on model.VAR

Files:

Legend:

Unmodified
Added
Removed
Modified
Copied
Moved

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

    r135 r140  
    7070        return pack.packwrap(method)(*possiblyPackedArgs) 
    7171 
    72     def likelihood(packedVector): 
    73         paramVector = pack.Unpacker(packedVector)() 
    74         L = M.likelihood(paramVector) 
     72    def likelihood(packedParams): 
     73        L = M.likelihood(list(pack.Unpacker(packedParams))) 
    7574        return str(float(L)) 
    7675 
    …… 
    148147        return d 
    149148 
    150     def likelihood(self, paramVector, remote=False, local=False): 
    151         """ 
    152         Returns the log likelihood of the supplied parameter vector for my 
     149    def likelihood(self, paramSequence, remote=False, local=False): 
     150        """ 
     151        Returns the log likelihood of the supplied parameter sequence for my 
    153152        model. 
    154153        """ 
    …… 
    161160         
    162161        if (remote or self.remoteMode) and not local: 
    163             packedVector = Packer(*paramVector) 
     162            packedParams = Packer(*paramSequence) 
    164163            d = self.client.run( 
    165                 'likelihood', packedVector(), timeout=self.timeout) 
     164                'likelihood', packedParams(), timeout=self.timeout) 
    166165            d.addCallback(gotLikelihood) 
    167166        else: 
    168             d = self.queue.call(self.modelObj.likelihood, paramVector) 
     167            d = self.queue.call(self.modelObj.likelihood, paramSequence) 
    169168        return d.addErrback(self._oops) 
    170169 
    …… 
    173172    """ 
    174173    You must define the name of an underlying dist object from the scipy.stats 
    175     L{distributions} module as the keyword 'distName'. 
    176  
    177     For all underlying dist objects, you must supply I{loc} and I{scale} as 
    178     keywords. Some of the objects also require I{a} and I{b}. 
    179     """ 
    180     paramNames = ['name', 'loc', 'scale', 'a', 'b'] 
     174    L{distributions} module via my parameter I{dname}. 
     175 
     176    For all underlying dist objects, you must also specify I{loc} and I{scale} as 
     177    additional parameters. Some of the objects also require I{a} and I{b}. 
     178    """ 
     179    paramNames = ['dname', 'loc', 'scale', 'a', 'b'] 
    181180     
    182181    def _get_distObj(self): 
    183182        if 'dist' not in self.cache: 
    184             if not hasattr(distributions, self.name): 
     183            if not hasattr(distributions, self.dname): 
    185184                raise ImportError( 
    186185                    "No distribution '%s' in scipy.stats.distributions" \ 
    187                     % self.name) 
    188             distClass = getattr(distributions, self.name) 
     186                    % self.dname) 
     187            distClass = getattr(distributions, self.dname) 
    189188            args = [] 
    190189            for xName in ('a', 'b'): 
    …… 
    245244 
    246245 
     246class VAR(Weaver): 
     247    """ 
     248    Construct me with a sequence of one or more time-series objects containing 
     249    observation data. 
     250 
     251    Then you can call me with a C{p x 1} drift vector I{a} and a C{p x p} 
     252    VAR(1) cross-correlation matrix I{b}, to get a C{p x n} vector of residuals 
     253    from my observation data after application of the VAR(1) model:: 
     254     
     255        y[:,t] = a + b*y[:,t-1] 
     256 
     257    @ivar y: Observations from I{n} intersecting dates of I{p} time series, 
     258      C{p,n} array. 
     259     
     260    """ 
     261    def __init__(self, tsList): 
     262        p = len(tsList) 
     263        tsObj = tsList[0] 
     264        self.y = s.empty(p, len(tsObj)) 
     265        for k in xrange(len(tsList)): 
     266            if k > 0: 
     267                # No need to intersect first object with itself 
     268                tsObj = tsObj.intersect(tsList[k]) 
     269            self.y[k,:] = tsObj() 
     270 
     271    def reverse(self, 'a', 'b'): 
     272        return self.c('r', 'y', 'a', 'b', r=s.empty_like(self.y), a=a, b=b) 
     273 
     274 
    247275class Model(Parameterized): 
    248276    """ 
    249     Stochastic Model for fitting a C{Kx1} vector of random variables to a 
    250     C{KxN} matrix of observations. 
    251     """ 
    252     keyAttrs = {'N_normParams':None, 'data':None, 'distribution':None} 
    253  
    254     def N_params(self): 
    255         if not hasattr(self, '_Np'): 
    256             Np = self.N_normParams + self.distribution.N_params() 
    257             self._Np = Np 
    258         return self._Np 
    259  
    260     def setParamVector(self, paramVector): 
    261         """ 
    262         Setting my parameters means setting my normalization parameters and 
    263         then setting my distribution with the parameters that are left. 
    264         """ 
    265         if len(paramVector) != self.N_params(): 
    266             raise ValueError("You must supply the exact number of parameters") 
    267         self.normParams = paramVector[:self.N_normParams] 
    268         self.distribution.setParamVector(paramVector[self.N_normParams:]) 
    269  
    270     def denormalize(self, values): 
    271         """ 
    272         Override this to do more sophisticated denormalization than mere drift 
    273         superposition. 
    274         """ 
    275         if len(self.normParams): 
    276             return values + sum(self.normParams) 
    277         return values 
     277    I am a sophisticated econometric model wherein the residuals of a 
     278    linear-drift, VAR(1) process follow a multivariate stochastic volatility 
     279    process with correlated volatility and return shocks. 
     280     
     281    """ 
     282    keyAttrs = {'tsList':None} 
     283    paramNames = ['a', 'b', 'c', 'd', 'e', 'f'] 
     284    support = resource_string(__name__, "%s_support.c" % name) 
     285 
     286    def _get_residualizer(self): 
     287        if not hasattr(self, '_residualizer'): 
     288            self._residualizer = Residualizer(self.tsList) 
     289        return self._residualizer 
     290    residualizer = property(_get_residualizer) 
     291     
     292    def dataToResiduals(self, a, b): 
     293        """ 
     294        """ 
     295        r = self.y 
     296        weave.inline( 
     297            """ 
     298             
     299            """, 
     300            ['r', 'a', 'b'], 
     301            extra_compile_args=['-w'], support_code=self.support) 
     302 
     303    def simulateVolatilities(self, svParams): 
     304        """ 
     305        """ 
     306 
     307    def  
    278308 
    279309    def rDistribution(self, N): 
    …… 
    320350        return Y, self.distribution.pdf(Yn) 
    321351     
    322     def likelihood(self, X): 
     352    def likelihood(self, paramSequence): 
    323353        """ 
    324354        Returns the log-likelihood of my normalized data given my distribution, 

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

    r139 r140  
    1515   this program; if not, write to the Free Software Foundation, Inc., 51 
    1616   Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA 
    17  
    18  
    19    Inline code for sample.NormalWalk 
    2017*/ 
    2118 
    2219 
     20// Resampler.__call__ 
     21// 
    2322// Supplied variables 
    24 // 'x', 'p', 'm', 'w' 
     23// ---------------------------------------------------------------------------- 
     24// 'x', 'a', 'b' 
    2525 
     26int ja = Na[0] - 1; 
     27int jb = Nb[0] - 1; 
     28int ka, kb; 
     29while (ja >= 0 && jb >= 0) { 
     30  ka = A1(ja); 
     31  kb = B1(jb); 
     32  X2(ka,1) = kb; 
     33  X2(kb,0) = X2(kb,0) + X2(ka,0) - 1.0; 
     34  if (X2(kb,0) < 1) { 
     35    A1(ja) = kb; 
     36    jb--; 
     37  } else { 
     38    ja--; 
     39  } 
     40} 
     41 
     42 
     43 
     44// NormalWalk.__call__ 
     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// w    2-D array of wiggle values for the uniform random walk proposals 
    2652 
    2753int i, j, k; 
    2854double xp, dp, p_xp; 
    2955 
     56// For each row... 
    3057for(i=0; i<Nx[0]; i++) { 
     58  // For each column... 
    3159  for(j=0; j<Nx[1]; j++) { 
     60    // For each oversampling iteration... 
    3261    for(k=0; k<m; k++) { 
     62      // Generate a uniform, symmetric, random walk proposal 
    3363      xp = X2(i,j) + W2(i,j) * (drand48() - 0.5); 
     64      // Do the ratio in log space 
    3465      p_xp = -0.5 * pow(xp, 2); 
    3566      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. 
    3670      if(dp > 0  || log(drand48()) < dp) { 
     71        // Update current value and its log probability when proposal accepted 
    3772        X2(i,j) = xp; 
    3873        P2(i,j) = p_xp; 

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

    r139 r140  
    2222import scipy as s 
    2323from scipy import random 
    24 from scipy.stats import distributions as dists 
    2524 
    26 from util import Weaver 
     25from weave import Weaver 
    2726 
    2827 
    …… 
    8180    updated array returned. 
    8281    """ 
    83     m = 30 
     82    m = 10 
    8483     
    8584    def __init__(self, rows, cols): 

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

    r139 r140  
    5353    def test_norm(self): 
    5454        # Now for real 
    55         N = 1000 
     55        N = 2000 
    5656        X = s.linspace(-6, +6, N) 
    5757        normDistObj = stats.distributions.norm() 
    …… 
    9191        walker = sample.NormalWalk(1, 1) 
    9292        w = s.array([[0.1]]) 
    93         x = s.array([walker(w)[0,0] for k in xrange(30000)]) 
    94         self.check(x[::30]) 
     93        x = s.array([walker(w)[0,0] for k in xrange(20000)]) 
     94        self.check(x[::100]) 
    9595             
     96    def test_multivariate(self): 
     97        N, T = 10, 20000 
     98        walker = sample.NormalWalk(N, 2) 
     99        w = s.column_stack([0.5*s.ones(N), 0.3*s.ones(N)]) 
     100        x1 = s.empty(T) 
     101        x2 = s.empty(T) 
     102        for k in xrange(T): 
     103            x = walker(w) 
     104            x1[k] = x[1,0] 
     105            x2[k] = x[8,1] 
     106        self.check(x1[::100]) 
     107        self.check(x2[::100]) 
     108        self.fig.add_subplot(111).plot(x1, x2) 
    96109         

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

    r133 r140  
    3333 
    3434from twisted.trial import unittest 
    35  
    36 for packageName in ('model',): 
    37     try: 
    38         exec "import %s" % packageName 
    39     except: 
    40         print "Error loading package '%s'..." % packageName 
    4135 
    4236 

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

    r137 r140  
    2020""" 
    2121 
     22import inspect, re 
    2223from pkg_resources import resource_string 
    23  
    24 from scipy import weave 
     24from scipy.weave import * 
    2525 
    2626 
    2727class Weaver(object): 
    2828    """ 
    29     Subclass from me to conveniently use a separate C code files for support 
    30     and an inline operation. 
     29    Subclass from me to conveniently use a separate C code file for inline 
     30    operations. 
    3131 
    32     The support file is for the entire module in which my sublcass resides, and 
    33     must be present in that module's directory under the name 
    34     C{<module>_support.c}. 
     32    The code file is for the entire module in which my sublcass resides, and 
     33    must be present in that module's directory under the name C{<module>.c}. It 
     34    is parsed into sections with headings as inline C comments in the following 
     35    format:: 
    3536 
    36     The inline code file is just for the subclass, and must be present in the 
    37     containing module's directory under the name 
    38     C{<module>_<classname>_support.c}. 
     37        // <subclass>.<method in which 'c' base class method is called> 
     38 
     39    The subclass name must start with a capital letter and the method name with 
     40    a lowercase letter (possibly followed by one or two underscores), as is 
     41    conventional. There must not be anything else but whitespace on the comment 
     42    line. 
    3943    """ 
     44    def _parseCode(self, text): 
     45        re_delimiter = re.compile( 
     46            "//\s+([A-Z][A-Za-z0-9_]+)\.(_{0,2}[a-z][A-Za-z0-9_]+)\s*$") 
     47        cName = self.__class__.__name__ 
     48        key = None 
     49        lines = [] 
     50        result = {'support':None} 
     51        for line in text.split("\n"): 
     52            m = re_delimiter.match(line) 
     53            if m is None: 
     54                if key: 
     55                    lines.append(line) 
     56                continue 
     57            if lines: 
     58                result[key] = "\n".join(lines) 
     59                lines = [] 
     60            if m.group(1) == cName: 
     61                key = m.group(2) 
     62            else: 
     63                key = None 
     64        if lines: 
     65            result[key] = "\n".join(lines) 
     66        return result 
     67     
    4068    def _get_code(self): 
    4169        if not hasattr(self, '_code'): 
    42             mName = self.__class__.__module__ 
    43             cName = self.__class__.__name__ 
    44             self._code = [ 
    45                 resource_string(mName, "%s_%s.c" % (mName, cName)), 
    46                 resource_string(mName, "%s_support.c" % mName)] 
     70            self._code = self._parseCode( 
     71                resource_string(__name__, "%s.c" % self.__class__.__module__)) 
    4772        return self._code 
    4873    code = property(_get_code) 
    …… 
    5075    def c(self, *args, **kw): 
    5176        """ 
    52         Call this method to run my inline code, supplying the names of 
    53         variables to use as arguments in order as arguments. The variables can 
    54         have their values set via keywords; otherwise they will be retrieved as 
    55         instance attributes. 
     77        Call this method to run the inline code for the calling method of the 
     78        subclass. 
     79 
     80        Supply the names of variables to use as arguments in order as 
     81        arguments. The variables can have their values set via keywords; 
     82        otherwise they will be retrieved as instance attributes. 
    5683 
    5784        A reference to the first variable is returned. 
    …… 
    6087            value = kw.get(varName, getattr(self, varName, None)) 
    6188            exec "%s = value" % varName 
    62         inline, support = self.code 
    63         weave.inline( 
    64             inline, args, extra_compile_args=['-w'], support_code=support) 
     89        methodName = inspect.currentframe().f_back.f_code.co_name 
     90        inline( 
     91            self.code[methodName], 
     92            args, extra_compile_args=['-w'], support_code=self.code['support']) 
    6593        return locals()[args[0]] 
    6694