Changeset 193
- Timestamp:
- 05/23/08 13:59:26 (8 months ago)
- Files:
-
- projects/AsynCluster/trunk/svpmc/params.py (modified) (9 diffs)
- projects/AsynCluster/trunk/svpmc/project.py (modified) (1 diff)
- projects/AsynCluster/trunk/svpmc/test/test_params.py (modified) (1 diff)
Legend:
- Unmodified
- Added
- Removed
- Modified
- Copied
- Moved
projects/AsynCluster/trunk/svpmc/params.py
r192 r193 241 241 The exception is that you can use 'constant' as a distribution name. Then 242 242 the only pertinent parameter is the constant, unvarying I{loc}. 243 """ 243 244 @ivar V: A 1-D array of all the jump deviations that are used in 245 computing probability densities of jumps 246 247 """ 248 keyAttrs = {'V':None} 244 249 paramNames = ['dname', 'loc', 'scale', 'a', 'b'] 250 N_draw = 1000 251 N_attempts = 30 245 252 246 253 class Constant(object): … … 253 260 def ppf(self, *args): 254 261 return 0.0 262 263 class Jump(object): 264 __slots__ = ('x', 'px', 'pV') 265 def __init__(self, x, px, pV): 266 self.x = x 267 self.px = px 268 self.pV = pV 255 269 256 270 def _get_distObj(self): … … 280 294 return self.cache['rds'] 281 295 rdsObj = property(_get_rdsObj) 296 297 def _newDraw(self): 298 rJumps = self.rdsObj.rvs(self.N_draw) 299 V_inv = 1.0 / self.V 300 pJumps = self.rdsObj.pdf( 301 s.kron(rJumps, V_inv).reshape(self.N_draw, len(V_inv))) 302 pJumps *= V_inv 303 return {'k':-1, 'N':self.N_draw, 'R':rJumps, 'P':pJumps} 304 305 def _rJumpDraw(self): 306 if 'draw' not in self.cache: 307 draw = self.cache['draw'] = self._newDraw() 308 else: 309 draw = self.cache['draw'] 310 draw['k'] += 1 311 if draw['k'] >= draw['N']: 312 draw = self.cache['draw'] = self._newDraw() 313 k = draw['k'] 314 return draw['R'][k], draw['P'][k,:] 282 315 283 316 def pdf(self, x): … … 293 326 return self.distObj.rvs(1)[0] 294 327 295 def rJump(self, wiggleIndex, wiggles): 296 """ 297 Returns a value that jumps from zero by a random amount that, before 298 scaling by the supplied I{wiggle}, is drawn from a normal distribution 299 with a pre-scaled standard deviation that is set to approximate the 70% 300 probability width of my distribution. (That is the +/- 1 standard 301 deviation width of a normal distribution.) 302 """ 303 return wiggle * self.rdsObj.rvs(1)[0] 304 305 def pJump(self, x, wiggle): 306 """ 307 Returns the probability density of the supplied jump value I{x}, given 308 that it was generated by a call to my L{rJump} method with the 309 specified I{wiggle}. 310 311 For my Gaussian-PDF deviate generator, the probability density of the 328 def rJump(self, startValue, vIndex): 329 """ 330 Call this method with a I{startValue} and the index I{vIndex} that 331 points to the value in my array I{V} for the jump deviation to use for 332 the jump. 333 334 Draws a jump value from a normal distribution with a standard deviation 335 that is set to approximate the 70% probability width of my prior 336 distribution, scaled by the jump deviation. (The 70% width is the +/- 1 337 standard deviation width of a normal distribution.) 338 339 Returns a tiny object with two attributes, conveniently accessible as 340 part of a L{FlexArray}: 341 342 - B{x}: the value after the jump, and 343 344 - B{pV}: a 1-D array of probabilities of the jump for each jump 345 deviation in my array I{V}. 346 347 For my Gaussian-PDF deviate generator, the probability density of a 312 348 jump is inversely proportional to the scaling to its standard deviation 313 imparted by I{wiggle}. Thus, the unit-normal PDF is divided by the 314 I{wiggle} in the result. 315 """ 316 return self.rdsObj.pdf(x/wiggle) / wiggle 349 imparted by the jump deviation. Thus, the unit-normal PDF is divided by 350 each possible jump deviation in the I{p} array of the result. 351 """ 352 for k in xrange(self.N_attempts): 353 jumpValue, pV = self._rJumpDraw() 354 jumpValue *= self.V[vIndex] 355 x = startValue + jumpValue 356 px = self.pdf(x) 357 if px > 0: 358 break 359 else: 360 # Could't come up with a jump from this point, error 361 raise ValueError("Couldn't come up with a valid jump") 362 # TODO: fall back to an independent draw from the prior 363 # distribution 364 return self.Jump(x, px, pV) 317 365 318 366 … … 332 380 333 381 """ 334 attempts = 30335 382 typeChecking = True 336 383 … … 363 410 return paramContainer 364 411 365 def proposal(self, paramContainer, vIndex, pV):412 def proposal(self, paramContainer, vIndex, vWeights): 366 413 """ 367 414 Returns a new instance of L{ParameterContainer} with a random-walk 368 415 proposal based on the supplied I{paramContainer}, my priors, the jump 369 416 deviation defined at the I{vIndex} position of my 1-D array of jump 370 deviations, and a 1-D array I{ pV} of mixture weights for the jump417 deviations, and a 1-D array I{vWeights} of mixture weights for the jump 371 418 deviations in the jump distribution. 372 419 … … 379 426 normals has a standard deviation defined by one element of my jump 380 427 deviation array with a corresponding mixture weight in the supplied 381 I{ pV} array.428 I{vWeights} array. 382 429 """ 383 430 if vIndex < 0 or vIndex >= len(self.V): 384 431 raise ValueError("Invalid jump deviation index %d" % vIndex) 385 if getattr( pV, 'shape', 0) != self.V.shape:432 if getattr(vWeights, 'shape', 0) != self.V.shape: 386 433 raise ValueError( 387 434 "Must supply a mixture weight array with one element per "+\ … … 392 439 "You must supply a ParameterContainer as the first "+\ 393 440 "argument, not a '%s'" % type(paramContainer)) 394 vJump = self.V[vIndex]395 441 newVersion = ParameterContainer( 396 442 paramNames=self.paramNames, z=s.array(paramContainer.z).copy()) … … 400 446 # the weight for each jump deviation, just start the array with the 401 447 # weights. 402 pJumps = pV.copy()448 pJumps = vWeights.copy() 403 449 for name in self.paramNames: 404 450 priorsFA = getattr(self, name) 405 jumpsFA = priorsFA.call('rJump', vIndex, pV) 406 paramArray = getattr(paramContainer, name) + jumpsFA.x 407 pPriors = priorFlexArray.call('pdf', paramArray) 408 Lp += s.log(pPriors).sum() 451 startValues = getattr(paramContainer, name) 452 jumpsFA = priorsFA.call('rJump', startValues, vIndex) 453 Lp += s.log(jumpsFA.px).sum() 409 454 for k in xrange(len(self.V)): 410 pJumps[k] *= jumpsFA.p [k].prod()411 setattr(newVersion, name, paramArray)455 pJumps[k] *= jumpsFA.pV.call('item', k).prod() 456 setattr(newVersion, name, jumpsFA.x) 412 457 newVersion.Lp = Lp 413 458 newVersion.Lj = s.log(pJumps.sum()) projects/AsynCluster/trunk/svpmc/project.py
r174 r193 183 183 def constructPriors(priorSpec, shape): 184 184 index = priorSpec[0] 185 kw = {'dname':priorSpec[1] }185 kw = {'dname':priorSpec[1], 'V':self.V} 186 186 for k, name in enumerate(('loc', 'scale', 'a', 'b')): 187 187 stringValue = priorSpec[k+2] projects/AsynCluster/trunk/svpmc/test/test_params.py
r176 r193 232 232 self.priorContainer.typeChecking = False 233 233 self.priorContainer.paramNames = self.MPC.paramNames 234 kw = {'dname':'norm', 'loc':0 }234 kw = {'dname':'norm', 'loc':0, 'V':self.V} 235 235 for i, name in enumerate(self.priorContainer.paramNames): 236 236 shape = self.MPC.paramShapes[i]
