Changeset 145

Show
Ignore:
Timestamp:
04/10/08 22:47:22 (9 months ago)
Author:
edsuom
Message:

Improvements to weave.Weaver; much better unit testing for VAR

Files:

Legend:

Unmodified
Added
Removed
Modified
Copied
Moved

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

    r143 r145  
    8383 
    8484 
    85 // Model.simulateVolatilities 
     85// Model.draw_h 
    8686//  
    8787// Supplied variables 
    8888// ---------------------------------------------------------------------------- 
    89 // h    Volatility values (initially empty), [pn] array 
     89// h    Volatility values (values to be overwritten), [pn] array 
    9090// v    Volatility shock values, [pn] array 
    9191// d    Volatility offset, [p] vector 

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

    r144 r145  
    266266     
    267267    def reverse(self, a, b): 
    268         return self.c( 
    269             'v', 'y', 'a', 'b', v=s.empty_like(self.y), a=a, b=b) 
     268        return self.c('v', 'y', 'a', 'b', v=s.empty_like(self.y)) 
    270269 
    271270    def forward(self, v, a, b): 
    272         return self.c( 
    273             'y', 'v', 'a', 'b', y=s.empty_like(v), v=v, a=a, b=b) 
     271        return self.c('y', 'v', 'a', 'b', y=s.empty_like(v)) 
    274272 
    275273 
    …… 
    340338 
    341339 
    342     #--- Methods -------------------------------------------------------------- 
    343  
    344          
    345                  
    346  
    347  
    348     def pDistribution(self, paramSequence, N, M=1): 
    349         """ 
    350         After application of the normalization and distribution parameters from 
    351         the supplied parameter sequence, returns a 1d array of I{N} points over 
    352         the range of my normalized data and a 1d array of probability densities 
    353         for each. 
    354  
    355         You can compute the probability densities over a multiple of the data 
    356         range by setting the keyword I{M} to something greater than 1. 
    357         """ 
    358         # TODO 
    359      
    360     def likelihood(self, paramSequence): 
    361         """ 
    362         Returns the log-likelihood of my normalized data given my distribution, 
    363         after application of the parameters from the supplied parameter 
    364         sequence. 
     340    #--- Supporting Methods --------------------------------------------------- 
     341 
     342    def draw_h(self, wiggle): 
     343        """ 
     344        Iterates my normal walkers by one step using the supplied fraction 
     345        I{wiggle} and returns a 2-D array of volatilities, one row for each of 
     346        my time series of observations. 
     347        """ 
     348        if not hasattr(self, 'h'): 
     349            self.h = s.empty_like(self.nw.x) 
     350        self.nw.step(wiggle) 
     351        return self.c('h', 'v', 'd', 'e', v=self.nw.correlate(self.f)) 
     352 
     353 
     354    #--- The API -------------------------------------------------------------- 
     355 
     356    def likelihood(self, wiggle, paramSequence): 
     357        """ 
     358        Returns a log-likelihood of my observations for one set of simulated 
     359        volatilies, given the supplied fractional I{wiggle} and the parameter 
     360        arrays in the supplied sequence I{paramSequence}. 
    365361 
    366362        This method returns the likelihood of the data given the parameters. It 
    367363        does not consider the prior probability of the parameters. You have to 
    368         do that yourself via a separate call to L{pPriors}, preferably before 
    369         calling this method. If the prior probability of any parameter is zero, 
    370         the posterior density for that parameter will also be zero, and a call 
    371         to this method will be unnecessary. 
    372         """ 
    373         # TODO 
     364        do that yourself, preferably before calling this method. If the prior 
     365        probability of any parameter is zero, the posterior density for that 
     366        parameter will also be zero, and a call to this method will be 
     367        unnecessary. 
     368        """ 
     369        for k, value in enumerate(paramSequence): 
     370            setattr(self, self.paramNames[k], value) 
     371        # Obtain innovations as residuals of the reversed VAR(1) 
     372        x = self.var.reverse(self.a, self.b) 
     373        # Draw a new set of simulated volatilities 
     374        h = self.draw_h(wiggle) 
     375        # Compute the conditional densities 
     376         

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

    r144 r145  
    4242 
    4343 
    44 // NormalWalk.__call__ 
     44// NormalWalk.step 
    4545// 
    4646// Supplied variables 
    4747// ---------------------------------------------------------------------------- 
    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 
     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 
    5252 
    5353int i, j, k; 
    …… 
    6161    for(k=0; k<m; k++) { 
    6262      // Generate a uniform, symmetric, random walk proposal 
    63       xp = X2(i,j) + W2(i,j) * (drand48() - 0.5); 
     63      xp = X2(i,j) + wiggle * (drand48() - 0.5); 
    6464      // Do the ratio in log space 
    6565      p_xp = -0.5 * pow(xp, 2); 
    …… 
    8383// Supplied variables 
    8484// ---------------------------------------------------------------------------- 
    85 // y    2-D array of correlated random values (initially empty) 
    86 // x    2-D array of independent random values 
    87 // p    2-D array of cross-correlations between values in each column of x 
     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 
    8888 
    8989int i, j, k; 
    9090double sum; 
    91   
    92 // For each column... 
    93 for(i=1; i<Nv[1]; i++) { 
    94   // For each time series... 
    95   for(i=0; i<Nv[0]; i++) { 
    96     // The offset plus the current shock... 
    97     sum = D1(i) + V2(i,j); 
    98     // ...plus the dot product for the VAR(1) term... 
    99     for(k=0; k<Nv[0]; k++) 
    100       sum += E2(i,k) * H2(k,j-1); 
    101     // ...is the modeled output 
    102     H2(i,j) = sum; 
     91 
     92// For each column of independent variates to be correlated... 
     93for(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; 
    103103  } 
    104104} 

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

    r144 r145  
    6363        x[:,0] = K*W[I0] 
    6464        a = s.less(x[:,0], 1.0).nonzero()[0] 
    65         x = self.c('x', 'a', 'b', x=x, a=a, b=s.setdiff1d(s.arange(0, K), a)) 
     65        x = self.c('x', 'a', 'b', b=s.setdiff1d(s.arange(0, K), a)) 
    6666        V = s.rand(N) 
    6767        RI = s.random.randint(0, K, N) 
    …… 
    7474    I run a 2-D array of values through a random walk within a zero-mean, unit 
    7575    variance normal distribution. 
    76      
    77     Call an instance of me with an array I{w} of fractional wiggle values in 
    78     the range 0.0 to +1.0, with one array element for each of my random 
    79     walkers. The elements of my value array will be perturbed accordingly and 
    80     updated array returned. 
    8176    """ 
    8277    m = 10 
    …… 
    8580        self.x = s.randn(rows, cols) 
    8681        self.p = s.zeros_like(self.x) 
     82        self.r = s.empty_like(self.x) 
    8783     
    88     def __call__(self, w): 
    89         if s.shape(w) != self.x.shape: 
    90             raise ValueError( 
    91                 "Wiggle array must have same dimensions as walker array") 
    92         return self.c('x', 'p', 'm', 'w', w=w) 
     84    def step(self, wiggle): 
     85        """ 
     86        Call with a fractional wiggle value in the range 0.0 to +1.0. The 
     87        elements of my value array will be perturbed via a random walk MCMC 
     88        step with a symmetric, uniform proposal over +/- half the wiggle value. 
     89 
     90        A reference to the updated array is returned. 
     91        """ 
     92        return self.c('x', 'p', 'm', 'wiggle') 
    9393 
    9494    def _covarMatrix(self, correlations): 
    …… 
    9898        """ 
    9999        i = 0 
    100         p = self.x.shape[0] 
    101         c = s.ones((p, p)) 
    102         for j in xrange(p-1): 
    103             for k in xrange(j+1, p): 
    104                 c[j, k] = c[k, j] = correlations[i] 
     100        n = self.x.shape[0] 
     101        cv = s.ones((n, n)) 
     102        for j in xrange(n - 1): 
     103            for k in xrange(j+1, n): 
     104                cv[j, k] = cv[k, j] = correlations[i] 
    105105                i += 1 
    106         return c 
     106        return cv 
    107107     
    108108    def correlate(self, correlations=None): 
    …… 
    112112        which are in the following form:: 
    113113 
    114             [r01, r02,..., r0p, r12,..., r1p,..] 
     114            [r01, r02,..., r0n, r12,..., r1n,..] 
    115115 
    116116        The values in each column of the result are correlated. The values from 
    …… 
    122122        """ 
    123123        if correlations: 
    124             self.c = linalg.cholesky( 
     124            # Generate the correlation matrix and save in case not supplied 
     125            # next time 
     126            self.r = linalg.cholesky( 
    125127                self._covarMatrix(correlations), lower=True) 
     128            self.y = s.empty_like(self.x) 
     129        n = self.x.shape[0] 
     130        if s.shape(getattr(self, 'r', None)) != (n, n): 
     131            raise ValueError( 
     132                "Must have a [%d x %d] cross-correlation matrix defined" \ 
     133                % (n, n)) 
     134        return self.c('y', 'x', 'r', y=s.empty_like(self.x)) 
    126135         
    127136         
     137         

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

    r143 r145  
    2121 
    2222import scipy as s 
    23 from scipy import random 
     23from scipy import stats, random 
    2424 
    2525from twisted.internet import reactor, defer 
    …… 
    3232import util 
    3333 
     34 
     35VERBOSE = True 
    3436 
    3537PRICES = """ 
    …… 
    168170 
    169171    def _plot(self, *vars): 
    170         N = len(vars) 
     172        if VERBOSE: 
     173            N = len(vars) 
     174            fig = self.fig 
     175            for k, var in enumerate(vars): 
     176                sp = fig.add_subplot(100*N+k+11) 
     177                sp.plot(var) 
     178                sp.grid(True) 
     179     
     180    def test_forward_indep_drift(self): 
     181        n = 1000 
     182        a = s.array([1.0, 1.0]) 
     183        v = s.zeros((2, n)) 
     184        for k in xrange(100): 
     185            b1, b2 = 0.8*s.rand(2) + 0.1 
     186            b = s.array([[b1, 0.0], [0.0, b2]]) 
     187            y = self.var.forward(v, a, b) 
     188            for k, bk in enumerate((b1, b2)): 
     189                self.failUnlessAlmostEqual(y[k,n/2:].mean(), 1.0/(1-bk), 3) 
     190 
     191    def test_forward_indep_shocks(self): 
     192        n = 10000 
     193        b1, b2 = 0.95, 0.8 
     194        a = s.array([0.0, 0.0]) 
     195        b = s.array([[b1, 0.0], [0.0, b2]]) 
     196        v = s.randn(2, n) 
     197        y = self.var.forward(v, a, b) 
     198        # IIR filtering 
    171199        fig = self.fig 
    172         for k, var in enumerate(vars): 
    173             sp = fig.add_subplot(100*N+k+11) 
    174             sp.plot(var) 
    175             sp.grid(True) 
    176  
    177     def test_forward_indep(self): 
    178         a = s.array([1.0, 1.0]) 
    179         b = s.array([[0.8, 0.0], [0.0, 0.2]]) 
    180         v = 0.1*s.randn(2, 100) 
    181         y = self.var.forward(v, a, b) 
    182         self._plot(y[0,:], y[1,:]) 
    183  
    184     def test_reverse_indep(self): 
    185         a = s.array([1.0, 1.0]) 
    186         b = s.array([[0.8, 0.0], [0.0, 0.2]]) 
    187         v1 = 0.1*s.randn(2, 100) 
    188         _y = self.var.y 
    189         self.var.y = self.var.forward(v1, a, b) 
    190         v2 = self.var.reverse(a, b) 
    191         self._plot(v1[0,:], v2[0,:], v1[1,:], v2[1,:]) 
    192         self.var.y = _y 
     200        for k, bk in enumerate((b1, b2)): 
     201            sp = fig.add_subplot(221 + 2*k) 
     202            pxx = sp.psd(y[k,n/2:], NFFT=32)[0] 
     203            ratio = ( pxx[-1]/pxx[0] ) / ( 2*((1-bk)/(1+bk))**2 ) 
     204            self.failUnless(ratio > 0.4, "Too little rolloff in PSD") 
     205            self.failUnless(ratio < 1.5, "Too much rolloff in PSD") 
     206            sp = fig.add_subplot(221 + 2*k + 1) 
     207            sp.plot(y[k,:]) 
     208        # Uncorrelated 
     209        r = stats.spearmanr(y[0,:], y[1,:])[0] 
     210        self.failUnless(abs(r) < 0.1, "Doesn't appear uncorrelated") 
    193211 
    194212    def test_forward_xcorr(self): 
     213        n = 1000 
    195214        a = s.array([0.0, 0.0]) 
    196215        b = s.array([[0.0, 0.5], [0.0, 0.0]]) 
    197         v = 0.1*s.randn(2, 100) 
     216        v = 0.1*s.randn(2, n) 
    198217        y = self.var.forward(v, a, b) 
    199         self._plot(v[0,:], v[1,:], y[0,:]) 
    200  
    201     def test_reverse_xcorr(self): 
    202         a = s.array([1.0, 1.0]) 
    203         b = s.array([[0.01, 0.5], [0.01, 0.01]]) 
    204         v1 = 0.1*s.randn(2, 100) 
     218        # Correlation 
     219        nd = 6 
     220        d = s.arange(-nd, nd+1) 
     221        r = s.array( 
     222            [s.correlate(y[0,nd+lag:n-nd+lag], y[1,nd:n-nd])[0] for lag in d]) 
     223        lagOneMetric = r[nd+1] / max(s.concatenate([r[:nd+1], r[nd+1:]])) 
     224        self.failUnless(lagOneMetric < 10) 
     225        if VERBOSE: 
     226            sp = self.fig.add_subplot(111) 
     227            sp.plot(d, r, 'o-') 
     228            sp.grid(True) 
     229 
     230    def test_reverse(self): 
     231        n = 200 
     232        t = 2*s.pi*s.arange(n)/n 
     233        v1 = 0.1*s.randn(2,n) 
    205234        _y = self.var.y 
    206         self.var.y = self.var.forward(v1, a, b) 
    207         v2 = self.var.reverse(a, b) 
    208         for k in (0,1): 
    209             self._plot(v1[0,:], v1[1,:], self.var.y[k,:], v2[k,:]) 
     235        for j in xrange(100): 
     236            a = s.randn(2) 
     237            # Keep all coefficients small enough to ensure stability 
     238            b = 0.6*s.rand(2,2) 
     239            # Forward... 
     240            self.var.y = self.var.forward(v1, a, b) 
     241            # ...then reverse 
     242            v2 = self.var.reverse(a, b) 
     243            # Now check 
     244            for k in (0,1): 
     245                diff = s.absolute(v1[k,:] - v2[k,:]) 
     246                if max(diff) > 1E-1: 
     247                    table = "\n".join([ 
     248                        ", ".join(["%6.5f" % x for x in row]) 
     249                        for row in diff.reshape(20, n/20)]) 
     250                    self.fail( 
     251                        "Run #%d: excess round trip error:\n\n%s" \ 
     252                        % (j+1, "b:\n%s\n\n|diff|:\n%s" % (str(b), table))) 
    210253        self.var.y = _y 
     254 

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

    r144 r145  
    8282            sp.hist(x, bins=20) 
    8383 
     84    def walkTwo(self, N, p=2): 
     85        M = 20 
     86        NM = N*M 
     87        for k in xrange(NM): 
     88            x = self.walker.step(0.5) 
     89            if k % M == 0: 
     90                yield x 
     91     
    8492    def test_univariate_bigJumps(self): 
    8593        walker = sample.NormalWalk(1, 1) 
    86         w = s.array([[1.0]]) 
    87         x = s.array([walker(w)[0,0] for k in xrange(200)]) 
     94        x = s.array([walker.step(1.0)[0,0] for k in xrange(200)]) 
    8895        self.check(x) 
    8996 
    9097    def test_univariate_smallJumps(self): 
    9198        walker = sample.NormalWalk(1, 1) 
    92         w = s.array([[0.1]]) 
    93         x = s.array([walker(w)[0,0] for k in xrange(20000)]) 
     99        x = s.array([walker.step(0.1)[0,0] for k in xrange(20000)]) 
    94100        self.check(x[::100]) 
    95101             
    96102    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) 
     103        N = 100 
     104        x1 = s.empty(N) 
     105        x2 = s.empty(N) 
     106        self.walker = sample.NormalWalk(2, 2) 
     107        for k, x in enumerate(self.walkTwo(N)): 
     108            x1[k] = x[0,0] 
     109            x2[k] = x[1,1] 
     110        self.check(x1) 
     111        self.check(x2) 
     112        self.failUnless(stats.pearsonr(x1, x2)[1] > 0.05) 
     113        if VERBOSE: 
     114            self.fig.add_subplot(111).plot(x1, x2) 
    109115         
    110116    def test_covarMatrix_2x2(self): 
    …… 
    132138            [0.14, 0.24, 0.34, 1.00]]) 
    133139        self.failUnless(s.equal(x, y).all()) 
     140 
     141    def test_correlate(self): 
     142        def vector(row, col): 
     143            return s.array([x[row, col] for x in yList]) 
     144         
     145        N = 500 
     146        yList = [] 
     147        self.walker = sample.NormalWalk(2, 2) 
     148        self.failUnlessEqual(self.walker.correlate([0.5]).shape, (2,2)) 
     149        for k, x in enumerate(self.walkTwo(N)): 
     150            y = self.walker.correlate() 
     151            yList.append(y.copy()) 
     152        # Correlated values in columns 
     153        for col in (0,1): 
     154            y1 = vector(0,col) 
     155            y2 = vector(1,col) 
     156            self.check(y1) 
     157            self.check(y2) 
     158            corrInfo = stats.pearsonr(y1, y2) 
     159            self.failUnless( 
     160                abs(corrInfo[0] - 0.5) < 0.1, 
     161                "Expected correlation of 0.5, not %f" % corrInfo[0]) 
     162            if VERBOSE: 
     163                self.fig.add_subplot(111).plot(y1, y2) 
     164        # Uncorrelated values in rows 
     165        for row in (0,1): 
     166            y1 = vector(row,0) 
     167            y2 = vector(row,1) 
     168            self.check(y1) 
     169            self.check(y2) 
     170            corrInfo = stats.pearsonr(y1, y2) 
     171            self.failUnless( 
     172                abs(corrInfo[0]) < 0.1, 
     173                "Expected correlation of 0, not %f" % corrInfo[0]) 
     174            if VERBOSE: 
     175                self.fig.add_subplot(111).plot(y1, y2) 

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

    r140 r145  
    5656 
    5757    def __call__(self, x): 
    58         return self.c('x', x=s.array([x]))[0] 
     58        return self.c('x')[0] 
    5959     
    6060    def foo(self, x): 
    …… 
    6868    def test_c_alone(self): 
    6969        self.weaver._code = {'support':None, 'foo':CODE} 
    70         self.failUnlessEqual(self.weaver.foo(1), 1164) 
     70        self.failUnlessEqual(self.weaver.foo(s.array([1])), 1164) 
    7171 
    7272    def test_parseCode(self): 
    …… 
    8080    def test_c_privateMethod(self): 
    8181        self.weaver._code = self.weaver._parseCode(CODE) 
    82         self.failUnlessEqual(self.weaver(7), 107) 
     82        self.failUnlessEqual(self.weaver(s.array([7])), 107) 

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

    r140 r145  
    1717 
    1818""" 
    19 Miscellaneous utility stuff 
     19Convenience class for doing inline C with SciPy.weave. 
     20 
     21You can import from this module instead of scipy.weave. Then you get the 
     22convenience class L{Weaver} in addition to the good stuff in scipy.weave. 
    2023""" 
    2124 
    …… 
    8083        Supply the names of variables to use as arguments in order as 
    8184        arguments. The variables can have their values set via keywords; 
    82         otherwise they will be retrieved as instance attributes. 
     85        otherwise they will be retrieved from the caller's namespace or, 
     86        failing that, from my instance's namespace as C{self.<varname>}. 
    8387 
    8488        A reference to the first variable is returned. 
    8589        """ 
    86         for varName in args: 
    87             value = kw.get(varName, getattr(self, varName, None)) 
    88             exec "%s = value" % varName 
    89         methodName = inspect.currentframe().f_back.f_code.co_name 
     90        frame = inspect.currentframe() 
     91        methodName = frame.f_back.f_code.co_name 
     92        callerVarNames = frame.f_back.f_locals.keys() 
     93        try: 
     94            for varName in args: 
     95                if varName in kw: 
     96                    value = kw[varName] 
     97                elif varName in callerVarNames: 
     98                    value = frame.f_back.f_locals[varName] 
     99                elif hasattr(self, varName): 
     100                    value = getattr(self, varName) 
     101                else: 
     102                    raise AttributeError( 
     103                        "Variable '%s' not defined " % varName +\ 
     104                        "in keyword, caller, or instance namespace") 
     105                exec "%s = value" % varName 
     106        finally: 
     107            # Avoid cycle problems, per Python Library Reference 3.11.4 
     108            del frame 
    90109        inline( 
    91110            self.code[methodName],