Changeset 183

Show
Ignore:
Timestamp:
05/16/08 16:40:27 (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

    r182 r183  
    107107struct sample 
    108108{ 
    109   // Log-density of the proposal on z 
    110   double Lz; 
    111   // Log-likelihood of the current innovation, given h 
    112   double Lx; 
    113   // Linear probability density of the proposal, given h 
    114   double pProb; 
    115   // Multivariate normal proposal 
     109  // Log-likelihood of the innovations up to the current one, given the history 
     110  // of h to this point 
     111  double Lx = 0; 
     112  // Multivariate normal proposal, which is just z for time-series samples 
     113  // after the first in each log-volatility computation 
    116114  double *z; 
    117115  // Multivariate log-volatility value based on the current proposal, or the 
    …… 
    124122 
    125123 
    126 // Uniformly distributed random variate in (a,b) 
     124// A selection of a proposal sample 
     125struct selection 
     126{ 
     127  int kp; 
     128  double Lp; 
     129} 
     130 
     131 
     132// Uniformly distributed random variate in [a,b) 
    127133double uniform(double a, double b) 
    128134{ 
    …… 
    202208  double *y; 
    203209  double val; 
    204   double L = 0; 
    205   const int N = (int) P->d->dimensions[0]; 
    206  
     210  const int N = (int) P->g->dimensions[0]; 
     211   
    207212  // Inverse-scale the innovation by the square root of the volatilities in 
    208213  // preparation for decorrelation 
    …… 
    213218  // log-likelihood computation for this multivariate log-volatility proposal 
    214219  matrixMultiply(P->ri, S->x0, S->x1); 
    215  
     220   
    216221  // Now compute log Pr(x[t]|h[t]), the log-likelihood of the decorrelated 
    217222  // multivariate innovation for this sample, conditional on the sample's 
    …… 
    230235    // scaling coefficient 
    231236    val = SPA1(S, x0, k) / SPP2(P, sc, k, 0); 
    232     L -= 0.5 * val*val + SPP2(P, sc, k, 1); // Fast squaring instead of pow() 
     237    // Fast squaring instead of pow() 
     238    S->Lx -= 0.5 * val*val + SPP2(P, sc, k, 1); 
    233239    // Since the innovation has been inverse-scaled by the square root of the 
    234240    // proposed volatility, the probability density needs to be scaled as 
    235241    // well. In log space, this is a subtraction of half the log-volatility. 
    236     L -= 0.5 * SPA1(S, h, k); 
    237   } 
    238  
    239   // All done 
    240   S->Lx = L; 
     242    S->Lx -= 0.5 * SPA1(S, h, k); 
     243  } 
    241244} 
    242245 
    …… 
    244247// Returns with zero only if the post-first proposals are close enough to 
    245248// finish likelihood computation 
    246 int notCloseEnough(struct sample *sp[2][], int n) 
     249int notCloseEnough(struct sample *sp[][2], int n) 
    247250{ 
    248251  int k; 
    …… 
    251254  double Lmax = -1E20;   // Anything will be more positive 
    252255  for(k=0; k<n; k++) { 
    253     L = sp[1][k]->Lx; 
     256    L = sp[k][1]->Lx; 
    254257    if(L < Lmin) 
    255258      Lmin = L; 
    …… 
    258261  } 
    259262  return (Lmax - Lmin > 0.001); 
     263} 
     264 
     265 
     266// Importance-samples one of the proposals based on a weight computed from the 
     267// log-density of z and the log-likelihood of the innovation given the 
     268// proposal's log-volatility. 
     269// 
     270// Returns a struct with an integer index to the sampled proposal's struct and 
     271// the log-probability of the proposal's having been selected. 
     272struct selection \ 
     273sample(struct sample *sp[][2], double Lz[], double Dp[], int n) 
     274{ 
     275  int k0, kpMax; 
     276  double val; 
     277  double pSum = 0; 
     278  double pMax = -1E20;  // Anything will be more positive 
     279  struct selection result; 
     280 
     281  // Compute linear probability density for each proposal 
     282  for(k0=0; k0<n; k0++) { 
     283    val = exp(sp[k0][1]->Lx + Lz[k0]); 
     284    Dp[k0] = val; 
     285    pSum += val; 
     286    if(val > pMax) { 
     287      pMax = val; 
     288      kpMax = k0; 
     289  } 
     290 
     291  // Accept-reject sampling with uniform random proposal selection 
     292  do { 
     293    // Randomly select one of the proposals... 
     294    result.kp = (int) uniform(0, n); 
     295    // ...until the selected proposal is accepted, with probability 
     296    // proportional to its relative weight 
     297  } while(result.kp != kpMax && uniform(0, pMax) > Dp[result.kp]); 
     298 
     299  // Only one division operation is needed! 
     300  result.Lp = log(Dp[result.kp] / (pMax*pSum)); 
     301  return result; 
    260302} 
    261303 
    …… 
    288330// 1    Lh: Log-likelihood of simulated log-volatilities 
    289331 
     332 
     333// The number of time series making up a single multivariate sample 
    290334const int Nt = Nz[0]; 
     335// The number of multivariate samples 
    291336const int Ns = Nz[1]; 
    292337 
    293338int notFirst; 
    294339int k0, k1, k2, k3, kp; 
    295 double zp; 
     340double val; 
     341double Lx = 0; 
    296342double Lp = 0; 
    297343 
     344// Multivariate innovation for one current time-series sample 
    298345double *xk; 
    299 struct sample sp[2][n]; 
     346// Log probability density of each proposal on z 
     347double Lz[n]; 
     348// Linear probability density of each log-volatility proposal 
     349double Dp[n]; 
     350 
     351// A struct for the result of the importance-sampling function 
     352struct selection spSelection; 
     353// Two sets of proposal sample structs, one for the first time-series sample in 
     354// each log-volatility computation (from which the innovation's log-volatility 
     355// is computed) and another for all time-series samples thereafter. 
     356struct sample sp[n][2]; 
     357// A struct for the model parameters 
    300358struct parameters P; 
     359 
     360// The results go here 
    301361py::tuple result(2); 
    302362 
    …… 
    309369xk = (double *) calloc(sizeof(double), Nt); 
    310370for(k0=0; k0<n; k0++) { 
    311   for(k1=0; k1<n; k1++) { 
     371  for(k1=0; k1<2; k1++) { 
    312372    sp[k0][k1].h = (double *) calloc(sizeof(double), Nt); 
    313373    sp[k0][k1].x0 = (double *) calloc(sizeof(double), Nt); 
    …… 
    316376} 
    317377 
    318 // For each sample... 
     378// For each time-series sample... 
    319379for(k0=0; k0<Ns; k0++) { 
    320380   
    …… 
    323383 
    324384  k1 = k0; 
     385  notFirst = 0; 
     386  for(k2=0; k2<n; k2++) 
     387    Lz[k2] = 0; 
     388   
    325389  do { 
    326390 
    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 
     391    // Prepare current-innovation vector xk 
    330392    for(k2=0; k2<Nt; k2++) 
    331393      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... 
     394     
     395    // Set z for the set of sample structs... 
    349396    for(k2=0; k2<n; k2++) { 
    350397      for(k3=0; k3<Nt; k3++) { 
    351         zp = Z2(k3, k1) + uniform(-w, +w); 
    352         SPA1(sp[notFirst][k2], z, k3) = zp; 
    353         sp[notFirst][k2].Lz = normLogDensity(zp); 
     398        val = Z2(k3, k1); 
     399        if(!notFirst) { 
     400          // First time-series sample, use a proposal on z instead of z itself 
     401          val += uniform(-w, +w); 
     402          Lz[k2] += normLogDensity(val); 
     403        } 
     404        SPA1(sp[k2][notFirst], z, k3) = val; 
    354405      } 
    355406    } 
    356407    // ...compute their log-volatilities 
    357408    for(k2=0; k2<n; k2++) 
    358       logVol(&sp[notFirst][k2], &P); 
     409      logVol(&sp[k2][notFirst], &P); 
    359410    // ...and the log-likelihood of the innovation for this time-series sample, 
    360411    // for each proposal on z and its log-volatility 
    361412    for(k2=0; k2<n; k2++) 
    362       logLikelihood(&sp[notFirst][k2], &P, xk); 
     413      logLikelihood(&sp[k2][notFirst], &P, xk); 
    363414 
    364415    // Ready for the next time series sample, if we go that far 
    365416    k1++; 
    366  
     417     
     418    if(k1 == k0+1) { 
     419      // Next will be the second time-series sample... 
     420      notFirst = 1; 
     421      // ...and will be starting with copy of h arrays from proposal structs 
     422      // for the just-completed first time-series sample 
     423      for(k2=0; k2<n; k2++) { 
     424        for(k3=0; k3<Nt; k3++) 
     425          SPA1(sp[k2][1], h, k3) = SPA1(sp[k2][0], h, k3); 
     426      } 
     427    } 
     428   
    367429  } while (k1==k0+1 || (k1<Ns && notCloseEnough(sp, n))); 
    368430 
    369431  // Importance sample proposal using pProb=exp(Lz+Lx) for each 
    370   kp = sample(sp, n); 
    371   Lp += log(sp[1][kp]->pProb); 
    372  
     432  spSelection = sample(sp, Lz, Dp, n); 
     433  Lp += log(spSelection.Lp); 
     434  for(k1=0; k1<Nt; k1++) 
     435    Z2(k1,k0) = SPA1(sp[spSelection.kp][1], z, k1); 
     436 
     437  // Set the initial log-volatility of the proposals for the next time-series 
     438  // sample to the log-volatility of the selected proposal for this time-series 
     439  // sample. 
     440  for(k1=0; k1<n; k1++) { 
     441    val = SPA1(sp[spSelection.kp][1], h, k1); 
     442    for(k2=0; k2<Nt; k2++) 
     443      SPA1(sp[k1][0], h, k2) = val; 
     444  } 
     445 
     446  // Done with log-volatility draw for this time-series sample 
    373447} 
    374448 
    …… 
    376450// sample 
    377451for(k0=0; k0<Nt; k0++) 
    378   H1(k0) = SPA1(sp[1] 
     452  H1(k0) = SPA1(sp[kp][1], h, k0); 
    379453 
    380454// Free array memory 
    381455free(x0); 
    382 for(k0=0; k0<2; k0++) { 
    383   for(k1=0; k1<n; k1++) { 
     456for(k0=0; k0<n; k0++) { 
     457  for(k1=0; k1<2; k1++) { 
    384458    free(sp[k0][k1].h); 
    385459    free(sp[k0][k1].x0); 
    …… 
    391465result[0] = 0; 
    392466for(k0=0; k0<Nt; k0++) { 
    393   result[0] += sp[0][k0]->Lx; 
     467  result[0] += sp[k0][0]->Lx; 
    394468result[1] = Lp; 
    395469return_val = result; 

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

    r181 r183  
    358358 
    359359    code = """ 
    360     int k; 
     360    int k0, k1; 
    361361    struct sample S; 
    362362    struct parameters P; 
    …… 
    367367    S.x0 = (double *) malloc(sizeof(double) * Nz[0]); 
    368368    S.x1 = (double *) malloc(sizeof(double) * Nz[0]); 
    369     for(k=0; k<Nz[0]; k++) 
    370       PA1(S.z, k) = Z1(k); 
    371  
    372     logVol(&S, &P); 
    373  
    374     for(k=0; k<Nz[0]; k++) 
    375       H1(k) = PA1(S.h, k); 
     369     
     370    for(k0=0; k0<Nz[1]; k0++) { 
     371      for(k1=0; k1<Nz[0]; k1++) 
     372        PA1(S.z, k1) = Z2(k1, k0); 
     373      // This is what's being tested 
     374      logVol(&S, &P); 
     375      for(k1=0; k1<Nz[0]; k1++) 
     376        H2(k1, k0) = PA1(S.h, k1); 
     377    } 
    376378    """ 
    377379     
    378     def _draw_h(self, z, d, e, fs, fr): 
     380    def _draw_h(self, d, e, fs, fr, z=None): 
     381        if z is None: 
     382            z = s.randn(len(d), self.N) 
    379383        h = s.zeros_like(z) 
    380384        rz = s.array(self._rv(fs, fr)) 
    381         for k, zk in enumerate(z.transpose()): 
    382             self.inline( 
    383                 self.code, z=zk, 
    384                 d=s.array(d), e=s.array(e), h=h[:,k], rz=rz) 
     385        self.inline( 
     386            self.code, z=z, d=s.array(d), e=s.array(e), h=h, rz=rz) 
    385387        return h 
    386  
     388     
    387389    def test_univariate_LPF(self): 
    388390        e = 0.95 
    389391        # Univariate model 
    390         z = s.randn(1, self.N) 
    391         h = self._draw_h(z, d=[1], e=[[e]], fs=[1], fr=[]) 
     392        h = self._draw_h([1], [[e]], [1], []) 
    392393        # IIR filtering 
    393394        self._check_psd(h, e) 
    394395 
    395396    def test_indep_offset(self): 
    396         modelObj = self.modelFactory( 
    397             d = s.array([1.0, 1.0]), 
    398             fs = [1.0, 1.0], 
    399             fr = [0.0]) 
     397        d = s.array([1.0, 1.0]) 
     398        fs = [1.0, 1.0] 
     399        fr = [0.0] 
    400400        for j in xrange(10): 
    401401            e1, e2 = 0.8*s.rand(2) + 0.1 
    402             modelObj.e = s.array([[e1, 0.0], [0.0, e2]]) 
     402            e = s.array([[e1, 0.0], [0.0, e2]]) 
    403403            for k in xrange(20): 
    404                 h = self._draw_h(modelObj) 
     404                h = self._draw_h(d, e, fs, fr) 
    405405            for m, em in enumerate((e1, e2)): 
    406406                ratio = h[m,self.N/2:].mean() /  ( 1.0/(1-em) ) 
    407407                self.failUnlessAlmostEqual(ratio, 1.0, 0) 
    408  
     408     
    409409    def test_indep_shocks(self): 
    410410        e1, e2 = 0.95, 0.8 
    411         modelObj = self.modelFactory( 
    412             d = s.array([0.0, 0.0]), 
    413             e = s.array([[e1, 0.0], [0.0, e2]]), 
    414             fs = [1.0, 1.0], 
    415             fr = [0.0]) 
    416         h = self._draw_h(modelObj) 
     411        d = s.array([0.0, 0.0]) 
     412        e = s.array([[e1, 0.0], [0.0, e2]]) 
     413        fs = [1.0, 1.0] 
     414        fr = [0.0] 
     415        h = self._draw_h(d, e, fs, fr) 
    417416        # Uncorrelated 
    418417        self._check_uncorr(h[0,:], h[1,:]) 
    …… 
    420419        self._check_psd(h, e1, e2) 
    421420 
    422     def _check_xcorr(self, dPeak, **kw): 
    423         modelObj = self.modelFactory(**kw) 
    424         h = self._draw_h(modelObj) 
     421    def _check_xcorr(self, dPeak, *args): 
     422        h = self._draw_h(*args) 
    425423        # Correlation 
    426424        nd = 6 
    …… 
    440438        self._check_xcorr( 
    441439            1, 
    442             d = [0.0, 0.0], 
    443             e = [[0.0, 0.9], [0.0, 0.0]], 
    444             fs = [1.0, 1.0], 
    445             fr = [0.0], 
     440            [0.0, 0.0], 
     441            [[0.0, 0.9], [0.0, 0.0]], 
     442            [1.0, 1.0], 
     443            [0.0], 
    446444            ) 
    447445 
    …… 
    449447        self._check_xcorr( 
    450448            0, 
    451             d = [0.0, 0.0], 
    452             e = [[0.0, 0.0], [0.0, 0.0]], 
    453             fs = [1.0, 1.0], 
    454             fr = [0.5], 
     449            [0.0, 0.0], 
     450            [[0.0, 0.0], [0.0, 0.0]], 
     451            [1.0, 1.0], 
     452            [0.5], 
    455453            ) 
    456454 
    457455 
    458 class Test_Model_logVolatilities_decorrelate(BaseCase): 
    459     N = 1000 
    460     y = s.randn(1, N) 
    461  
    462     def setUp(self): 
    463         self.model = self.modelFactory() 
    464      
    465     def _pGenerator(self, low, high): 
    466         for p in xrange(low, high): 
    467             # Tweak the model object to make it testable for this p value 
    468             self.model.y = s.zeros((p, self.N)) 
    469             self.model.d = s.zeros(p) 
    470             self.model.e = s.zeros((p, p)) 
    471             self.model.g = s.zeros(p) 
    472             # Generate some bogus arrays 
    473             z = s.zeros((p, 1)) 
    474             v = s.zeros((p, 1)) 
    475             h0 = s.zeros((p, 1)) 
    476             # Yield a container with p and room for one or more 2-tuples 
    477             # containing a correlation parameter and arrays of innovations to 
    478             # decorrelate and test using that parameter value 
    479             container = [p] 
    480             yield container 
    481             # Use the logVolatilities method to decorrelate the innovations 
    482             # that the caller put into the container I just yielded 
    483             for self.model.cr, x in container[1:]: 
    484                 xd = s.zeros_like(x) 
    485                 for k in xrange(x.shape[1]): 
    486                     self.model.logVolatilities(z, v, x[:,k], h0) 
    487                     y = self.model.cache['lv_work'][0] 
    488                     xd[:,k] = y[:,3] 
    489                 # Check for non-correlation 
    490                 for j in xrange(p): 
    491                     for k in xrange(p): 
    492                         if j == k: 
    493                             continue 
    494                         self._check_uncorr(xd[j,:], xd[k,:], plot=False) 
    495      
    496     def test_basic(self): 
    497         for pc in self._pGenerator(2, 4): 
    498             p = pc[0] 
    499             for correlation in s.linspace(0.1, 0.9, 5): 
    500                 # Generate some multivariate normal test data 
    501                 r = s.maximum(s.identity(p), correlation*s.ones((p,p))) 
    502                 if VERBOSE: 
    503                     print "\n\n", p 
    504                     print r 
    505                 c = correlation * s.ones(sum([p-k-1 for k in xrange(p)])) 
    506                 x = random.multivariate_normal( 
    507                     s.zeros(p), r, self.N).transpose() 
    508                 pc.append((c, x)) 
    509  
    510     def test_complex(self): 
    511         for pc in self._pGenerator(2, 7): 
    512             p = pc[0] 
    513             # Stuff specific to this number of time series 
    514             N_correlations = sum([p-k-1 for k in xrange(p)]) 
    515             cMax = 0.90 / s.sqrt(p/1.5) 
    516             for i in xrange(10): 
    517                 # Randomly generate a set of cross-correlations and set the 
    518                 # model's innovation shock correlations parameter 
    519                 c = 0.05 + cMax * s.rand(N_correlations) 
    520                 # Generate some multivariate normal test data having the 
    521                 # generated correlations 
    522                 r = self._covarMatrix(s.ones(p), c) 
    523                 if VERBOSE: 
    524                     print "\n\n", p 
    525                     print r 
    526                 x = random.multivariate_normal( 
    527                     s.zeros(p), r, self.N).transpose() 
    528                 pc.append((c, x)) 
    529  
    530  
    531 class Test_Model_logVolatilities(BaseCase): 
     456class Test_Model_logLikelihood(BaseCase): 
    532457    N = 100 
    533  
    534     def setUp(self): 
    535         self.model = self.modelFactory() 
     458     
     459    code = """ 
     460    int k0, k1; 
     461    struct sample S; 
     462    struct parameters P; 
     463 
     464    P.g = g_array; P.ri = ri_array; P.sc = sc_array; 
     465    S.z = (double *) malloc(sizeof(double) * Nz[0]); 
     466    S.h = (double *) calloc(sizeof(double), Nz[0]); 
     467    S.x0 = (double *) malloc(sizeof(double) * Nz[0]); 
     468    S.x1 = (double *) malloc(sizeof(double) * Nz[0]); 
     469     
     470    for(k0=0; k0<Nz[1]; k0++) { 
     471      for(k1=0; k1<Nz[0]; k1++) 
     472        PA1(S.z, k1) = Z2(k1, k0); 
     473      // This is what's being tested 
     474      logLikelihood(&S, &P, x); 
     475      for(k1=0; k1<Nz[0]; k1++) 
     476        H2(k1, k0) = PA1(S.h, k1); 
     477    } 
     478    """ 
     479     
     480    def _draw_h(self, d, e, fs, fr, z=None): 
     481        if z is None: 
     482            z = s.randn(len(d), self.N) 
     483        h = s.zeros_like(z) 
     484        rz = s.array(self._rv(fs, fr)) 
     485        self.inline( 
     486            self.code, z=z, d=s.array(d), e=s.array(e), h=h, rz=rz) 
     487        return h 
    536488 
    537489    def _check_L(self, distObj, hValue):