Changeset 153

Show
Ignore:
Timestamp:
04/21/08 21:26:43 (9 months ago)
Author:
edsuom
Message:

Working on Model.likelihood, some really hairy C code

Files:

Legend:

Unmodified
Added
Removed
Modified
Copied
Moved

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

    r152 r153  
    8787// Supplied variables 
    8888// ---------------------------------------------------------------------------- 
    89 // L    Log-likelihood of innovations, given params & simulated volatilities 
     89// p    Log-likelihoods [2n] array (overwritten) 
    9090 
    9191//--- Supplied arrays --- 
    92 // z    Independent normal variates, [pn] array 
     92// z    Independent normal variates for log-volatilities, [pn] array (updated) 
     93// h    Log-volatilities, [pn] (overwritten) 
    9394// x    Innovation values, [pn] array 
    94 // r    Concurrent cross-correlations between volatility shocks, [pp] 
    95 // s    Inverse, concurrent cross-correlations between innovation shocks, [pp] 
    96  
    97 //--- Externally declared scratchpad --- 
    98 // v    Volatility shock values, [p] vector for one multivariate sample 
    99 // h    Current and previous log-volatilities, [p2] array 
    100 // y    Volatility-normalized and decorrelated innovation values, [p2] array 
     95// rvc  Concurrent cross-correlations between volatility shocks, [pp] 
     96// rvi  Inverse, concurrent cross-correlations between volatility shocks, [pp] 
     97// rii  Inverse, concurrent cross-correlations between innovation shocks, [pp] 
     98// y    Miscellaneous multivariate scratchpad, [p8] array (overwritten) 
    10199 
    102100//--- Model parameters --- 
    …… 
    108106// n0   Number of volatility draws per likelihood evaluation.   
    109107// n3   Number of random-walk normal draws per volatility draw. 
    110 // wiggle = Uniform proposal +/- extent 
    111  
    112  
    113 // THIS IS WHERE THE VAST MAJORITY OF THE CPU TIME IS EXPENDED! 
     108// sigma = Standard deviation of random-walk normal proposal 
     109 
     110 
     111// The log-likelihood array is used as follows, by column: 
     112// ---------------------------------------------------------------------------- 
     113// 0    Log Pr(x[t]|hp[t]), the sum of which for all t is the log-likelihood of 
     114//      the innovation values given the simulated log-volatilities 
     115// 1    Log Pr(hp[t]|x[t],h[t+1]), which is used in each M-H accept/reject step 
     116//      for the log-volatility simulation 
     117 
     118 
     119// The multivariate scratchpad y is used as follows, by column: 
     120// ---------------------------------------------------------------------------- 
     121// 0    2*(1-sigma**2) for bivariate normal probability density function 
     122// 1    0.5*log(pi*2*sigma**2) offset for bivariate normal pdf 
     123// 2    Independent normal variate random-walk proposals 
     124// 3    Correlated volatility shocks 
     125// 4    Innovation values after inverse-scaling by log-volatility 
     126// 5    Inverse-scaled innovation values after decorrelation 
     127// 6    Next-sample log-volatility values after VAR term removal 
     128// 7    VAR-term removed log-volatility values after decorrelation 
     129 
     130 
     131// THIS CODE IS WHERE THE VAST MAJORITY OF THE TOTAL CPU TIME IS EXPENDED! 
    114132 
    115133 
    116134int km, kh; 
    117135int k0, k1, k2, k3; 
    118 double sum, normp, dL, L_normp, Lp; 
     136double sum, normp, dL, L0, L1; 
     137 
     138double wiggle = sigma * pow(n3, -0.5); 
     139 
     140 
     141// Initialize  constants for bivariate normal probability density function 
     142for(k2=0; k2<Nz[0]; k2++) { 
     143  Y2(k2, 0) = 2*(1 - pow(G1(k2), 2)); 
     144  Y2(k2, 1) = 0.5 * log(3.1415926535897931 * Y2(k2, 0)); 
     145} 
    119146 
    120147 
    …… 
    122149for(k0=0; k0<n0; k0++) { 
    123150 
    124   // The log-likelihood of the innovations given the log-volatility proposal to 
    125   // be constructed 
    126   Lp = 0; 
    127  
    128   // Initialize log-volatility scratchpad: the "previous" volatility value 
    129   // for the first multivariate sample is just the offset 
    130   for(k2=0; k2<Nz[0]; k2++) 
    131       H2(k2, 1) = D1(i); 
    132  
    133151  // for each multivariate sample... 
    134152  for(k1=0; k1<Nz[1]; k1++) { 
    135153 
    136     // Update this column of independent random-walk normal variates to produce 
    137     // a multivariate i.i.d. normal proposal. It is this on which we build a 
    138     // multivariate log-volatility proposal for this multivariate sample. 
    139  
    140     // For each time series... 
    141     for(k2=0; k2<Nz[0]; k2++) { 
    142       // Get the log-likelihood of the current value 
    143       L_norm = -0.5 * pow(Z2(k2, k1), 2); 
    144       // For each oversampling iteration... 
    145       for(k3=0; k3<n3; k3++) { 
    146         // Generate a uniform, symmetric, random walk proposal 
    147         normp = Z2(k2, k1) + wiggle * (drand48() - 0.5); 
    148         // Do the ratio in log space 
    149         L_normp = -0.5 * pow(normp, 2); 
    150         dL = L_normp - L_norm; 
    151         // The uniform random variate thus needs to be in logspace as well, but 
    152         // it's not always even needed. Thus the log operation is done just 
    153         // once, and only some of the time. 
    154         if(dL > 0  || log(drand48()) < dL) { 
    155           // Update current value and its log probability when proposal 
    156           // accepted 
    157           Z2(k2, k1) = normp; 
    158           L_norm = L_normp; 
     154    // Unless this is the first, initialization iteration, update this column 
     155    // of independent random-walk normal variates to produce a multivariate 
     156    // i.i.d. normal proposal. It is this on which we build a multivariate 
     157    // log-volatility proposal for this multivariate sample. 
     158 
     159    // For each time series... 
     160    for(k2=0; k2<Nz[0]; k2++) { 
     161      // Initialize the scratchpad with the current proposal value and compute 
     162      // its log-likelihood 
     163      Y2(k2, 2) = Z2(k2, k1); 
     164      L0 = -0.5 * pow(Y2(k2, 2), 2); 
     165      // Proceed with the random walk only if this isn't the first, 
     166      // initialization iteration... 
     167      if(k0 > 0) { 
     168        // For each oversampling iteration... 
     169        for(k3=0; k3<n3; k3++) { 
     170          // Generate a uniform, symmetric, random walk proposal 
     171          normp = Y2(k2, 2) + wiggle * (drand48() - 0.5); 
     172          // Do the ratio in log space 
     173          L1 = -0.5 * pow(normp, 2); 
     174          dL = L1 - L0; 
     175          // The uniform random variate thus needs to be in logspace as well, 
     176          // but it's not always even needed. Thus the log operation is done 
     177          // just once, and only some of the time. 
     178          if(dL > 0  || log(drand48()) < dL) { 
     179            // Update current value and its log probability when proposal 
     180            // accepted 
     181            Y2(k2, 2) = normp; 
     182            L0 = L1; 
     183          } 
    159184        } 
    160185      } 
    161186    } 
    162  
     187     
    163188    // Now correlate the i.i.d. proposal to form a proposed column of 
    164189    // volatility shocks 
    165  
    166     // For each time series... 
    167     for(k2=0; k2<Nz[0]; k2++) { 
    168       // The matrix product of the cross-correlations and the independent 
    169       // variates... 
     190     
     191    // For each time series... 
     192    for(k2=0; k2<Nz[0]; k2++) { 
     193      // The matrix product of the cross-correlations and the proposed 
     194      // independent variates... 
    170195      sum = 0; 
    171196      for(km=0; km<Nz[0]; km++) 
    172         sum += R2(k2, km) * Z2(k2, k1); 
    173       // ...is the correlated volatility shock for this time series 
    174       V1(k2) = sum; 
    175     } 
    176  
    177     // Now do the VAR(1) to get a multivariate log-volatility proposal 
    178  
    179     // We alternate between the two columns of the log-volatility scratchpad 
    180     // for current & previous multivariate log-volatility values 
    181     kh = k1 % 2; 
    182  
    183     // For each time series... 
    184     for(k2=0; k2<Nv[0]; k2++) { 
    185       // The offset plus the current shock... 
    186       sum = D1(i) + V1(k2); 
    187       // ...plus the dot product for the VAR(1) term... 
    188       for(km=0; km<Nv[0]; km++) 
    189         sum += E2(k2, km) * H2(k2, 1-kh); 
    190       // ...is the modeled output 
    191       H2(k2, kh) = sum; 
    192     } 
    193  
    194     // Now inverse-scale the innovations for this sample by the 
    195     // log-volatilities in preparation for decorrelating them 
    196     for(k2=0; k2<Nv[0]; k2++) 
    197       Y2(k2, 0) = exp(-0.5 * H2(k2, kh)) * Z2(k2, k1); 
    198      
    199     // Now decorrelate the equi-variance innovations in preparation for 
    200     // (finally) the log-likelihood computation for this multivariate sample 
     197        sum += RVC2(k2, km) * Y2(k2, 2); 
     198      // ...is the correlated volatility shock proposal for this time series 
     199      Y2(k2, 3) = sum; 
     200    } 
     201     
     202    // Now do the VAR(1) to get a multivariate log-volatility proposal. The 
     203    // proposals are generated from "left to right" based solely on the 
     204    // random-walk normal proposals, so there's no need to worry about 
     205    // overwriting anything. 
     206 
     207    // For each time series... 
     208    for(k2=0; k2<Nz[0]; k2++) { 
     209      if(k1 == 0) { 
     210        // The log-volatility value of the first multivariate sample is just the 
     211        // offset scaled by its VAR gain. 
     212        H2(k2, 0) = D1(k2) / (1 - E2(k2, k2)); 
     213      } else { 
     214        // For all other samples, the offset plus the current proposed shock... 
     215        sum = D1(k2) + Y2(k2, 3); 
     216        // ...plus the dot product for the VAR(1) term... 
     217        for(km=0; km<Nz[0]; km++) 
     218          sum += E2(k2, km) * H2(k2, k1-1); 
     219        // ...is the proposed log-volatility value 
     220        H2(k2, k1) = sum; 
     221      } 
     222    } 
     223     
     224    // Now inverse-scale the innovations for this sample by the square root of 
     225    // the volatilities in preparation for decorrelating the innovations 
     226    for(k2=0; k2<Nz[0]; k2++) 
     227      Y2(k2, 4) = exp(-0.5 * H2(k2, k1)) * X2(k2, k1); 
     228     
     229    // Now decorrelate the equi-variance innovations in preparation for the 
     230    // log-likelihood computation for this multivariate log-volatility proposal 
    201231     
    202232    // For each time series... 
    …… 
    206236      sum = 0; 
    207237      for(km=0; km<Nz[0]; km++) 
    208         sum += S2(k2, km) * Y2(k2, 0); 
    209       // ...is the deccorrelated multivariate innovation for this sample 
    210       Y2(k2, 1) = sum; 
    211     } 
    212  
    213     // Finally, compute the log-likelihood of this deccorrelated multivariate 
    214     // innovation 
    215  
    216     // TODO 
    217      
     238        sum += RII2(k2, km) * Y2(k2, 4); 
     239      // ...is the decorrelated multivariate innovation for this sample, given 
     240      // its proposed log-volatility 
     241      Y2(k2, 5) = sum; 
     242    } 
     243 
     244    // Now compute log Pr(x[t]|hp[t]), the log-likelihood of the decorrelated 
     245    // multivariate innovation, conditional on the proposed log-volatility 
     246    // hp[t], which has been generated conditional on the log-volatility of the 
     247    // previous sample. Consider the proposal as a draw from Pr(hp[t]|h[t-1]). 
     248     
     249    sum = 0; 
     250    // For each time series... 
     251    for(k2=0; k2<Nz[0]; k2++) { 
     252      // Offset the decorrelated innovation by the mean, which is the 
     253      // correlation between the normal variates for innovation and proposed 
     254      // volatility shocks, times the value of the normal variate value for the 
     255      // proposed volatility shock for this sample 
     256      Y2(k2, 5) = Y2(k2, 5)  -  G1(k2) * H2(k2, k1); 
     257      // The log-likelihood is simply the argument of the exponential of the 
     258      // normal pdf, with the reciprocal scaling of the exponential being done 
     259      // in log space by subtraction. 
     260      sum -= pow(Y2(k2, 5) / Y2(k2, 0), 2) + Y2(k2, 1); 
     261      // Since the innovation has been inverse-scaled by the square root of the 
     262      // proposed volatility, the probability density needs to be scaled as 
     263      // well. In log space, this is a subtraction of half the proposed 
     264      // log-volatility. 
     265      sum -= 0.5 * H2(k2, k1); 
     266    } 
     267    P2(0, k1) = sum; 
     268     
     269    // The probability density that we need for our Metropolis-Hastings 
     270    // accept/reject decision is Pr(hp[t]|x[t],h[t+1]), where h[t+1] is the 
     271    // current log-volatility of the next sample. 
     272    // 
     273    // According to Bayes' Rule, the required density is proportional to 
     274    // Pr(x[t]|hp[t],h[t+1]) * Pr(hp[t]|h[t+1]). So, we need to compute and 
     275    // scale by the probability density of the log-volatility proposal given 
     276    // the current log-volatility of the next sample, which is, by Bayes rule: 
     277    // 
     278    // Pr(h[t+1]|hp[t]) * Pr(hp[t]) / Pr(h[t+1]) 
     279    //                    <-----+--------------> 
     280    //                          | 
     281    //                          +--- But this = 1 because both are equiprobable 
     282    // 
     283    // So, we now compute Pr(h[t+1]|hp[t]), which is Pr(h[t+1] - E2*hp[t]) 
     284    // unless we're at the last sample, in which case it is = 1. 
     285 
     286    if(k1 < Nz[1]) { 
     287      // First we need to remove the offset and the current proposal's VAR term 
     288      // from the next sample's current log-volatility 
     289 
     290      // For each time series... 
     291      for(k2=0; k2<Nz[0]; k2++) { 
     292        // ...the offset 
     293        sum = D1(k2); 
     294        // ...plus the matrix product of the Lag-1 cross-correlations and the 
     295        // proposed multivariate log-volatility components... 
     296        for(km=0; km<Nz[0]; km++) 
     297          sum += E2(k2, km) * H2(k2, k1); 
     298        // ...is the VAR term to subtract from the next sample's log-volatility 
     299        // value 
     300        Y2(k2, 6) = H2(k2, k1+1) - sum; 
     301      } 
     302 
     303      // Then we need to decorrelate the residuals to get independent normal 
     304      // variates 
     305 
     306      // For each time series... 
     307      for(k2=0; k2<Nz[0]; k2++) { 
     308        // The product of the inverted cross-correlation matrix and the 
     309        // correlated variates... 
     310        sum = 0; 
     311        for(km=0; km<Nz[0]; km++) 
     312          sum += RVI2(k2, km) * Y2(k2, 6); 
     313        // ...is the decorrelated multivariate normal that would generate the 
     314        // next sample's log-volatility given the one proposed for the current 
     315        // sample. 
     316        Y2(k2, 7) = sum; 
     317      } 
     318 
     319      // Then we need to compute the joint log-density of the decorrelated 
     320      // multivariate normal values 
     321      sum = - Nz[0] * 0.918938533205; // N*log(sqrt(2*pi)) 
     322      // For each time series... 
     323      for(k2=0; k2<Nz[0]; k2++) 
     324        sum -= 0.5 * pow(Y2(k2, 7), 2); 
     325        
     326      // Now do the scaling 
     327      // TODO 
     328 
     329    } 
     330 
     331 
     332           
     333        
     334       
     335       
     336 
     337     
     338 
     339    // Do a Metropolis-Hastings accept/reject based on the log-likelihood of 
     340    // the proposal, now in the "sum" variable, versus the log-likelihood of 
     341    // the previously accepted normal value for this sample. 
     342 
     343    // If this is the first, initialization iteration, we always accept the 
     344    // null "proposal," which was only obtained to (re)compute the 
     345    // log-likelihood given the current normal variates. 
     346    if(k0 == 0) { 
     347      P2(0, k1) = L0; 
     348      P2(1, k1) = L1; 
     349    } else { 
     350      // Do the ratio of likelihoods for hp[t] vs current h[t] in log space 
     351      dL = L1 - P2(0, k1); 
     352      // The uniform random variate thus needs to be in logspace as well, but 
     353      // it's not always even needed. Thus the log operation is done just once, 
     354      // and only some of the time. 
     355      if(dL > 0 || log(drand48()) < dL) { 
     356        P2(0, k1) = L0; 
     357        P2(1, k1) = L1; 
     358        // For each time series... 
     359        for(k2=0; k2<Nz[0]; k2++) 
     360          Z2(k2, k1) = Y2(k2, 2); 
     361      } 
     362    } 
     363 
     364    // End of Metropolis-Hastings step for this multivariate sample 
    218365  } 
     366 
     367  // End of this iteration 
    219368} 
     369 
     370// End of MCMC run. The array Z2 has been overwritten with updated 
     371// i.i.d. normal variates underlying an updated set of log-volatilities, and 
     372// the vector P1 has been overwritten with the log-likelihood of each 
     373// multivariate innovation sample given those log-volatilities. 
     374 
     375 

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

    r152 r153  
    219219    var = property(_get_var) 
    220220 
    221     def _get_nw(self): 
    222         if 'nw' not in self.cache: 
    223             self.cache['nw'] = sample.NormalWalk(self.p, self.n) 
    224         return self.cache['nw'] 
    225     nw = property(_get_nw) 
    226  
    227  
    228     #--- Supporting Methods --------------------------------------------------- 
    229  
    230     def draw_h(self, wiggle): 
    231         """ 
    232         Iterates my normal walkers by one step using the supplied fraction 
    233         I{wiggle} and returns a 2-D array of volatilities, one row for each of 
    234         my time series of observations. 
    235         """ 
    236         if not hasattr(self, 'h'): 
    237             self.h = s.empty_like(self.nw()) 
    238         self.nw.step(wiggle) 
    239         return self.inline('h', 'v', 'd', 'e', v=self.nw.correlate(self.f)) 
    240  
    241     def decorrelate(self, x, h=None): 
    242         """ 
    243         Call this method with a C{[pn]} array I{x} containing one row of C{n} 
    244         innovations for each of C{p} Wiener processes, with unit variance 
    245         assumed and correlated by a covariance matrix constructed from my 
    246         correlation parameters I{c}. The result will be the decorrelated 
    247         innovations of the Wiener processes. 
    248  
    249         You can supply a C{[pn]} array I{h} of log-volatilities for each of the 
    250         innovations. Then the innovations I{x} are assumed to be scaled by 
    251         their respective volatilities, and will be inversely scaled to meet the 
    252         assumption of unit variance. 
    253         """ 
    254         xs = x.shape 
    255         if xs[0] != self.p: 
    256             raise ValueError("Data array must have %d rows" % self.p) 
    257         hs = s.shape(h) 
    258         if hs: 
    259             if hs != xs: 
    260                 raise ValueError("Need one log-volatility per innovation") 
    261             x = s.exp(-0.5*h) * x.copy() 
    262         if 'decorrelate' not in self.cache: 
    263             self.cache['decorrelate'] = linalg.inv( 
    264                 linalg.cholesky(self.nw.covarMatrix(self.c), lower=True)) 
    265         r = self.cache['decorrelate'] 
    266         return self.inline('y', 'x', 'r', y=s.empty_like(x)) 
    267  
    268     def likelihoodOfInnovations(self, y): 
    269         """ 
    270         Call this method with a C{[pn]} array I{y} of uncorrelated 
    271         innovations. 
    272  
    273         The result will be the total log-likelihood of the innovations, given 
    274         my volatility/innovation shock correlations C{g} and the volatility 
    275         shocks currently in my random walker C{nw}. 
    276         """ 
    277         mu = +self.g * self.nw() 
    278         sigma2 = 1 - self.g**2 
    279         logProbs = -(y - mu)**2 / (2*sigma2) -  0.5 * s.log(2*s.pi*sigma2) 
    280         return logProbs.sum() 
     221 
     222    #--- Support -------------------------------------------------------------- 
    281223 
    282224    def covarMatrix(self, correlations): 
    …… 
    338280        x = self.var.reverse(self.a, self.b) 
    339281        # Concurrent cross-correlations between volatility shocks 
    340         r = linalg.cholesky(self.covarMatrix(self.f), lower=True)) 
     282        rvc = linalg.cholesky(self.covarMatrix(self.f), lower=True) 
     283        # Inverse of concurrent cross-correlations between volatility shocks 
     284        rii = linalg.inv(rvc) 
    341285        # Inverse of concurrent cross-correlations between innovation shocks 
    342         s = linalg.inv(linalg.cholesky(self.covarMatrix(self.c), lower=True)) 
    343         # Declare some scratchpad arrays 
    344         v = s.empty(self.p) 
    345         h = s.empty((self.p, 2)) 
    346         y = s.empty((self.p, 2)) 
     286        rii = linalg.inv(linalg.cholesky(self.covarMatrix(self.c), lower=True)) 
     287        # Declare working array with zeros 
     288        y = s.zeros((self.p, 8)) 
    347289        # Finally, run the hybrid Gibbs simulation of log-volatilities to get a 
    348290        # log-likelihood that is not (very) conditional on the latent 
    349291        # parameters that the volatilities represent 
     292        p = s.zeros((2, self.n)) 
     293        h = s.zeros(self.p, self.n)) 
    350294        L = self.inline( 
    351             'L', 
    352             'v', 'x', 'r', 's', 
    353             'z', 'h', 'y', 
     295            'p', 
     296            'z', 'x', 'rvc', 'rvi', 'rii', 'h', 'y', 
    354297            'd', 'e', 'g', 
    355             'n0', 'n3', 'wiggle', 
    356             L=0, wiggle=sigma/s.sqrt(self.n3)) 
    357         return L, v 
    358      
    359      
     298            'n0', 'n3', 'sigma')[1,:].sum() 
     299        return L, z, p 
     300     
     301    def foo(self): 
     302        pass 

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

    r151 r153  
    369369 
    370370 
    371 class Test_Model_other(Multivariate_BaseCase): 
    372     N = 100 
     371class Test_Model(Multivariate_BaseCase): 
     372    N = 10000 
    373373    corr = 0.0 
    374374    params = { 
    …… 
    377377        'c': [], 
    378378        'd': [0.01], 
    379         'e': [[0.8]], 
     379        'e': [[0.0]], 
    380380        'f': [], 
    381381        'g': [corr], 
    …… 
    400400        self.iv = random.multivariate_normal( 
    401401            [0, 0], [[1.0, self.corr], [self.corr, 1.0]], self.N).transpose() 
     402        #--- DEBUG, sinusoidal variates for log-volatility  -------------- 
     403        self.iv[1,:] = s.sin(s.linspace(0, 2*s.pi, self.N)) 
     404        #----------------------------------------------------------------- 
    402405        h = signaltools.lfilter( 
    403406            [1], [1.0, self.params['e'][0][0]], self.iv[1,:]) 
    404         self.x = s.exp(0.5*h) * self.iv[1,:] 
     407        self.x = s.exp(0.5*h) * self.iv[0,:] 
    405408        # Instantiate a model to be tested and set its parameters 
    406409        self.model = model.Model(tsList=[tseries.TimeSeries( 
    …… 
    408411        for name, value in self.params.iteritems(): 
    409412            setattr(self.model, name, value) 
    410      
    411     def test_likelihoodOfInnovations(self): 
    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 
    424  
    425     def test_logUniform(self): 
    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) 
    439413         
    440414    def test_likelihood(self): 
    441         N_plot = 100 
    442         wiggle = 0.001 
    443         self.model.M = 10000 
     415        N_plot = self.N 
     416        sigma = 0.05 
     417        self.model.n0 = 1000 
    444418        paramContainer = self.pcFactory() 
    445419        # 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() 
     420        paramContainer.latent = 0.0*self.iv[1,:] + 0.0*s.randn(1, self.N) 
     421        L, z, p = self.model.likelihood(paramContainer, sigma) 
    449422        if VERBOSE: 
    450423            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]) 
     424            vectors = (self.iv[1,:], self.x, p, z[0]) 
     425            for k, vector in enumerate(vectors): 
     426                sp = fig.add_subplot(100*len(vectors) + k + 11) 
     427                sp.plot(vector[:N_plot]) 
     428            #self._check_corr(self.iv[1,:], z[0]) 
    456429