Changeset 188
- Timestamp:
- 05/21/08 13:20:34 (7 months ago)
- Files:
-
- projects/AsynCluster/trunk/svpmc/model.c (modified) (10 diffs)
- projects/AsynCluster/trunk/svpmc/model.py (modified) (2 diffs)
- projects/AsynCluster/trunk/svpmc/test/test_model.py (modified) (6 diffs)
Legend:
- Unmodified
- Added
- Removed
- Modified
- Copied
- Moved
projects/AsynCluster/trunk/svpmc/model.c
r187 r188 110 110 // of h to this point 111 111 double Lx; 112 // Change in Lx from this time-series sample to the previous one 113 double L_diff; 112 114 // Multivariate normal proposal, which is just z for time-series samples 113 115 // after the first in each log-volatility computation … … 225 227 226 228 // Now the log-likelihoood. For each time series... 229 S->L_diff = 0; 227 230 for(k=0; k<N; k++) { 228 231 // Offset the decorrelated innovation by the mean, which is the correlation … … 237 240 val = SPA1(S, x0, k) / SPP2(P, sc, k, 0); 238 241 // Fast squaring instead of pow() 239 S->L x-= 0.5 * val*val + SPP2(P, sc, k, 1);242 S->L_diff -= 0.5 * val*val + SPP2(P, sc, k, 1); 240 243 // Since the innovation has been inverse-scaled by the square root of the 241 244 // proposed volatility, the probability density needs to be scaled as 242 245 // well. In log space, this is a subtraction of half the log-volatility. 243 S->Lx -= 0.5 * SPA1(S, h, k); 244 } 246 S->L_diff -= 0.5 * SPA1(S, h, k); 247 } 248 S->Lx += S->L_diff; 245 249 } 246 250 … … 255 259 double Lmax = -HUGE_VAL; // Anything will be more positive 256 260 for(k=0; k<n; k++) { 257 L = sp[k][1].L x;261 L = sp[k][1].L_diff; 258 262 if(L < Lmin) 259 263 Lmin = L; … … 349 353 int notFirst; 350 354 int k0, k1, k2, k3; 351 double val ;355 double val, L_diff; 352 356 double Lx = 0; 353 357 double Lh = 0; … … 369 373 // The scalar results go here 370 374 py::tuple result(2); 371 372 375 373 376 // Generate parameter struct … … 425 428 // ...and the log-likelihood of the innovation for this time-series sample, 426 429 // for each proposal on z and its log-volatility 427 for(k2=0; k2<n; k2++) 430 val = 0; 431 for(k2=0; k2<n; k2++) { 428 432 logLikelihood(&sp[k2][notFirst], &P, xk); 429 430 // Ready for the next time series sample, if we go that far 431 k1++; 433 } 432 434 433 if(k1 == k0+1) {435 if(k1==k0) { 434 436 // Next will be the second time-series sample... 435 437 notFirst = 1; … … 441 443 } 442 444 } 443 444 } while (k1<Ns); //k1==k0+1 || (k1<Ns && notCloseEnough(sp, n))); 445 445 446 // Ready for the next time series sample, assuming there is one 447 k1++; 448 449 } while (k1<Ns && (k1==k0+1 || notCloseEnough(sp, n))); 450 446 451 // Importance sample proposal using pProb=exp(Lz+Lx) for each... 447 452 spSelection = importanceSample(sp, Lz, n); … … 454 459 Z2(k1,k0) = PA1(sp[spSelection.kp][0].z, k1); 455 460 461 //--- DEBUG ----------------------------------------------------------------- 462 printf("\n%d %d\n---------------------------\n",k0, spSelection.kp); 463 for(k1=0; k1<n; k1++) 464 printf("%d %f %f\n", k1, sp[k1][1].Lx, Lz[k1]); 465 //--------------------------------------------------------------------------- 466 456 467 // Set the initial log-volatility of the proposals for the next time-series 457 468 // sample to the log-volatility of the selected proposal for this time-series … … 464 475 465 476 // Add the log-likelihood of this time-series sample's innovation to the 466 // running total 467 for(k1=0; k1<Nt; k1++) 477 // running total and clear the Lx accumulators for the next time-series 478 // sample 479 for(k1=0; k1<Nt; k1++) { 468 480 Lx += sp[k1][0].Lx; 481 sp[k1][0].Lx = 0; 482 sp[k1][1].Lx = 0; 483 } 469 484 470 485 // Done with log-volatility draw for this time-series sample projects/AsynCluster/trunk/svpmc/model.py
r187 r188 287 287 try: 288 288 rz = linalg.cholesky( 289 self.covarMatrix(self.fs, self.fr), lower=True) 290 ri = linalg.inv(linalg.cholesky( 291 s elf.covarMatrix(s.ones(self.p), self.cr), lower=True))289 self.covarMatrix(self.fs, self.fr), lower=True).astype('d') 290 ri = linalg.inv(linalg.cholesky(self.covarMatrix( 291 s.ones(self.p), self.cr), lower=True)).astype('d') 292 292 except linalg.LinAlgError: 293 293 return … … 301 301 Lx = 0 302 302 Lh = 0 303 print "\n\n", z[0,:]304 303 for k in xrange(self.N2): 305 304 result = self.inline( 306 305 'z', 'h', 'x', 'w', 'n', 'd', 'e', 'g', 'rz', 'ri', 'sc', 307 306 w=self.wiggle, n=self.N1) 308 print "\n", k, result 309 print z[0,:] 307 print k, result 310 308 Lx += result[0] 311 309 Lh += result[1] projects/AsynCluster/trunk/svpmc/test/test_model.py
r187 r188 511 511 const int n = NX[0]; 512 512 struct sample sp[n][2]; 513 for(k=0; k<n; k++) 513 for(k=0; k<n; k++) { 514 sp[k][1].L_diff = X1(k); 514 515 sp[k][1].Lx = X1(k); 516 } 515 517 """ 516 518 … … 581 583 'd': [0], 582 584 'e': [[0]], 583 'fs': [ 0.5],585 'fs': [1.0], 584 586 'fr': [], 585 587 'g': [corr], 586 588 } 587 589 588 def _zToh(self, z, eValue, fsValue): 590 def _z_to_h(self, z, eValue, fsValue): 591 """ 592 Runs the volatility shocks through a VAR(1) to get simulated 593 log-volatilities. 594 """ 589 595 return signal.lfilter([1], [1.0, -eValue], fsValue*z) 596 597 def _x_to_y(self, x, aValue, bValue): 598 """ 599 Runs the innovations through a VAR(1) to get simulated 600 observations. 601 """ 602 return signal.lfilter([1], [1.0, -bValue], x) + aValue * s.ones_like(x) 590 603 591 604 def _simData(self, Ns, eValue, gValue=0, fsValue=1): … … 593 606 [0, 0], 594 607 [[1.0, gValue], [gValue, 1.0]], Ns).transpose() 595 h = self._z Toh(iv[1,:], eValue, fsValue)608 h = self._z_to_h(iv[1,:], eValue, fsValue) 596 609 x = s.exp(0.5 * h) * iv[0,:] 597 610 return iv, h, x 598 611 599 def _modelAndParamFactory(self, Ns, N1, N2, **kw): 612 def _modelAndParamFactory(self, wiggle, Ns, N1, N2, **kw): 613 """ 614 Returns a parameter container populated from my I{params} dict, 615 optionally overrided by keywords, along with an instance of 616 L{model.Model} to be tested. 617 """ 600 618 # Create a parameter container with the parameters set by my params 601 619 # dict, optionally overridden by keywords, and a zero log-volatility … … 604 622 pc.paramNames = util.PARAM_NAMES 605 623 for name, value in self.params.iteritems(): 606 setattr(pc, name, s.array(kw.get(name, self.params[name]))) 607 print pc.paramSequence() 608 pc.z = 0*s.ones((1, Ns)) 609 # Simulate a single time series of innovation data per the parameters 610 self.iv, self.h, self.x = self._simData( 611 Ns, pc.e[0][0], pc.g[0], pc.fs[0]) 612 # Run the innovations through the VAR(1) to get the simulated 613 # observations 614 y = signal.lfilter([1], [1.0, -pc.b[0][0]], self.x) + \ 615 pc.a * s.ones_like(self.x) 624 setattr( 625 pc, name, 626 s.array(kw.get(name, self.params[name]), dtype='d')) 616 627 # Return the parameter container along with a model to be tested 617 628 modelObj = model.Model(N1=N1, N2=N2) … … 620 631 return pc, modelObj 621 632 622 def _check(self, x, zReal, zSim, e , fs):633 def _check(self, x, zReal, zSim, eValue, fsValue=1.0): 623 634 # Reconstruct simulated h from z 624 hSim = self._z Toh(zSim, e, fs)625 hReal = self._z Toh(zReal, e, fs)635 hSim = self._z_to_h(zSim, e, fs) 636 hReal = self._z_to_h(zReal, e, fs) 626 637 if VERBOSE: 627 638 fig = self.fig … … 634 645 for v in vector: 635 646 sp.plot(v) 636 self._check_corr(hReal, hSim) 637 638 def _likelihood(self, z, x, eValue, wiggle, N1, N2): 639 kw = {'h':s.empty(1), 640 'sc':s.array([1, s.log(s.sqrt(2*s.pi))]), 641 'z':s.row_stack([z]), 642 'x':s.row_stack([x]), 643 'w':wiggle, 644 'n':N1, 645 'd':s.array([0]), 646 'e':s.array([[eValue]]), 647 'g':s.array([0]), 648 'rz':s.array([[1]]), 649 'ri':s.array([[1]]) 650 } 647 self._check_corr(hReal, hSim) 648 649 def _likelihood(self, x, eValue, wiggle, N1, N2): 650 """ 651 Runs C code for L{model.Model.likelihood} independent of the method 652 itself. 653 """ 654 z = s.zeros((1, len(x))) 655 zAll = s.empty((N2, len(x))) 656 kw = { 657 'h':s.empty(1), 658 'sc':s.array([1, s.log(s.sqrt(2*s.pi))]), 659 'z':z, 660 'x':s.row_stack([x]), 661 'w':wiggle, 662 'n':N1, 663 'd':s.array([0.0]), 664 'e':s.array([[float(eValue)]]), 665 'g':s.array([0.0]), 666 'rz':s.array([[1.0]]), 667 'ri':s.array([[1.0]]) 668 } 651 669 modelObj = model.Model() 652 670 Lx = s.empty(N2) 653 671 Lh = s.empty(N2) 654 print "\n\n", z 655 for k in xrange(self.N2): 672 for k in xrange(N2): 656 673 result = modelObj.inlineDirect(modelObj.code['likelihood'], **kw) 657 print "\n", k, result658 print z659 674 Lx[k] = result[0] 660 675 Lh[k] = result[1] 661 return Lx, Lh, z 662 663 def _modelLikelihood(self, wiggle, Ns, N1, N2, **kw): 664 pc, modelObj = self._modelAndParamFactory(Ns, N1, N2, **kw) 665 modelObj.wiggle = wiggle 676 zAll[k,:] = z.copy() 677 return Lx, Lh, zAll 678 679 def test_basic(self): 680 wiggle = 0.1 681 eValue = 0.90 682 Ns, N1, N2 = 3, 5, 10 683 x = 2.0*s.ones(Ns) 684 Lx, Lh, zSim = self._likelihood(x, eValue, wiggle, N1, N2) 685 #fig = self.fig 686 #fig.add_subplot(211).plot(Lx) 687 #fig.add_subplot(212).plot(Lh) 688 print "\nk Lx Lh Simulated Volatility Shocks" 689 print "-"*100 690 for k, z in enumerate(zSim): 691 zLine = " ".join(["%+4.2f" % xx for xx in z]) 692 print "%d: %+5.1f %+6.1f%s%s" % (k, Lx[k], Lh[k], " "*5, zLine) 693 694 def test_highCorrelation_model(self): 695 wiggle = 0.1 696 eValue = 0.95 697 Ns, N1, N2 = 1000, 10, 100 698 iv, h, x = self._simData(Ns, eValue) 699 pc, modelObj = self._modelAndParamFactory(wiggle, Ns, N1, N2, e=[[0.95]]) 700 modelObj.y = self._x_to_y(x) 666 701 Lx, Lh, z, hn = modelObj.likelihood(pc) 667 self._check(self.x, self.iv[1,:], z[0,:], pc.e, pc.fs) 668 669 def test_basic(self): 670 eValue = 0.95 671 672 673 def test_highCorrelation_model(self): 674 self._modelLikelihood(0.1, 200, 10, 10, e=[[0.95]]) 702 self._check(x, iv[1,:], z[0,:], eValue) 675 703 676 704 def test_hybridGibbs_lowCorrelation(self):
