Changeset 182

Show
Ignore:
Timestamp:
05/15/08 13:35:24 (7 months ago)
Author:
edsuom
Message:

Grinding along with C code for new model

Files:

Legend:

Unmodified
Added
Removed
Modified
Copied
Moved

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

    r181 r182  
    111111  // Log-likelihood of the current innovation, given h 
    112112  double Lx; 
     113  // Linear probability density of the proposal, given h 
     114  double pProb; 
    113115  // Multivariate normal proposal 
    114116  double *z; 
    …… 
    240242 
    241243 
    242 // Returns with zero only if the proposals are close enough to finish 
    243 // likelihood computation 
    244 int notCloseEnough(struct sample *sp[], int n) 
     244// Returns with zero only if the post-first proposals are close enough to 
     245// finish likelihood computation 
     246int notCloseEnough(struct sample *sp[2][], int n) 
    245247{ 
    246248  int k; 
    …… 
    249251  double Lmax = -1E20;   // Anything will be more positive 
    250252  for(k=0; k<n; k++) { 
    251     L = sp[k]->Lx; 
     253    L = sp[1][k]->Lx; 
    252254    if(L < Lmin) 
    253255      Lmin = L; 
    …… 
    270272// x    Innovation values, [pn] array 
    271273// w    Wiggle value for +/- proposals (float) 
    272 // n    Number of proposals to try per sample (int) 
     274// n    Number of proposals to try per time-series sample (int) 
    273275 
    274276//--- Model parameters, direct --- 
    …… 
    289291const int Ns = Nz[1]; 
    290292 
    291 int k0, k1, k2, k3; 
    292 double zp, Lx, Lh; 
    293  
    294 struct sample sp[n]; 
     293int notFirst; 
     294int k0, k1, k2, k3, kp; 
     295double zp; 
     296double Lp = 0; 
     297 
     298double *xk; 
     299struct sample sp[2][n]; 
    295300struct parameters P; 
    296301py::tuple result(2); 
    …… 
    302307 
    303308// Initialize various arrays 
    304 x0 = (double *) calloc(sizeof(double), Nt); 
     309xk = (double *) calloc(sizeof(double), Nt); 
    305310for(k0=0; k0<n; k0++) { 
    306   sp[k0].h = (double *) calloc(sizeof(double), Nt); 
    307   sp[k0].x0 = (double *) calloc(sizeof(double), Nt); 
    308   sp[k0].x1 = (double *) calloc(sizeof(double), Nt); 
     311  for(k1=0; k1<n; k1++) { 
     312    sp[k0][k1].h = (double *) calloc(sizeof(double), Nt); 
     313    sp[k0][k1].x0 = (double *) calloc(sizeof(double), Nt); 
     314    sp[k0][k1].x1 = (double *) calloc(sizeof(double), Nt); 
     315  } 
    309316} 
    310317 
    …… 
    314321  // Simulate n sets of h samples from separate proposals on z until a 
    315322  // substantially common-valued sample is reached 
     323 
    316324  k1 = k0; 
    317325  do { 
    318     // Generate the next set of sample structs for the proposals... 
    319     for(k2=0; k2<n, k2++) { 
     326 
     327    // Prepare current-innovation vector xk and a set of sample structs for 
     328    // proposals, one set for the first time-series sample and another set for 
     329    // any remaining time-series samples to be computed 
     330    for(k2=0; k2<Nt; k2++) 
     331      PA1(xk, k2) = X2(k2, k1); 
     332    if(k1 == k0) { 
     333      // First time-series sample 
     334      notFirst = 0; 
     335    } else if(k1 == k0+1) { 
     336      // Second time-series sample, copy h arrays from proposal structs for 
     337      // first time-series sample 
     338      notFirst = 1; 
     339      for(k2=0; k2<n; k2++) { 
     340        for(k3=0; k3<Nt; k3++) 
     341          SPA1(sp[1][k2], h, k3) = SPA1(sp[0][k2], h, k3); 
     342      } 
     343    } else { 
     344      // Third or later time-series sample 
     345      notFirst = 1; 
     346    } 
     347 
     348    // Make proposals on z in the set of sample structs... 
     349    for(k2=0; k2<n; k2++) { 
    320350      for(k3=0; k3<Nt; k3++) { 
    321351        zp = Z2(k3, k1) + uniform(-w, +w); 
    322         SPA1(sp[k2], z, k3) = zp; 
    323         sp[k2].Lz = normLogDensity(zp); 
     352        SPA1(sp[notFirst][k2], z, k3) = zp; 
     353        sp[notFirst][k2].Lz = normLogDensity(zp); 
    324354      } 
    325355    } 
    326356    // ...compute their log-volatilities 
    327     for(k2=0; k2<n, k2++) 
    328       logVol(&sp[k2], &P); 
    329     // ...and the log-likelihood of the innovation for this sample, for each 
    330     // proposal on z and its log-volatility 
    331     for(k2=0; k2<Nt; k2++) 
    332       PA1(x0, k2) = X2(k2, k1); 
    333357    for(k2=0; k2<n; k2++) 
    334       logLikelihood(&sp[k2], &P, x0); 
    335  
    336     // TODO 
     358      logVol(&sp[notFirst][k2], &P); 
     359    // ...and the log-likelihood of the innovation for this time-series sample, 
     360    // for each proposal on z and its log-volatility 
     361    for(k2=0; k2<n; k2++) 
     362      logLikelihood(&sp[notFirst][k2], &P, xk); 
    337363 
    338364    // Ready for the next time series sample, if we go that far 
    339365    k1++; 
    340366 
    341   } while (notCloseEnough(sp, n)); 
    342  
    343   // TODO: Importance sample proposal using total Lz+Lh 
    344  
    345 } 
     367  } while (k1==k0+1 || (k1<Ns && notCloseEnough(sp, n))); 
     368 
     369  // Importance sample proposal using pProb=exp(Lz+Lx) for each 
     370  kp = sample(sp, n); 
     371  Lp += log(sp[1][kp]->pProb); 
     372 
     373} 
     374 
     375// Save the multivariate log-volatility simulated for the last time-series 
     376// sample 
     377for(k0=0; k0<Nt; k0++) 
     378  H1(k0) = SPA1(sp[1] 
    346379 
    347380// Free array memory 
    348381free(x0); 
    349 for(k0=0; k0<n; k0++) { 
    350   free(sp[k0].h); 
    351   free(sp[k0].x0); 
    352   free(sp[k0].x1); 
     382for(k0=0; k0<2; k0++) { 
     383  for(k1=0; k1<n; k1++) { 
     384    free(sp[k0][k1].h); 
     385    free(sp[k0][k1].x0); 
     386    free(sp[k0][k1].x1); 
     387  } 
    353388} 
    354389 
    355390// All done 
    356 result[0] = Lx; 
    357 result[1] = Lh; 
     391result[0] = 0; 
     392for(k0=0; k0<Nt; k0++) { 
     393  result[0] += sp[0][k0]->Lx; 
     394result[1] = Lp; 
    358395return_val = result;