Changeset 192

Show
Ignore:
Timestamp:
05/23/08 11:43:54 (8 months ago)
Author:
edsuom
Message:

Got C code working nicely and integrated with pmc; working on making proposals more efficient

Files:

Legend:

Unmodified
Added
Removed
Modified
Copied
Moved

  • projects/AsynCluster/trunk/svpmc/model.c

    r191 r192  
    257257  int k; 
    258258  double val; 
     259  double Dp_sum = 0; 
    259260  double Lp[n], Dp[n]; 
    260261  double L_max = -HUGE_VAL; // Anything will be more positive 
    …… 
    271272  // Compute linear probability density for each proposal, relative to the 
    272273  // maximum one 
    273   for(k=0; k<n; k++) 
    274     Dp[k] = exp(Lp[k] - L_max); 
     274  for(k=0; k<n; k++) { 
     275    val = exp(Lp[k] - L_max); 
     276    Dp[k] = val; 
     277    Dp_sum += val; 
     278  } 
    275279 
    276280  // Accept-reject sampling with uniform random proposal selection 
    …… 
    282286  } while(uniform(0, 1.0) > Dp[result.kp]); 
    283287 
    284   // The probability density of the selected log-volatility proposal is just 
    285   // the linear, unnormalized weight of the selected sample 
    286   result.Lp = Lp[result.kp]; 
     288  // The probability density of the selected log-volatility proposal 
     289  result.Lp = log(Dp[result.kp] / Dp_sum); 
    287290  return result; 
    288291} 

  • projects/AsynCluster/trunk/svpmc/model.py

    r190 r192  
    127127                # Unpack string result from remote call 
    128128                result = list(pack.Unpacker(result)) 
    129             for k, name in enumerate(('Lx', 'Lp', 'z', 'h')): 
     129            for k, name in enumerate(('Lx', 'Lh', 'z', 'h')): 
    130130                setattr(paramContainer, name, result[k]) 
    131131            return paramContainer 
    …… 
    178178 
    179179    @ivar wiggle: The +/- range of the random-walk uniform proposal on z for 
    180       each hybrid Gibbs proposal. 
     180      each hybrid Gibbs proposal. Use +/-1.0 for balance between fast 
     181      convergence and reasonable accept-reject efficiency in z proposal 
     182      generator. 
    181183 
    182184    @ivar N1: Number of log-volatility proposals to compute for the importance 
    183       sampling of each hybrid Gibbs step. 
    184  
    185     @ivar N2: Number of hybrid-Gibbs iterations per likelihood evaluation. 
    186  
    187     """ 
    188     keyAttrs = {'y':None, 'wiggle':0.3, 'N1': 8, 'N2':1} 
     185      sampling of each hybrid Gibbs step. For univariate, N1=10 appears most 
     186      cost-effective with wiggle=2.0. Scale for smaller wiggle and to account 
     187      for multivariate. 
     188 
     189    @ivar N2: Number of hybrid-Gibbs iterations per likelihood evaluation. For 
     190      univariate, N2=2 appears most cost-effective, but set to 4 anyays. Maybe 
     191      exploration of multivariate space needs more iterations. 
     192 
     193    @ivar Ne: The number of successive time-series samples to include in the 
     194      evaluation of each log-volatility proposal. For some weird reason, Ne=3 
     195      is optimal, at least for univariate. 
     196     
     197    """ 
     198    keyAttrs = {'y':None, 'wiggle':1.0, 'N1': 20, 'N2':6, 'Ne':3} 
    189199 
    190200    #--- Properties ----------------------------------------------------------- 
    …… 
    303313            'w', 'ni', 'np', 'ne', 
    304314            'd', 'e', 'g', 'rz', 'ri', 'sc', 
    305             w=self.wiggle, np=self.N1, ni=self.N2, ne=2) 
    306         return Lx, Lh, z, h 
     315            w=self.wiggle, np=self.N1, ni=self.N2, ne=self.Ne) 
     316        return Lx/self.N2, Lh/(self.N2*self.Ne), z, h 

  • projects/AsynCluster/trunk/svpmc/params.py

    r176 r192  
    224224          resulting in my parameters. 
    225225 
    226     """ 
    227     keyAttrs = {'z':[], 'Lx':None, 'Lp':None, 'Lj':None} 
     226        - A scalar I{Lh} with the log-probability density of the simulated 
     227          log-volatilities, reflected in the simulated values of I{z}. 
     228     
     229    """ 
     230    keyAttrs = {'z':[], 'Lx':None, 'Lp':None, 'Lj':None, 'Lh':None} 
    228231 
    229232 
    …… 
    290293        return self.distObj.rvs(1)[0] 
    291294 
    292     def rJump(self, wiggle): 
     295    def rJump(self, wiggleIndex, wiggles): 
    293296        """ 
    294297        Returns a value that jumps from zero by a random amount that, before 
    …… 
    399402        pJumps = pV.copy() 
    400403        for name in self.paramNames: 
    401             priorFlexArray = getattr(self, name) 
    402             for attempt in xrange(self.attempts): 
    403                 jumps = priorFlexArray.call('rJump', vJump) 
    404                 paramArray = getattr(paramContainer, name) + jumps 
    405                 pPriors = priorFlexArray.call('pdf', paramArray) 
    406                 if s.greater(pPriors, 0).all(): 
    407                     break 
    408             else: 
    409                 #print "\nWARNING: Failed to generate a valid proposal" 
    410                 return paramContainer 
     404            priorsFA = getattr(self, name) 
     405            jumpsFA = priorsFA.call('rJump', vIndex, pV) 
     406            paramArray = getattr(paramContainer, name) + jumpsFA.x 
     407            pPriors = priorFlexArray.call('pdf', paramArray) 
    411408            Lp += s.log(pPriors).sum() 
    412             for k, v in enumerate(self.V): 
    413                 pJumps[k] *= priorFlexArray.call('pJump', jumps, v).prod() 
     409            for k in xrange(len(self.V)): 
     410                pJumps[k] *= jumpsFA.p[k].prod() 
    414411            setattr(newVersion, name, paramArray) 
    415412        newVersion.Lp = Lp 

  • projects/AsynCluster/trunk/svpmc/pmc.py

    r178 r192  
    141141    """ 
    142142    chunkSize = 500 
    143     initialSigma = 0.3 
    144143     
    145144    def __init__(self, projectManager, socket=None): 
    …… 
    204203                    "%+12.2f" % (getattr(paramContainer, x),) 
    205204                    for x in ('Lx', 'Lp', 'Lj', 'Lh')]) 
    206                 info += "\t\t%02d%%" % (100*paramContainer.acceptanceRate,) 
    207205            else: 
    208206                info = "" 
    …… 
    214212        N = len(X) 
    215213        W = s.empty(len(X)) 
    216         # TODO: Adapt sigma over iterations for optimum acceptance rate in 
    217         # volatility simulations 
    218         sigma = self.initialSigma 
    219214        wfd = defer.waitForDeferred(self.queue.call(self.proposals, X, vIndex)) 
    220215        yield wfd 
    …… 
    224219            # Here's where the vast majority of the CPU time is expended! 
    225220            dList = [ 
    226                 self.mm.likelihood(XP[k], sigma).addCallback(weight) 
     221                self.mm.likelihood(XP[k]).addCallback(weight) 
    227222                for k in xrange(j, j+N_this)] 
    228223            wfd = defer.waitForDeferred(defer.gatherResults(dList)) 

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

    r191 r192  
    506506 
    507507 
    508 class Test_Model_spArray(BaseCase): 
    509     codeBase = """ 
     508class Test_Model_importanceSample(BaseCase): 
     509    code = """ 
    510510    int k; 
    511511    const int n = NX[0]; 
    512512    struct sample sp[n]; 
    513     for(k=0; k<n; k++) { 
    514       sp[k].L_diff = X1(k); 
    515       sp[k].Lx = X1(k); 
    516     } 
    517     """ 
    518      
    519     codeNCE = codeBase + """ 
    520     return_val = notCloseEnough(sp, n); 
    521     """ 
    522      
    523     codeIS = codeBase + """ 
    524513    double Lz[n]; 
    525514    py::tuple result(2); 
    526515    struct selection spSelection; 
    527     for(k=0; k<n; k++) 
     516 
     517    for(k=0; k<n; k++) { 
     518      sp[k].Lx = X1(k); 
    528519      Lz[k] = Z1(k); 
     520    } 
    529521    spSelection = importanceSample(sp, Lz, n); 
    530522    result[0] = spSelection.kp; 
    …… 
    533525    """ 
    534526 
    535     def _notCloseEnough(self, Lx_array): 
    536         return self.inline(self.codeNCE, X=Lx_array) 
    537  
    538527    def _importanceSample(self, Lx_array, Lz_array=None): 
    539528        if Lz_array is None: 
    540529            Lz_array = s.zeros_like(Lx_array) 
    541         return self.inline(self.codeIS, X=Lx_array, Z=Lz_array) 
    542  
    543     def test_notCloseEnough_basic(self): 
    544         N = 8 
    545         Lx_array = s.zeros(N) 
    546         self.failUnlessEqual(self._notCloseEnough(Lx_array), 0) 
    547         Lx_array[0] = 1 
    548         self.failUnlessEqual(self._notCloseEnough(Lx_array), 1) 
    549  
    550     def test_importanceSample_basic(self): 
     530        return self.inline(self.code, X=Lx_array, Z=Lz_array) 
     531 
     532    def test_basic(self): 
    551533        N = 8 
    552534        N_iter = 1000 
    …… 
    556538            self.failUnless(k in I) 
    557539 
    558     def test_importanceSample_normal(self): 
    559         N = 10 
     540    def _check(self, z0, wiggle): 
     541        N = 20 
    560542        N_iter = 1000 
    561         distObj = stats.distributions.norm() 
    562         X = 8.0*s.rand(N_iter, N) - 4.0 
     543        distObj = stats.distributions.truncnorm(z0-wiggle, z0+wiggle) 
     544        X = 2*wiggle*(s.rand(N_iter, N) - 0.5) + z0 
    563545        Y = s.empty(N_iter) 
    564546        Lp = s.empty(N_iter) 
    …… 
    571553            I = Y.argsort() 
    572554            sp = self.fig.add_subplot(111) 
    573             sp.plot(Y, s.exp(Lp), '.', Y[I], distObj.pdf(Y[I])) 
     555            sp.plot(Y, s.exp(Lp), '.', Y[I], s.pi/N*distObj.pdf(Y[I])) 
    574556        self.failUnlessNormal(Y) 
     557 
     558    def test_full(self): 
     559        return self._check(0, 4.0) 
     560 
     561    def test_closePortion(self): 
     562        return self._check(1.5, 1.0) 
     563 
     564    def test_farPortion(self): 
     565        return self._check(2.5, 0.5) 
    575566 
    576567 
    …… 
    712703        inlineVars = self._inlineVars(wiggle, eValue, N1, N2) 
    713704        self._simRounds(Nr, Ns, inlineVars, 'ne', (2,3,4,5,6)) 
    714         #self._check(x, iv[1,:], zSim[-1,:], eValue) 
    715705 
    716706    def test_optimize_N1(self): 
    …… 
    724714        inlineVars = self._inlineVars(wiggle, eValue, 1, N2) 
    725715        self._simRounds(Nr, Ns, inlineVars, 'np', (2,5,10,20)) 
    726         #self._check(x, iv[1,:], zSim[-1,:], eValue) 
    727716 
    728717    def test_optimize_N2(self): 
    …… 
    736725        inlineVars = self._inlineVars(wiggle, eValue, N1, 1) 
    737726        self._simRounds(Nr, Ns, inlineVars, 'ni', (1,2,3,4)) 
    738         #self._check(x, iv[1,:], zSim[-1,:], eValue) 
    739727     
    740728    def test_highCorrelation_c(self): 
    …… 
    748736     
    749737    def test_highCorrelation_model(self): 
    750         wiggle = 3.0 
     738        wiggle = 2.0 
    751739        eValue = 0.95 
    752         Ns, N1, N2 = 1000, 2, 1 
     740        Ns, N1, N2 = 1000, 10, 2 
    753741        iv, h, x = self._simData(Ns, eValue) 
    754742        pc, modelObj = self._modelAndParamFactory(wiggle, N1, N2, e=[[eValue]])