Changeset 177
- Timestamp:
- 05/10/08 22:37:06 (7 months ago)
- Files:
-
- projects/AsynCluster/trunk/svpmc/model.c (modified) (3 diffs)
- projects/AsynCluster/trunk/svpmc/model.py (modified) (6 diffs)
- projects/AsynCluster/trunk/svpmc/test/test_model.py (modified) (2 diffs)
- projects/AsynCluster/trunk/svpmc/test/util.py (modified) (2 diffs)
Legend:
- Unmodified
- Added
- Removed
- Modified
- Copied
- Moved
projects/AsynCluster/trunk/svpmc/model.c
r158 r177 83 83 84 84 85 // Model.support 86 // ---------------------------------------------------------------------------- 87 88 // Normal (zero mean, unit variance) probability density 89 double npdf(double x) 90 { 91 // 1/sqrt(2*pi) = 0.39894228... 92 return 0.398942280401 * exp(-0.5 * pow(x, 2)); 93 } 94 95 96 // Normal (zero mean, unit variance) cumulative density 97 double ncdf(double x) 98 { 99 double ax, t, sum; 100 // The cdf is computed from the absolute value of x and corrected at the end 101 // if necessary 102 ax = fabs(x); 103 // Compute t, the variable used in the series summation 104 t = 1.0 / (1.0 + 0.2316419*ax); 105 // Efficiently compute the series summation 106 sum = t * +0.31938; // t 107 sum += (t *= t) * -0.35656; // t^2 108 sum += (t *= t) * +1.78148; // t^3 109 sum += (t *= t) * -1.82125; // t^4 110 sum += (t *= t) * +1.33027; // t^5 111 // Compute the cdf, assuming x to be positive 112 sum *= (1.0 - npdf(ax)); 113 // Correct if x < 0 114 if(x < 0) 115 sum = 1 - sum; 116 // All done 117 return sum; 118 } 119 120 121 // Truncated normal (zero mean, unit variance) probability density 122 double tnpdf(double x, double a, double b) 123 { 124 double y; 125 if(x < a || x > b) { 126 y = 0; 127 } else { 128 y = npdf(x) / (ncdf(b) - ncdf(a)); 129 } 130 return y; 131 } 132 133 134 85 135 // Model.zProposal 86 136 // … … 88 138 // ---------------------------------------------------------------------------- 89 139 // z Independent normal variates, [pn] array (updated) 90 // M Number of uniform-proposal MCMC iterations per IID normal draw91 // wiggle = Width of random-walk uniform proposal 92 93 94 int k0, k1 , k2;95 double norm, normp, L, Lp, dL;140 // wiggle = half-width of independent uniform proposal 141 142 #define w (double)wiggle 143 144 int k0, k1; 145 double x, zk, c, cp, L; 96 146 97 147 // Seed the PRNG on the first run only … … 111 161 } 112 162 163 // The joint log-density of all the jumps 164 L = 0; 113 165 // for each multivariate sample... 114 166 for(k0=0; k0<Nz[1]; k0++) { 115 167 116 // Do a random walk from this column of independent random-walk normal117 // variates to produce a multivariate iid normal proposal. It is this on118 // which we build a multivariate log-volatility proposal for this119 // multivariate sample.168 // Do a random walk (a single uniformly distributed step) from this column of 169 // independent random-walk normal variates to produce a multivariate iid 170 // normal proposal. It is this on which we build a multivariate 171 // log-volatility proposal for this multivariate sample. 120 172 121 173 // For each time series... 122 174 for(k1=0; k1<Nz[0]; k1++) { 123 // Initialize the current normal value with the current proposal value and 124 // compute its log-likelihood 125 norm = Z2(k1, k0); 126 L = -0.5 * pow(norm, 2); 127 // For each oversampling iteration... 128 for(k2=0; k2<M; k2++) { 129 // Generate a uniform, symmetric, random walk proposal 130 normp = norm + (double) wiggle * (drand48() - 0.5); 131 // Do the ratio in log space 132 Lp = -0.5 * pow(normp, 2); 133 dL = Lp - L; 134 // The uniform random variate thus needs to be in logspace as well, but 135 // it's not always even needed. Thus the log operation is done just 136 // once, and only some of the time. 137 if(dL > 0 || log(drand48()) < dL) { 138 // Update current value and its log probability when proposal 139 // accepted 140 norm = normp; 141 L = Lp; 142 } 175 // Compute the scale value c as a bit more* than the maximum pdf (regular 176 // normal distribution) at the endpoints and middle of the truncated range 177 c = 0; 178 zk = Z2(k1, k0); 179 for(x = -w; x < w; x += w) { 180 cp = 1.12 * npdf(zk + x); // *Scaling by 1.12 works for w <= 1.0 181 if(cp > c) 182 c = cp; 143 183 } 144 // Update this element of the IID normal array in place 145 Z2(k1, k0) = norm; 146 } 147 } 184 // Accept-reject sampling for the truncated normal distribution 185 do { 186 // Generate a proposal within the truncated range 187 x = zk + 2*w*(drand48() - 0.5); 188 // Compute the pdf of the proposal, using the regular normal distribution 189 cp = npdf(x); 190 // Accept (and exit the loop) when the scaled uniform variate is less 191 // than the proposal's pdf 192 } while(c*drand48() > cp); 193 // Update this element of the IID normal array in place with the accepted 194 // proposal 195 Z2(k1, k0) = x; 196 // Add the log-density of the accepted proposal 197 L += log(tnpdf(x, zk-w, zk+w)); 198 } 199 } 200 201 return_val = L; 148 202 149 203 projects/AsynCluster/trunk/svpmc/model.py
r176 r177 21 21 22 22 import scipy as s 23 from scipy import linalg 23 from scipy import linalg, stats 24 24 25 25 from twisted.internet import defer … … 257 257 return logU[self._logU_index] 258 258 259 def zProposal(self, z, sigma): 260 """ 261 Does a random walk from the supplied 2-D array I{z} of normal random 262 variates, using a symmetric normal proposal distribution with standard 263 deviation I{sigma}. 259 def zProposal(self, z, wiggle): 260 """ 261 Does a random walk step from the supplied 2-D array I{z} of zero-mean, 262 unit-variance normal random variates with a jump from a uniform 263 distribution having width +/- I{wiggle}. 264 265 Returns a copy of the array with all elements changed from the original 266 values, independently, along with the total log-density of the jumps. 264 267 """ 265 268 z = z.copy() 266 self.inline('z', 'M', 'wiggle', M=self.Mz, wiggle=sigma)267 return z 269 L = self.inline('z', 'M', 'wiggle', M=self.Mz) 270 return z, L 268 271 269 272 def logVolatilities(self, z, v, x, h0): … … 322 325 """ 323 326 L_total = 0 324 L_result = 0325 327 acceptances = 0 326 328 N = self.decaySamples(tol) … … 332 334 # ...and a set of random-walk proposals for the IID variates, along 333 335 # with their volatility shocks 334 zp = self.zProposal(z, sigma)336 zp, L_result = self.zProposal(z, sigma) 335 337 vp = s.array(rv*zp) 336 338 # Do conditional draws of each sample in turn … … 412 414 3. The IID normal variates underlying the log-volatilities. 413 415 414 4. The log-probability of the simulated result. 416 4. The log-probability of the simulated result, normalized to be 417 independent of the number of hybrid Gibbs steps. 415 418 416 419 5. The mean acceptance rate of the Metropolis-Hastings proposals. … … 447 450 except linalg.LinAlgError: 448 451 return 449 return info[0], info[1], z, L_result , sum(acceptRates)/self.Mv450 452 return info[0], info[1], z, L_result/self.Mv, sum(acceptRates)/self.Mv 453 projects/AsynCluster/trunk/svpmc/test/test_model.py
r175 r177 331 331 self.failUnless(percent < 40) 332 332 333 def test_zProposal (self):333 def test_zProposal_dist(self): 334 334 N = 1000 335 335 z = s.zeros((2,3)) … … 338 338 z2 = s.empty(N) 339 339 for k in xrange(N): 340 z = modelObj.zProposal(z, 0.5)340 z, L_result = modelObj.zProposal(z, 0.8) 341 341 z1[k] = z[0,0] 342 342 z2[k] = z[1,1] 343 self.failUnlessNormal(z1 [0.7*N:])344 self.failUnlessNormal(z2 [0.7*N:])343 self.failUnlessNormal(z1) 344 self.failUnlessNormal(z2) 345 345 self._check_uncorr(z1, z2) 346 346 347 def test_zProposal_L(self): 348 N = 1000 349 z0 = s.ones((1,1)) 350 modelObj = self.modelFactory() 351 z = s.empty(N) 352 L = s.empty(N) 353 for k in xrange(N): 354 zk, L_result = modelObj.zProposal(z0, 0.8) 355 z[k] = zk[0,0] 356 L[k] = L_result 357 fig = self.fig 358 sp = fig.add_subplot(211) 359 sp.hist(z) 360 sp = fig.add_subplot(212) 361 sp.plot(z, s.exp(L), '.') 362 347 363 def test_decaySamples(self): 348 364 x = s.zeros((5, 1000)) projects/AsynCluster/trunk/svpmc/test/util.py
r168 r177 268 268 self.failUnlessAlmostEqual(error, 0.0, places) 269 269 270 def failUnlessNormal(self, x ):270 def failUnlessNormal(self, x, sdev=None): 271 271 p = stats.normaltest(x)[1] 272 272 if Mock.verbose(): … … 278 278 self.failUnless( 279 279 p > 0.05, "Deviation from normality with significance p=%f" % p) 280 if sdev: 281 if Mock.verbose(): 282 print "sdev = %f, expected %f" % (x.std(), sdev) 283 df = len(x) - 1 284 chi2 = df * x.var() / sdev**2 285 p = stats.chisqprob(chi2, df) 286 self.failUnless( 287 p > 0.05, 288 "Variance greater than expected with significance p=%f" % p) 280 289 281 290 @defer.deferredGenerator
