Changeset 176

Show
Ignore:
Timestamp:
05/09/08 22:13:35 (7 months ago)
Author:
edsuom
Message:

Redoing Rao-Blackwellization as the sum of the joint proposal-jump distributions

Files:

Legend:

Unmodified
Added
Removed
Modified
Copied
Moved

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

    r175 r176  
    123123                # An error from the remote likelihood method call is treated as 
    124124                # infinitely low likelihood 
    125                 result = (-s.inf, [], [], 0) 
     125                result = (-s.inf, [], [], 0, 0) 
    126126            elif isinstance(result, str): 
    127127                # Unpack string result from remote call 
    128128                result = list(pack.Unpacker(result)) 
    129             for k, name in enumerate(('Lx', 'h', 'z', 'acceptanceRate')): 
     129            for k, name in enumerate(('Lx', 'h', 'z', 'Lh', 'acceptanceRate')): 
    130130                setattr(paramContainer, name, result[k]) 
    131131            return paramContainer 
    …… 
    312312        Updates the IID normal variates in place. Returns the likelihood of the 
    313313        innovations in I{x}, given the simulated log-volatilities, along with 
    314         the last sample's multivariate log-volatility value and the fractional 
    315         acceptance rate of the MCMC steps. 
    316  
     314        the last sample's multivariate log-volatility value, the 
     315        log-probability of the simulated result, and the fractional acceptance 
     316        rate of the MCMC steps. 
     317         
    317318        The method will account for the effect of log-volatility proposals on 
    318319        downstream samples. You can specify the maximum fractional effect to 
    …… 
    321322        """ 
    322323        L_total = 0 
     324        L_result = 0 
    323325        acceptances = 0 
    324326        N = self.decaySamples(tol) 
    …… 
    343345            # Metropolis-Hastings accept/reject of the proposal 
    344346            dL = LVp[1] - LV[1] 
    345             if dL > 0 or self.logUniform() < dL: 
     347            if dL < -600 or not s.isfinite(dL): 
     348                # Bogus likelihood or infintesimally small probability of 
     349                # acceptance, just reject 
     350                accept = False 
     351            elif dL > 0: 
     352                # Always accept if the proposal's likelihood is better 
     353                accept = True 
     354            elif self.logUniform() < dL: 
     355                # Proposal likelihood worse, but accepted 
     356                L_result += dL 
     357                accept = True 
     358            else: 
     359                # Proposal likelihood worse and rejected 
     360                L_result += s.log(1 - s.exp(dL)) 
     361                accept = False 
     362            if accept: 
    346363                LV = LVp 
    347364                z[:,k] = zp[:,k] 
    …… 
    354371            # guaranteed not to be resampled, just bail out now and save time 
    355372            if L_total < -1E6: 
    356                 return -s.inf, s.zeros_like(h0), 0.0 
    357         return L_total, h0, float(acceptances)/self.n 
     373                return -s.inf, s.zeros_like(h0), 0.0, 0.0 
     374        return L_total, h0, L_result, float(acceptances)/self.n 
    358375     
    359376 
    …… 
    395412            3. The IID normal variates underlying the log-volatilities. 
    396413 
    397             4. The mean acceptance rate of the Metropolis-Hastings proposals. 
     414            4. The log-probability of the simulated result. 
     415 
     416            5. The mean acceptance rate of the Metropolis-Hastings proposals. 
    398417         
    399418        The returned likelihood does not consider the prior probability of the 
    …… 
    419438        # one along with the state of the IID variate vector at that point, 
    420439        # bailing out with no result if invalid parameter raises exception 
     440        L_result = 0 
    421441        acceptRates = [] 
    422442        try: 
    423             for k in xrange(self.Mv-1): 
    424                 acceptRates.append(self.hybridGibbs(z, x, rv, sigma)[2]) 
     443            for k in xrange(self.Mv): 
     444                info = self.hybridGibbs(z, x, rv, sigma) 
     445                L_result += info[2] 
     446                acceptRates.append(info[3]) 
    425447        except linalg.LinAlgError: 
    426448            return 
    427         L, h0, thisRate = self.hybridGibbs(z, x, rv, sigma) 
    428         return L, h0, z, sum(acceptRates + [thisRate])/self.Mv 
    429      
     449        return info[0], info[1], z, L_result, sum(acceptRates)/self.Mv 
     450     

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

    r174 r176  
    290290        return self.distObj.rvs(1)[0] 
    291291 
    292     def rJump(self): 
    293         """ 
    294         Returns a value that jumps from zero by a random amount that is drawn 
    295         from a normal distribution with a pre-scaled standard deviation that is 
    296         set to approximate the 70% probability width of my distribution. (That 
    297         is the +/- 1 standard deviation width of a normal distribution.) 
    298         """ 
    299         return self.rdsObj.rvs(1)[0] 
    300  
    301     def pJump(self, x, V, pV): 
    302         """ 
    303         Returns the probability density of the supplied jump value I{x} given a 
    304         mixture distribution of normals with standard deviations in the 
    305         supplied 1-D array I{V} of jump deviations and like-dimensioned array 
    306         I{pV} of mixture weights for normals having those deviations. 
    307         """ 
    308         densities = self.rdsObj.pdf(x/V) / V 
    309         return (pV * densities).sum() 
     292    def rJump(self, wiggle): 
     293        """ 
     294        Returns a value that jumps from zero by a random amount that, before 
     295        scaling by the supplied I{wiggle}, is drawn from a normal distribution 
     296        with a pre-scaled standard deviation that is set to approximate the 70% 
     297        probability width of my distribution. (That is the +/- 1 standard 
     298        deviation width of a normal distribution.) 
     299        """ 
     300        return wiggle * self.rdsObj.rvs(1)[0] 
     301 
     302    def pJump(self, x, wiggle): 
     303        """ 
     304        Returns the probability density of the supplied jump value I{x}, given 
     305        that it was generated by a call to my L{rJump} method with the 
     306        specified I{wiggle}. 
     307 
     308        For my Gaussian-PDF deviate generator, the probability density of the 
     309        jump is inversely proportional to the scaling to its standard deviation 
     310        imparted by I{wiggle}. Thus, the unit-normal PDF is divided by the 
     311        I{wiggle} in the result. 
     312        """ 
     313        return self.rdsObj.pdf(x/wiggle) / wiggle 
    310314 
    311315 
    …… 
    360364        Returns a new instance of L{ParameterContainer} with a random-walk 
    361365        proposal based on the supplied I{paramContainer}, my priors, the jump 
    362         deviation defined at the I{vIndex} position of my 1-D arrayt of jump 
     366        deviation defined at the I{vIndex} position of my 1-D array of jump 
    363367        deviations, and a 1-D array I{pV} of mixture weights for the jump 
    364368        deviations in the jump distribution. 
    …… 
    376380        if vIndex < 0 or vIndex >= len(self.V): 
    377381            raise ValueError("Invalid jump deviation index %d" % vIndex) 
    378         if len(pV) != len(self.V): 
     382        if getattr(pV, 'shape', 0) != self.V.shape: 
    379383            raise ValueError( 
    380                 "Mixture weight sequence must have one element per "+\ 
     384                "Must supply a mixture weight array with one element per "+\ 
    381385                "jump deviation") 
    382386        if self.typeChecking and \ 
    …… 
    388392        newVersion = ParameterContainer( 
    389393            paramNames=self.paramNames, z=s.array(paramContainer.z).copy()) 
    390         Lp, Lj = 0, 0 
     394        Lp = 0 
     395        # Define a 1-D array holding the probability density for each jump, 
     396        # joint across all parameters. Since the densities will be scaled by 
     397        # the weight for each jump deviation, just start the array with the 
     398        # weights. 
     399        pJumps = pV.copy() 
    391400        for name in self.paramNames: 
    392401            priorFlexArray = getattr(self, name) 
    393402            for attempt in xrange(self.attempts): 
    394                 jumps = vJump * priorFlexArray.call('rJump') 
     403                jumps = priorFlexArray.call('rJump', vJump) 
    395404                paramArray = getattr(paramContainer, name) + jumps 
    396405                pPriors = priorFlexArray.call('pdf', paramArray) 
    …… 
    401410                return paramContainer 
    402411            Lp += s.log(pPriors).sum() 
    403             pMixtures = priorFlexArray.call( 
    404                 'pJump', jumps, self.V, pV, arrayArgs=(1,2)) 
    405             Lj += s.log(pMixtures).sum() 
     412            for k, v in enumerate(self.V): 
     413                pJumps[k] *= priorFlexArray.call('pJump', jumps, v).prod() 
    406414            setattr(newVersion, name, paramArray) 
    407415        newVersion.Lp = Lp 
    408         newVersion.Lj = Lj 
     416        newVersion.Lj = s.log(pJumps.sum()) 
    409417        return newVersion 

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

    r175 r176  
    199199        def weight(paramContainer): 
    200200            if s.isfinite(paramContainer.Lx): 
    201                 L = paramContainer.Lx + paramContainer.Lp - paramContainer.Lj 
     201                L = paramContainer.Lx + paramContainer.Lp \ 
     202                    - paramContainer.Lj #- paramContainer.Lh 
    202203                info = " ".join([ 
    203204                    "%+12.2f" % (getattr(paramContainer, x),) 
    204                     for x in ("Lx", "Lp", "Lj")]) 
     205                    for x in ('Lx', 'Lp', 'Lj', 'Lh')]) 
    205206                info += "\t\t%02d%%" % (100*paramContainer.acceptanceRate,) 
    206207            else: 
    207208                info = "" 
    208209                L = -s.inf 
    209             print "%6.4f\t%s" % (self.pm.V[vIndex], info) 
     210            print "%6.4f\t%+12.2f\t%s" % (self.pm.V[vIndex], L, info) 
    210211            return L 
    211212         
    …… 
    256257              % (N_iter, N_members) 
    257258        for i in xrange(N_iter): 
    258             print "\n%03d\n%s" % (i+1, "-"*79) 
     259            print "\n%03d\n%s" % (i+1, "-"*100) 
    259260            dList, resultList = [], [] 
    260261            for vIndex, I, d in self.allocator.assembler(resultList): 

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

    r174 r176  
    290290    def test_proposal_basic(self): 
    291291        old = util.Mock_ParameterContainer(0.0) 
    292         new = self.priorContainer.proposal(old, 0, [1.0, 0.0, 0.0]) 
     292        new = self.priorContainer.proposal(old, 0, s.array([1.0, 0.0, 0.0])) 
     293        self.failUnless(s.isfinite(new.Lj)) 
     294        #print new.Lj, sum([ 
     295        #    (x**2).sum() 
     296        #    for x in [ 
     297        #    getattr(new, y) for y in self.priorContainer.paramNames]]) 
    293298        return self._check_basic(new) 
    294299     
    …… 
    300305        #X = self.priorContainer.new() 
    301306        for k in xrange(N): 
    302             XP = self.priorContainer.proposal(X, 1, [0.0, 1.0, 0.0]) 
     307            XP = self.priorContainer.proposal(X, 1, s.array([0.0, 1.0, 0.0])) 
    303308            if k == 0 or XP.Lp > X.Lp or s.log(s.rand()) < (XP.Lp - X.Lp): 
    304309                print "+",