Changeset 198

Show
Ignore:
Timestamp:
05/27/08 00:19:36 (6 months ago)
Author:
edsuom
Message:

Tried Metropolis-within-Gibbs, but ran into the difficulty of getting the proposal density on z for weighting

Files:

Legend:

Unmodified
Added
Removed
Modified
Copied
Moved

  • projects/AsynCluster/branches/metropolis-within-gibbs/doc/svpmc/example/svpmc.conf

    r196 r198  
    2727Variable    Value           Description 
    2828------------------------------------------------------------------------------- 
    29 D           6               Number of jump deviations in the D-kernel sim 
    30 Df          0.2             Fractional value of each dev from prev (or 1.0) 
     29D           4               Number of jump deviations in the D-kernel sim 
     30Ds          0.1             Starting jump deviation value 
     31Df          0.2             Fractional value of each dev from prev (or start) 
    3132wiggle      0.5             Range +/- of uniform random walks on z for h sim 
    32 N1          30              Proposals on z for each importance-sampling of h 
    33 N2          20              Hybrid-Gibbs iterations for h sim 
     33Ni          50              Hybrid-Gibbs rounds for simulation of h 
     34 
    3435 
    3536 
    …… 
    6768a           Distribution    Loc         Scale       a       b 
    6869------------------------------------------------------------------------------- 
    69 :           beta            0.00        0.005       2       5 
     70:           beta            0.00        0.002       2       3 
    7071 
    7172 

  • projects/AsynCluster/branches/metropolis-within-gibbs/svpmc/model.c

    r196 r198  
    104104 
    105105 
    106 // A proposal sample 
     106// An evaluation sample 
    107107struct sample 
    108108{ 
     109  // Log-likelihood of the time-series sample in the current evaluation, given 
     110  // h at this point 
     111  double Lxk; 
    109112  // Log-likelihood of the innovations up to the current one, given the history 
    110113  // of h to this point 
    111114  double Lx; 
    112   // Multivariate normal proposal, which is just z for time-series samples 
    113   // after the first in each log-volatility computation 
     115  // Multivariate normal proposal, which is just z for the current evaluation, 
     116  // and for all time-series samples after the first in each log-volatility 
     117  // computation 
    114118  double *z; 
    115   // Multivariate log-volatility value based on the current proposal, or the 
     119  // Multivariate log-volatility value based on the current evaluation, or the 
    116120  // previous proposal to use this struct if the current proposal hasn't been 
    117121  // crunched yet 
    …… 
    119123  // Scratchpad 1-D arrays of the same length as z, h 
    120124  double *x0, *x1; 
    121 }; 
    122  
    123  
    124 // A selection of a proposed log-volatility simulation 
    125 struct selection 
    126 { 
    127   int kp; 
    128   double Lp; 
    129125}; 
    130126 
    …… 
    212208  // Inverse-scale the innovation by the square root of the volatilities in 
    213209  // preparation for decorrelation 
     210 
    214211  for(k=0; k<N; k++) 
    215     // TODO: Should we use a stripped-down version of exp here for speed? 
     212    // We could use a stripped-down version of exp here for speed; doing it 
     213    // with the stdlib exp function takes about 70 CPU cycles (20 nS) on an 
     214    // Intel QX9760 running at 3.8 GHz, though some of that is the unavoidable 
     215    // array lookups and products in the equation. 
    216216    SPA1(S, x0, k) = exp(-0.5 * SPA1(S, h, k)) * PA1(x, k); 
    217      
     217 
    218218  // Now decorrelate the equi-variance innovations in preparation for the 
    219219  // log-likelihood computation for this multivariate log-volatility proposal 
    …… 
    223223  // multivariate innovation for this sample, conditional on the sample's 
    224224  // log-volatility. 
    225    
    226   // Now the log-likelihoood. For each time series... 
     225 
     226  S->Lxk = 0; 
     227 
     228  // For each time series... 
    227229  for(k=0; k<N; k++) { 
    228230    // Offset the decorrelated innovation by the mean, which is the correlation 
    …… 
    237239    val = SPA1(S, x0, k) / SPP2(P, sc, k, 0); 
    238240    // Fast squaring instead of pow() 
    239     S->Lx -= 0.5 * val*val + SPP2(P, sc, k, 1); 
     241    S->Lxk -= 0.5 * val*val + SPP2(P, sc, k, 1); 
    240242    // Since the innovation has been inverse-scaled by the square root of the 
    241243    // proposed volatility, the probability density needs to be scaled as 
    242244    // well. In log space, this is a subtraction of half the log-volatility. 
    243     S->Lx -= 0.5 * SPA1(S, h, k); 
    244   } 
    245 } 
    246  
    247  
    248 // Importance-samples one of the proposals based on a weight computed from the 
    249 // log-density of z and the log-likelihood of the innovation given the 
    250 // proposal's log-volatility. 
    251 // 
    252 // Returns a struct with an integer index to the sampled proposal's struct and 
    253 // the log-probability of the proposal's having been selected. 
    254 struct selection \ 
    255 importanceSample(struct sample sp[], double Lz[], int n) 
    256 { 
    257   int k; 
    258   double val; 
    259   double Dp_sum = 0; 
    260   double Lp[n], Dp[n]; 
    261   double L_max = -HUGE_VAL; // Anything will be more positive 
    262   struct selection result; 
    263  
    264   // Normalize log-densities to prevent overflow 
    265   for(k=0; k<n; k++) { 
    266     val = sp[k].Lx + Lz[k]; 
    267     if(val > L_max) 
    268       L_max = val; 
    269     Lp[k] = val; 
    270   } 
    271  
    272   // Compute linear probability density for each proposal, relative to the 
    273   // maximum one 
    274   for(k=0; k<n; k++) { 
    275     val = exp(Lp[k] - L_max); 
    276     Dp[k] = val; 
    277     Dp_sum += val; 
    278   } 
    279  
    280   // Accept-reject sampling with uniform random proposal selection 
    281   do { 
    282     // Randomly select one of the proposals... 
    283     result.kp = (int) uniform(0, n); 
    284     // ...until the selected proposal is accepted, with probability 
    285     // proportional to its relative weight 
    286   } while(uniform(0, 1.0) > Dp[result.kp]); 
    287  
    288   // The probability density of the selected log-volatility proposal 
    289   result.Lp = log(Dp[result.kp] / Dp_sum); 
    290   return result; 
     245    S->Lxk -= 0.5 * SPA1(S, h, k); 
     246  } 
     247  // Add the current log-likelihood to the accumulated log-likelihoods 
     248  S->Lx += S->Lxk; 
    291249} 
    292250 
    …… 
    304262// w    Wiggle value for +/- proposals (float) 
    305263// ni   Number of hybrid-Gibbs iterations (int) 
    306 // np   Number of proposals to try per log-volatility simulation (int) 
    307 // ne   Max number of time-series samples to evaluate per proposal (int) 
    308264 
    309265//--- Model parameters, direct --- 
    …… 
    322278 
    323279 
     280// Run evaluations until odds of error in selection between current and 
     281// proposed z will be less than 1/100. 
     282#define MIN_DIFF 0.01 
     283 
    324284// The number of time series making up a single multivariate sample 
    325285const int Nt = Nx[0]; 
    …… 
    327287const int Ns = Nx[1]; 
    328288 
    329 int ki, k0, k1, k2, k3; 
     289int ki, ke; 
     290int k0, k1, k2, k3; 
    330291double val; 
    331292double Lx = 0; 
    …… 
    334295// Multivariate innovation for one current time-series sample 
    335296double *xk; 
    336 // The value of z and h for each proposal, at the first time-series sample of 
    337 // the evaluation interval 
    338 double zp[np], hp[np][Nt]; 
    339 // Log probability density of the value of z for each proposal 
    340 double Lz[np]; 
    341 // Log-likelihood of the innovation for the first time-series sample in 
    342 // proposal evalutions, for each proposal 
    343 double Lxk[np]; 
    344  
    345 // A struct for the result of the importance-sampling function 
    346 struct selection spSelection; 
    347 // A set of proposal sample structs 
    348 struct sample sp[np]; 
     297// The multivariate proposal on z 
     298double zp[Nt]; 
     299// The multivariate value of h for current and proposal evaluations, at the 
     300// first time-series sample of the evaluation interval 
     301double he[2][Nt]; 
     302// Log probability density of the value of z for current and proposal 
     303// evaluations 
     304double Lz[2]; 
     305// Log-likelihood of the innovation for the first time-series sample in current 
     306// and proposal evalutions 
     307double Lxk[2]; 
     308 
     309// A set of evaluation sample structs 
     310struct sample se[2]; 
    349311// A struct for the model parameters 
    350312struct parameters P; 
    …… 
    359321// Initialize various arrays 
    360322xk = (double *) malloc(sizeof(double) * Nt); 
    361 for(k0=0; k0<np; k0++) { 
    362   sp[k0].z = (double *) malloc(sizeof(double) * Nt); 
    363   sp[k0].h = (double *) malloc(sizeof(double) * Nt); 
    364   sp[k0].x0 = (double *) malloc(sizeof(double) * Nt); 
    365   sp[k0].x1 = (double *) malloc(sizeof(double) * Nt); 
     323for(k0=0; k0<2; k0++) { 
     324  se[k0].z = (double *) malloc(sizeof(double) * Nt); 
     325  se[k0].h = (double *) malloc(sizeof(double) * Nt); 
     326  se[k0].x0 = (double *) malloc(sizeof(double) * Nt); 
     327  se[k0].x1 = (double *) malloc(sizeof(double) * Nt); 
    366328} 
    367329 
    …… 
    373335  // the log-volatility offset plus a simulated, latent-parameter value. 
    374336  // TODO 
    375   for(k0=0; k0<np; k0++) { 
     337  for(k0=0; k0<2; k0++) { 
    376338    for(k1=0; k1<Nt; k1++) 
    377       PA1(sp[k0].h, k1) = D1(k1); 
     339      PA1(se[k0].h, k1) = D1(k1); 
    378340  } 
    379341 
    …… 
    381343  for(k0=0; k0<Ns; k0++) { 
    382344     
    383     // Simulate np sets of h samples from separate proposals on z from the 
    384     // current time-series sample to some subsequent time-series sample at which 
    385     // the log-likelihood of the innovation for that sample becomes substantially 
    386     // the same for all proposals. 
     345    // Simulate two sets of h samples from the current z and a proposal on z, 
     346    // with a likelihood evaluation that starts from the current time-series 
     347    // sample to some subsequent time-series sample at which the log-likelihood 
     348    // of the innovation for that sample becomes substantially the same for the 
     349    // current and proposal value of z. 
    387350   
    388351    k1 = k0; 
    389352   
    390     // Clear the log-likelihood accumulator of each proposal struct for the next 
    391     // proposal evalution 
    392     for(k2=0; k2<np; k2++) 
    393       sp[k2].Lx = 0; 
    394    
     353    // Clear the log-likelihood and normal log-density accumulators of each 
     354    // evaluation struct for the next pair of evaluations. 
     355    for(k2=0; k2<2; k2++) { 
     356      Lz[k2] = 0; 
     357      se[k2].Lx = 0; 
     358    } 
     359     
    395360    do { 
    396361   
    …… 
    399364        PA1(xk, k2) = X2(k2, k1); 
    400365       
    401       // For each proposal... 
    402       for(k2=0; k2<np; k2++) { 
     366      // For each evaluation... 
     367      for(k2=0; k2<2; k2++) { 
    403368   
    404369        // Set z... 
    …… 
    406371          val = Z2(k3, k1); 
    407372          if(k1==k0) { 
    408             // First time-series sample, use a proposal on z instead of z itself 
    409             val += uniform(-w, +w); 
    410             zp[k2] = val; 
    411             Lz[k2] = normLogDensity(val); 
     373            if(k2==1) { 
     374              // First time-series sample for the second evaluation, use a 
     375              // proposal on z instead of z itself 
     376              val += uniform(-w, +w); 
     377              zp[k3] = val; 
     378            } 
     379            // Store the log-density of z for the first time-series sample to 
     380            // compare current vs. proposal 
     381            Lz[k2] += normLogDensity(val); 
    412382          } 
    413           PA1(sp[k2].z, k3) = val; 
     383          PA1(se[k2].z, k3) = val; 
    414384        } 
    415385        // ...compute the multivariate log-volatility 
    416         logVol(&sp[k2], &P); 
     386        logVol(&se[k2], &P); 
    417387        // ...and the log-likelihood of the innovation for this time-series 
    418388        // sample, given the log-volatility 
    419         logLikelihood(&sp[k2], &P, xk); 
     389        logLikelihood(&se[k2], &P, xk); 
    420390         
    421391        // Save the first log-volatility and the log-likelihood based on it (and 
    …… 
    424394        if(k1==k0) { 
    425395          for(k3=0; k3<Nt; k3++) 
    426             hp[k2][k3] = PA1(sp[k2].h, k3); 
    427           Lxk[k2] = sp[k2].Lx; 
     396            he[k2][k3] = PA1(se[k2].h, k3); 
     397          Lxk[k2] = se[k2].Lx; 
    428398        } 
    429399      } 
    …… 
    431401      // Ready for the next time series sample, assuming there is one 
    432402      k1++; 
     403 
     404      // Continue only if difference in log-likelihood of the current 
     405      // time-series sample of the evaluations between current and proposed z 
     406      // is great enough to warrant doing so. 
     407      val = se[0].Lxk - se[1].Lxk; 
    433408     
    434     } while (k1<Ns && k1<k0+ne); 
    435      
    436     // Importance sample proposal using pProb=exp(Lz+Lx) for each... 
    437     spSelection = importanceSample(sp, Lz, np); 
    438     // ...add the log-likelihood of the innovation given the selected 
    439     // log-volatility proposal to what will be the total log-likelihood of the 
    440     // innovations given the entire simulated h... 
    441     Lx += Lxk[spSelection.kp]; 
    442     // ...add the log-density of the selected log-volatility proposal to what 
    443     // will be the total log-density of the entire simulated h... 
    444     Lh += spSelection.Lp; 
    445     // ...and update the normal variate underlying the current time-series sample 
    446     // to that on which the selected log-volatility proposal was based. 
    447     for(k1=0; k1<Nt; k1++) 
    448       Z2(k1, k0) = zp[spSelection.kp]; 
    449    
    450     // Set the initial log-volatility of the proposals for the next time-series 
    451     // sample to the log-volatility of the selected proposal for this time-series 
    452     // sample. 
     409    } while (k1<Ns && (val > MIN_DIFF || val < -MIN_DIFF)); 
     410 
     411    // Metropolis-hastings selection of current or proposal using Lz+Lx for 
     412    // each evaluation 
     413    val = Lz[1] + se[1].Lx  -  Lz[0] - se[0].Lx; 
     414    if(val > 0) { 
     415      // Proposal accepted with probability p=1 
     416      ke = 1; 
     417    } else if(log(uniform(0, 1)) < val) { 
     418      // Proposal accepted with probability Pr(zp) / Pr(z) 
     419      ke = 1; 
     420      Lh += val; 
     421    } else { 
     422      // Proposal rejected 
     423      ke = 0; 
     424    } 
     425 
     426    // If the proposal was accepted, overwrite z for this sample with the 
     427    // proposal on z 
     428    if(ke == 1) { 
     429      for(k1=0; k1<Nt; k1++) 
     430        Z2(k1, k0) = zp[k1]; 
     431    } 
     432 
     433    // Set the initial log-volatility of the evaluations for the next 
     434    // time-series sample to the log-volatility of the selected evaluation for 
     435    // this time-series sample. 
    453436    for(k1=0; k1<Nt; k1++) { 
    454       val = hp[spSelection.kp][k1]; 
    455       for(k2=0; k2<np; k2++) 
    456         PA1(sp[k2].h, k1) = val; 
     437      val = he[ke][k1]; 
     438      for(k2=0; k2<2; k2++) 
     439        PA1(se[k2].h, k1) = val; 
    457440    } 
     441 
     442    // Accumulate the log-likelihood of the innovation given the selected 
     443    // log-volatility evaluation 
     444    Lx += Lxk[ke]; 
    458445   
    459446    // Done with log-volatility draw for this time-series sample 
    …… 
    466453// sample in the last iteration 
    467454for(k0=0; k0<Nt; k0++) 
    468   H1(k0) = PA1(sp[spSelection.kp].h, k0); 
     455  H1(k0) = PA1(se[ke].h, k0); 
    469456 
    470457// Free array memory 
    471458free(xk); 
    472 for(k0=0; k0<np; k0++) { 
    473     free(sp[k0].z); 
    474     free(sp[k0].h); 
    475     free(sp[k0].x0); 
    476     free(sp[k0].x1); 
    477 } 
    478  
    479 // The pair of results is Lx, the log-likelihood of each innovation given its 
    480 // simulated log-volatility and Lh, the total log-density of all of the 
    481 // simulated log-volatilities 
    482 result[0] = Lx; 
     459for(k0=0; k0<2; k0++) { 
     460    free(se[k0].z); 
     461    free(se[k0].h); 
     462    free(se[k0].x0); 
     463    free(se[k0].x1); 
     464} 
     465 
     466// The pair of results is the integrated log-likelihood of the innovations 
     467// given the log-volatilities simulated over the iterations... 
     468result[0] = Lx / ni; 
     469// ...and the integrated log-density of the log-volatilities simulated over the 
     470// iterations, conditional on the starting z 
    483471result[1] = Lh; 
    484472return_val = result; 

  • projects/AsynCluster/branches/metropolis-within-gibbs/svpmc/model.py

    r196 r198  
    184184      generator. 
    185185 
    186     @ivar N1: Number of log-volatility proposals to compute for the importance 
    187       sampling of each hybrid Gibbs step. For univariate, N1=10 appears most 
    188       cost-effective with wiggle=2.0. Scale for smaller wiggle and to account 
    189       for multivariate. 
    190  
    191     @ivar N2: Number of hybrid-Gibbs iterations per likelihood evaluation. For 
    192       univariate, N2=2 appears most cost-effective, but set to 4 anyays. Maybe 
    193       exploration of multivariate space needs more iterations. 
    194  
    195     @ivar Ne: The number of successive time-series samples to include in the 
    196       evaluation of each log-volatility proposal. For some weird reason, Ne=3 
    197       is optimal, at least for univariate. 
    198      
    199     """ 
    200     keyAttrs = {'y':None, 'wiggle':1.0, 'N1': 30, 'N2':10, 'Ne':3} 
     186    @ivar Ni: Number of hybrid Gibbs iterations per likelihood evaluation. 
     187     
     188    """ 
     189    keyAttrs = {'y':None, 'wiggle':1.0, 'Ni':30} 
    201190 
    202191    #--- Properties ----------------------------------------------------------- 
    …… 
    264253        If a L{linalg.LinAlgError} is raised due to an invalid correlation 
    265254        matrix parameter, the method returns with no result. Otherwise, it 
    266         returns a tuple containing the following values reached at the end of 
    267         the log-volatility simulation: 
     255        returns a tuple containing the following: 
    268256         
    269257            1. The log-likelihood of my observations given the simulated 
    270                log-volatility values. 
    271  
    272             2. The log-probability of the simulated log-volatilities. 
     258               log-volatility values, integrated over the simulation rounds. 
     259 
     260            2. The log-probability of the simulated log-volatilities, 
     261               integrated over the simulation rounds. 
    273262 
    274263            3. The IID normal variates underlying the simulated 
    275                log-volatilities. 
     264               log-volatilities after the last simulation round. 
    276265 
    277266            4. The multivariate log-volatility value for the last time-series 
    278                sample, useful for extrapolating from the fitted model. 
     267               sample after the last simulation round. (Useful for 
     268               extrapolating from the fitted model.) 
    279269         
    280270        The returned likelihood does not consider the prior probability of the 
    …… 
    309299        sc[:,1] = s.log(s.sqrt(2*s.pi) * sc[:,0]) 
    310300 
    311         print self.wiggle, self.N1, self.N2, self.Ne 
    312         # Run the hybrid Gibbs rounds, returning the likelihood from the last 
    313         # one along with the state of the IID variate vector at that point. 
     301        # Run the hybrid Gibbs rounds 
    314302        Lx, Lh = self.inline( 
    315             'z', 'h', 'x', 
    316             'w', 'ni', 'np', 'ne', 
     303            'z', 'h', 'x', 'w', 'ni', 
    317304            'd', 'e', 'g', 'rz', 'ri', 'sc', 
    318             w=self.wiggle, np=self.N1, ni=self.N2, ne=self.Ne) 
    319         return Lx/self.N2, Lh/(self.N2*self.Ne), z, h 
     305            w=self.wiggle, ni=self.Ni) 
     306        return Lx, Lh, z, h 

  • projects/AsynCluster/branches/metropolis-within-gibbs/svpmc/project.py

    r196 r198  
    102102                value = int(value) 
    103103            setattr(self, name, value) 
    104         V = [1.0] 
     104        V = [self.Ds] 
    105105        for k in xrange(self.D-1): 
    106106            V.append(self.Df * V[-1]) 
    …… 
    112112        paramTitles, dimensions = self._setupParams() 
    113113        self.mgr = model.ModelManager( 
    114             self, tsData, wiggle=self.wiggle, N1=self.N1, N2=self.N2) 
     114            self, tsData, wiggle=self.wiggle, Ni=self.Ni) 
    115115        self._setupCDF( 
    116116            os.path.join(specDir, ncFileName), 

  • projects/AsynCluster/branches/metropolis-within-gibbs/svpmc/test/test_model.py

    r192 r198  
    506506 
    507507 
    508 class Test_Model_importanceSample(BaseCase): 
    509     code = """ 
    510     int k; 
    511     const int n = NX[0]; 
    512     struct sample sp[n]; 
    513     double Lz[n]; 
    514     py::tuple result(2); 
    515     struct selection spSelection; 
    516  
    517     for(k=0; k<n; k++) { 
    518       sp[k].Lx = X1(k); 
    519       Lz[k] = Z1(k); 
    520     } 
    521     spSelection = importanceSample(sp, Lz, n); 
    522     result[0] = spSelection.kp; 
    523     result[1] = spSelection.Lp; 
    524     return_val = result; 
    525     """ 
    526  
    527     def _importanceSample(self, Lx_array, Lz_array=None): 
    528         if Lz_array is None: 
    529             Lz_array = s.zeros_like(Lx_array) 
    530         return self.inline(self.code, X=Lx_array, Z=Lz_array) 
    531  
    532     def test_basic(self): 
    533         N = 8 
    534         N_iter = 1000 
    535         Lx_array = s.zeros(N) 
    536         I = [self._importanceSample(Lx_array)[0] for k in xrange(N_iter)] 
    537         for k in xrange(N): 
    538             self.failUnless(k in I) 
    539  
    540     def _check(self, z0, wiggle): 
    541         N = 20 
    542         N_iter = 1000 
    543         distObj = stats.distributions.truncnorm(z0-wiggle, z0+wiggle) 
    544         X = 2*wiggle*(s.rand(N_iter, N) - 0.5) + z0 
    545         Y = s.empty(N_iter) 
    546         Lp = s.empty(N_iter) 
    547         for k, Xk in enumerate(X): 
    548             Lx_array = s.log(distObj.pdf(Xk)) 
    549             kp, Lpk = self._importanceSample(Lx_array) 
    550             Y[k] = Xk[kp] 
    551             Lp[k] = Lpk 
    552         if VERBOSE: 
    553             I = Y.argsort() 
    554             sp = self.fig.add_subplot(111) 
    555             sp.plot(Y, s.exp(Lp), '.', Y[I], s.pi/N*distObj.pdf(Y[I])) 
    556         self.failUnlessNormal(Y) 
    557  
    558     def test_full(self): 
    559         return self._check(0, 4.0) 
    560  
    561     def test_closePortion(self): 
    562         return self._check(1.5, 1.0) 
    563  
    564     def test_farPortion(self): 
    565         return self._check(2.5, 0.5) 
    566  
    567  
    568508class Test_Model_likelihood(BaseCase): 
    569509    corr = 0.0 
    …… 
    602542        return iv, h, x 
    603543 
    604     def _inlineVars(self, wiggle, eValue, N1, N2=1, Ne=3): 
     544    def _inlineVars(self, wiggle, eValue, Ni): 
    605545        return { 
    606546            'h':s.empty(1), 
    607547            'sc':s.array([1, s.log(s.sqrt(2*s.pi))]), 
    608548            'w':wiggle, 
    609             'np':N1, 
    610             'ni':N2, 
    611             'ne':Ne, 
     549            'ni':Ni, 
    612550            'd':s.array([0.0]), 
    613551            'e':s.array([[float(eValue)]]), 
    …… 
    617555            } 
    618556     
    619     def _modelAndParamFactory(self, wiggle, N1, N2, **kw): 
     557    def _modelAndParamFactory(self, wiggle, Ni, **kw): 
    620558        """ 
    621559        Returns a parameter container populated from my I{params} dict, 
    …… 
    633571                s.array(kw.get(name, self.params[name]), dtype='d')) 
    634572        # Return the parameter container along with a model to be tested 
    635         modelObj = model.Model(N1=N1, N2=N2, wiggle=wiggle) 
     573        modelObj = model.Model(Ni=Ni, wiggle=wiggle) 
    636574        modelObj.paramNames = util.PARAM_NAMES 
    637575        return pc, modelObj 
    …… 
    651589        """ 
    652590        Runs C code for L{model.Model.likelihood} independent of the method 
    653         itself. 
     591        itself, repeating the run I{Nr} times with different simulated data. 
    654592        """ 
    655593        z = s.zeros((1, Ns)) 
    …… 
    667605        wiggle = 0.1 
    668606        eValue = 0.90 
    669         Ns, N1, N2 = 10, 5, 20 
    670         inlineVars = self._inlineVars(wiggle, eValue, N1, N2) 
    671         inlineVars.update({'x':2*s.ones(Ns), 'z':s.zeros((1, Ns))}) 
     607        Ns, Nr = 6, 50 
     608        inlineVars = self._inlineVars(wiggle, eValue, 1) 
     609        inlineVars.update({'x':2*s.ones((1,Ns)), 'z':s.zeros((1, Ns))}) 
    672610        modelObj = model.Model() 
    673611        print "\n  k      Lx       Lh     Simulated Volatility Shocks" 
    674612        print "-"*100 
    675         for k in xrange(N2): 
     613        for k in xrange(Nr): 
    676614            Lx, Lh = modelObj.inlineDirect( 
    677615                modelObj.code['likelihood'], **inlineVars) 
    678616            zLine = "   ".join(["%+4.2f" % xx for xx in inlineVars['z'][0,:]]) 
    679617            print "%3d: %+6.1f  %+7.1f%s%s" % (k, Lx, Lh, " "*5, zLine) 
    680  
     618     
    681619    def _simRounds(self, Nr, Ns, inlineVars, name, values): 
    682620        fig = self.fig 
    …… 
    694632            sp.set_title("%s=%d" % (name, value)) 
    695633 
    696     def test_optimize_Ne(self): 
    697         # Empirically, it looks like Ne=3 is best, just slightly better than 
    698         # Ne=2 and Ne=4. Beyond Ne=4 things decline significantly. 
     634 
     635    def test_optimize_Ni(self): 
     636        """ 
     637        With univariate, Ni=20 is significantly better than Ni=10, but not 
     638        significantly worse than Ni=40. 
     639        """ 
    699640        Nr = 1000 
    700         wiggle = 3.0 
     641        wiggle = 1.0 
    701642        eValue = 0.95 
    702         Ns, N1, N2 = 1000, 10, 2 
    703         inlineVars = self._inlineVars(wiggle, eValue, N1, N2) 
    704         self._simRounds(Nr, Ns, inlineVars, 'ne', (2,3,4,5,6)) 
    705  
    706     def test_optimize_N1(self): 
    707         # Empirically, it looks like np=10 is best. Going to np=20 barely gives 
    708         # any improvement, and obviously doubles the amount of computation 
    709         # time. Dropping to np=5 is significantly worse. 
    710         Nr = 1000 
    711         wiggle = 3.0 
    712         eValue = 0.95 
    713         Ns, N2 = 1000, 2 
    714         inlineVars = self._inlineVars(wiggle, eValue, 1, N2) 
    715         self._simRounds(Nr, Ns, inlineVars, 'np', (2,5,10,20)) 
    716  
    717     def test_optimize_N2(self): 
    718         # Empirically, it looks like ni=2 is best. It's noticeably better than 
    719         # ni=1 (though not dramatically so), but going to anything larger 
    720         # doesn't show any improvement. 
    721         Nr = 1000 
    722         wiggle = 2.0 
    723         eValue = 0.95 
    724         Ns, N1 = 1000, 10 
    725         inlineVars = self._inlineVars(wiggle, eValue, N1, 1) 
    726         self._simRounds(Nr, Ns, inlineVars, 'ni', (1,2,3,4)) 
     643        Ns = 1000 
     644        inlineVars = self._inlineVars(wiggle, eValue, 1) 
     645        self._simRounds(Nr, Ns, inlineVars, 'ni', (5,10,20,40)) 
    727646     
    728647    def test_highCorrelation_c(self): 
    729         Ne = 3 
    730         wiggle = 2.0 
     648        wiggle = 1.0 
    731649        eValue = 0.97 
    732         Ns, N1, N2 = 200, 10, 2 
    733         inlineVars = self._inlineVars(wiggle, eValue, N1, N2, Ne) 
     650        Ns, Ni = 200, 20 
     651        inlineVars = self._inlineVars(wiggle, eValue, Ni) 
    734652        for hReal, hSim in self._likelihood(Ns, 1, inlineVars): 
    735653            self._plot((hReal, hSim)) 
    736654     
    737655    def test_highCorrelation_model(self): 
    738         wiggle = 2.0 
    739         eValue = 0.95 
    740         Ns, N1, N2 = 1000, 10, 2 
     656        wiggle = 1.0 
     657        eValue = 0.97 
     658        Ns, Ni = 1000, 20 
    741659        iv, h, x = self._simData(Ns, eValue) 
    742         pc, modelObj = self._modelAndParamFactory(wiggle, N1, N2, e=[[eValue]]) 
     660        pc, modelObj = self._modelAndParamFactory(wiggle, Ni, e=[[eValue]]) 
    743661        pc.z = s.zeros((1, Ns)) 
    744662        modelObj.y = self._x_to_y(x) 

  • projects/AsynCluster/branches/metropolis-within-gibbs/svpmc/weave.py

    r179 r198  
    3030 
    3131class my_info(custom_info): 
    32     _extra_compile_args = ['-w'] 
     32    _extra_compile_args = ['-w', '-O3', '-march=native'] 
    3333    __macros = [ 
    3434        # Access an element i of a 1-D array whose pointer is pa 

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

    r179 r198  
    3030 
    3131class my_info(custom_info): 
    32     _extra_compile_args = ['-w'] 
     32    _extra_compile_args = ['-w', '-O3', '-march=native'] 
    3333    __macros = [ 
    3434        # Access an element i of a 1-D array whose pointer is pa