Changeset 187
- Timestamp:
- 05/20/08 18:09:32 (7 months ago)
- Files:
-
- projects/AsynCluster/trunk/svpmc/model.c (modified) (11 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
r186 r187 276 276 int k, kpMax; 277 277 double val; 278 double Dp[n]; 278 double Lp[n], Dp[n]; 279 double L_min = +HUGE_VAL; // Anything will be more negative 279 280 double pMax = -HUGE_VAL; // Anything will be more positive 280 281 struct selection result; 281 282 283 // Normalize log-densities to prevent overflow 284 for(k=0; k<n; k++) { 285 val = sp[k][1].Lx + Lz[k]; 286 if(val < L_min) 287 L_min = val; 288 Lp[k] = val; 289 } 290 282 291 // Compute linear probability density for each proposal and their maximum 283 292 for(k=0; k<n; k++) { 284 val = exp( sp[k][1].Lx + Lz[k]);293 val = exp(Lp[k] - L_min); 285 294 Dp[k] = val; 286 295 if(val > pMax) { … … 299 308 300 309 // The probability density of the selected log-volatility proposal is just 301 // the linear weight of the selected sample302 result.Lh = log(Dp[result.kp]);310 // the linear, unnormalized weight of the selected sample 311 result.Lh = Lp[result.kp]; 303 312 return result; 304 313 } … … 339 348 340 349 int notFirst; 341 int k0, k1, k2, k3 , kp;350 int k0, k1, k2, k3; 342 351 double val; 343 352 double Lx = 0; … … 358 367 struct parameters P; 359 368 360 // The results go here369 // The scalar results go here 361 370 py::tuple result(2); 362 371 … … 376 385 sp[k0][k1].x1 = (double *) malloc(sizeof(double) * Nt); 377 386 } 378 // The "previous" log-volatility of the first time-series sample is zero 387 // The "previous" log-volatility of the first time-series sample is set to 388 // the log-volatility offset. 389 // 390 // TODO: Treat it as another latent parameter and simulate its value 379 391 for(k1=0; k1<Nt; k1++) 380 PA1(sp[k0][0].h, k1) = 0;392 PA1(sp[k0][0].h, k1) = D1(k1); 381 393 } 382 394 … … 389 401 k1 = k0; 390 402 notFirst = 0; 391 for(k2=0; k2<n; k2++)392 Lz[k2] = 0;393 403 394 404 do { … … 405 415 // First time-series sample, use a proposal on z instead of z itself 406 416 val += uniform(-w, +w); 407 Lz[k2] += normLogDensity(val);417 Lz[k2] = normLogDensity(val); 408 418 } 409 419 PA1(sp[k2][notFirst].z, k3) = val; … … 432 442 } 433 443 434 } while (k1 ==k0+1 || (k1<Ns && notCloseEnough(sp, n)));444 } while (k1<Ns); //k1==k0+1 || (k1<Ns && notCloseEnough(sp, n))); 435 445 436 446 // Importance sample proposal using pProb=exp(Lz+Lx) for each... … … 453 463 } 454 464 465 // Add the log-likelihood of this time-series sample's innovation to the 466 // running total 467 for(k1=0; k1<Nt; k1++) 468 Lx += sp[k1][0].Lx; 469 455 470 // Done with log-volatility draw for this time-series sample 456 471 } … … 459 474 // sample 460 475 for(k0=0; k0<Nt; k0++) 461 H1(k0) = PA1(sp[ kp][1].h, k0);476 H1(k0) = PA1(sp[spSelection.kp][1].h, k0); 462 477 463 478 // Free array memory … … 473 488 474 489 // The pair of results is Lx, the log-likelihood of each innovation given its 475 // simulated log-volatility... 476 for(k0=0; k0<Nt; k0++) 477 Lx += sp[k0][0].Lx; 478 // ..and Lh, the total log-density of all of the simulated log-volatilities 490 // simulated log-volatility and Lh, the total log-density of all of the 491 // simulated log-volatilities 479 492 result[0] = Lx; 480 493 result[1] = Lh; projects/AsynCluster/trunk/svpmc/model.py
r186 r187 276 276 # ...including the latent parameter of IID normal variates 277 277 z = paramContainer.z.copy() 278 print self.fs, self.fr 279 280 # Set up arary variables for the inline call, including innovations as 278 279 # Set up array variables for the inline call, including innovations as 281 280 # residuals of my observations run through the reversed VAR(1)... 282 281 x = self.var.reverse(self.a, self.b) … … 302 301 Lx = 0 303 302 Lh = 0 304 print "\n\n", z[0,:] , rz, ri303 print "\n\n", z[0,:] 305 304 for k in xrange(self.N2): 306 305 result = self.inline( projects/AsynCluster/trunk/svpmc/test/test_model.py
r186 r187 527 527 spSelection = importanceSample(sp, Lz, n); 528 528 result[0] = spSelection.kp; 529 result[1] = spSelection.L p;529 result[1] = spSelection.Lh; 530 530 return_val = result; 531 531 """ … … 574 574 575 575 class Test_Model_likelihood(BaseCase): 576 N = 6577 576 corr = 0.0 578 577 params = { … … 587 586 } 588 587 589 def _zToh(self, z, e, fs): 590 return signal.lfilter([1], [1.0, -e[0][0]], fs[0]*z) 591 592 def _modelAndParamFactory(self, N1, N2, **kw): 588 def _zToh(self, z, eValue, fsValue): 589 return signal.lfilter([1], [1.0, -eValue], fsValue*z) 590 591 def _simData(self, Ns, eValue, gValue=0, fsValue=1): 592 iv = random.multivariate_normal( 593 [0, 0], 594 [[1.0, gValue], [gValue, 1.0]], Ns).transpose() 595 h = self._zToh(iv[1,:], eValue, fsValue) 596 x = s.exp(0.5 * h) * iv[0,:] 597 return iv, h, x 598 599 def _modelAndParamFactory(self, Ns, N1, N2, **kw): 593 600 # Create a parameter container with the parameters set by my params 594 601 # dict, optionally overridden by keywords, and a zero log-volatility … … 599 606 setattr(pc, name, s.array(kw.get(name, self.params[name]))) 600 607 print pc.paramSequence() 601 pc.z = 3*s.ones((1, self.N))608 pc.z = 0*s.ones((1, Ns)) 602 609 # Simulate a single time series of innovation data per the parameters 603 self.iv = random.multivariate_normal( 604 [0, 0], 605 [[1.0, pc.g[0]], [pc.g[0], 1.0]], self.N).transpose() 606 self.h = self._zToh(self.iv[1,:], pc.e, pc.fs) 607 self.x = s.exp(0.5 * self.h) * self.iv[0,:] 610 self.iv, self.h, self.x = self._simData( 611 Ns, pc.e[0][0], pc.g[0], pc.fs[0]) 608 612 # Run the innovations through the VAR(1) to get the simulated 609 613 # observations … … 615 619 modelObj.y = s.row_stack([y]) 616 620 return pc, modelObj 617 618 def _check(self, wiggle, N1, N2, **kw): 619 pc, modelObj = self._modelAndParamFactory(N1, N2, **kw) 620 modelObj.wiggle = wiggle 621 Lx, Lh, z, hn = modelObj.likelihood(pc) 621 622 def _check(self, x, zReal, zSim, e, fs): 622 623 # Reconstruct simulated h from z 623 h = self._zToh(z[0,:], pc.e, pc.fs) 624 hSim = self._zToh(zSim, e, fs) 625 hReal = self._zToh(zReal, e, fs) 624 626 if VERBOSE: 625 627 fig = self.fig 626 vectors = ( self.x, (self.h, h))628 vectors = (x, (zReal, zSim), (hReal, hSim)) 627 629 for k, vector in enumerate(vectors): 628 630 spNumber = 100*len(vectors) + k + 11 … … 632 634 for v in vector: 633 635 sp.plot(v) 634 self._check_corr(self.h, h) 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 } 651 modelObj = model.Model() 652 Lx = s.empty(N2) 653 Lh = s.empty(N2) 654 print "\n\n", z 655 for k in xrange(self.N2): 656 result = modelObj.inlineDirect(modelObj.code['likelihood'], **kw) 657 print "\n", k, result 658 print z 659 Lx[k] = result[0] 660 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 666 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 635 671 636 def test_highCorrelation(self): 637 self._check(0.1, 10, 10, e=[[0.95]]) 672 673 def test_highCorrelation_model(self): 674 self._modelLikelihood(0.1, 200, 10, 10, e=[[0.95]]) 638 675 639 676 def test_hybridGibbs_lowCorrelation(self):
