Changeset 207

Timestamp:

06/02/08 19:04:49 (7 months ago)

Author:

edsuom

Message:

Got Rao-Blackwellized version working locally & remotely, with better but still disappointing degeneracy

Files:

• projects/AsynCluster/trunk/asyncluster/ndm/worker.py (modified) (1 diff)
• projects/AsynCluster/trunk/svpmc/model.c (modified) (3 diffs)
• projects/AsynCluster/trunk/svpmc/model.py (modified) (8 diffs)
• projects/AsynCluster/trunk/svpmc/params.py (modified) (9 diffs)
• projects/AsynCluster/trunk/svpmc/pmc.py (modified) (8 diffs)
• projects/AsynCluster/trunk/svpmc/project.py (modified) (6 diffs)
• projects/AsynCluster/trunk/svpmc/test/test_model.py (modified) (4 diffs)
• projects/AsynCluster/trunk/svpmc/test/test_params.py (modified) (9 diffs)
• projects/AsynCluster/trunk/svpmc/test/util.py (modified) (3 diffs)

projects/AsynCluster/trunk/asyncluster/ndm/worker.py

@@ -91 +91 @@
-    Manager()
-
+
+def runForeground(debug=False):
+    """
+    Runs a child worker L{Manager} in the foreground. Useful for debugging,
+    because you can see stdout & stderr.
+    """
+    if debug:
+        import pdb
+        try:
+            Manager()
+        except:
+            pdb.pm()
+    else:
+        Manager()

projects/AsynCluster/trunk/svpmc/model.c

@@ -182 +182 @@
-  double c, cp;
-  const double diff = b-a;
-  const double scaleUp = 1.06; // Increase scale value to account for EXP error
-  // Compute the scale value c as the maximum possible pdf (regular normal
-  // distribution) within the truncated range
-  if(a < 0) {
-    if(b > 0) {
-its pdf
-scaleUp;
+the maximum pdf at center
+1.02; // Higher than 1.0 due to max EXP +error of 1.02
-    } else {
-
-scaleUp * EXP(-0.5*b*b);
+maximum 
+1.0625 * EXP(-0.5*b*b); // Scaled up by 1.02 / 0.96
-    }
-  } else {
-    // Range is above center, so use pdf at its left edge
-scaleUp * EXP(-0.5*a*a);
+1.0625 * EXP(-0.5*a*a); // Scaled up by 1.02 / 0.96
-  }
-  // Accept-reject sampling for the truncated normal distribution
@@ -315 +314 @@
-// z    Independent normal variates, [pn] array (updated)
-// x    Innovation values, [pn] array
-// h    Simulated last-sample multivariate log-volatility, [p] (overwritten)
-// w    Wiggle value for +/- proposals (float)
-// ni   Number of hybrid-Gibbs iterations (int)
@@ -469 +467 @@
-}
-
-// Save the multivariate log-volatility simulated for the last time-series
-// sample in the last iteration
-for(k0=0; k0<Nt; k0++)
-  H1(k0) = PA1(se[ke].h, k0);
-
-// Free array memory
-free(xk);

projects/AsynCluster/trunk/svpmc/model.py

@@ -89 +89 @@
-            self._z[member,:,:] = z
-        if len(self.failures) + len(self.successes) == self.m:
-
+
+
-            k = 0
-            N = len(zPacked)
@@ -125 +126 @@
-        global MODEL
-        MODEL = modelObj
+        MODEL.z = None
-
-    def newZ(zPacked, k, N):
@@ -130 +132 @@
-        if k == 0:
-            CHUNKS = []
+            MODEL.z = None
-        CHUNKS.append(zPacked)
-
-
+
+
+
-    
-    def hybridGibbs(paramContainer):
-
+
+
+
-    
-    def likelihood(paramContainer):
+        if MODEL is None:
+            raise RuntimeError('No model available')
+        if MODEL.z is None:
+            raise RuntimeError(
+                'Model not initialized with simulated log-volatility data')
-        return MODEL.likelihood(paramContainer)
-
@@ -220 +231 @@
-        def unpackResult(packedResult):
-            if packedResult:
-list(pack.Unpacker(packedResult)
-
+pack.Unpacker(packedResult)(
+
-        def simulationDone(result):
-            if result is None:
-                update()
-                return False
-
-
+
-            return True
-        
@@ -257 +267 @@
-            return paramContainer
-
+        # DEBUG
-        if (remote or self.remoteMode) and not local:
-            d = self.client.run(
@@ -332 +343 @@
-    #--- Support --------------------------------------------------------------
-
-    def covarMatrix(self, cs, cr):
-        """
-        Returns a covariance matrix from the vector or sequence of I{cs}
-        standard deviations and I{cr} cross-correlations.
-
-        The I{cs} argument is in the form s0, s1,..., sp. It must have C{p}
-        elements.
-
-        The I{cr} argument must have C{0.5*(p**2 + p)} elements and is in the
-        form r00, r01, r02,..., r0p, r11, r12,..., r1p,..., rpp
-        """
-        i = 0
-        p = self.p
-        cv = s.zeros((p, p))
-        for j in xrange(p):
-            cv[j,j] = cs[j]**2
-            for k in xrange(j+1, p):
-                cv[j, k] = cv[k, j] = cr[i] * cs[j] * cs[k]
-                i += 1
-        return cv
-
-    def setParams(self, paramContainer):
-        """
-
-
+
+
+
+
-        self.setParamSequence(paramContainer.paramSequence())
-
-        # Set up array variables for the inline call, including innovations as
-        # residuals of my observations run through the reversed VAR(1)...
-        self.x = self.var.reverse(self.a, self.b)
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
+
-    
-    #--- The API --------------------------------------------------------------
@@ -396 +371 @@
-        If a L{linalg.LinAlgError} is raised due to an invalid correlation
-        matrix parameter, the method returns with no result. Otherwise, it
-
-
-
-
-
-
-
-
+
+
-
-        """
@@ -410 +379 @@
-        # every log-volatility simulation
-        z = s.zeros_like(self.x)
-        # This will contain the simulated log-volatility of the last
-        # time-series sample
-        h = s.empty(self.p)
-        self.inline(
-h', '
+
-            'd', 'e', 'g', 'rz', 'ri', 'sc',
-            w=self.wiggle, ni=self.Ni)
-, h
+
-
-    def likelihood(self, paramContainer):

projects/AsynCluster/trunk/svpmc/params.py

@@ -22 +22 @@
-import os.path
-import scipy as s
+from scipy import linalg
-from scipy.stats import distributions
-
@@ -248 +249 @@
-        - A scalar I{Lj} with the log-probability density of the proposal jump
-          resulting in my parameters.
-
-
-
-
-
-
+
+
+
-
-
@@ -268 +266 @@
-
-    """
-    keyAttrs = {'V':None}
-    paramNames = ['dname', 'loc', 'scale', 'a', 'b']
-    N_draw = 1000
-    N_attempts = 200
-
-    class Constant(object):
@@ -318 +314 @@
-
-    def _newDraw(self):
-        D = len(self.V)
-        rJumps = self.rdsObj.rvs(self.N_draw)
-        pJumps = self.rdsObj.pdf(rJumps)
@@ -348 +343 @@
-        return self.distObj.rvs(1)[0]
-
-
-
-
-
-
+
+
+
+
-
-        Draws a jump value from a normal distribution with a standard deviation
@@ -375 +369 @@
-        each possible jump deviation in the I{p} array of the result.
-        """
-        v = self.V[vIndex]
-        vInverse = 1.0 / v
-self.
+
-            rJump, pJump = self._rJumpDraw()
-            rJump *= v
@@ -404 +397 @@
-    @ivar p: The number of time series.
-
-
+
+
+
-    
-    """
-    typeChecking = True
+    N_attempts = 200
-    
-    def __init__(self, p, n):
@@ -420 +416 @@
-        return array
-
+    def covarMatrix(self, cs, cr):
+        """
+        Returns a covariance matrix from the vector or sequence of I{cs}
+        standard deviations and I{cr} cross-correlations.
+
+        The I{cs} argument is in the form s0, s1,..., sp. It must have C{p}
+        elements.
+
+        The I{cr} argument must have C{0.5*(p**2 + p)} elements and is in the
+        form r00, r01, r02,..., r0p, r11, r12,..., r1p,..., rpp
+        """
+        i = 0
+        p = len(cs)
+        cv = s.zeros((p, p))
+        for j in xrange(p):
+            cv[j,j] = cs[j]**2
+            for k in xrange(j+1, p):
+                cv[j, k] = cv[k, j] = cr[i] * cs[j] * cs[k]
+                i += 1
+        return cv
+    
+    def derivedParams(self, paramContainer):
+        """
+        Sets derived parameters in the supplied I{paramContainer} based on its
+        primary parameters, unless there is an error in generating the derived
+        parameters.
+
+        Returns C{True} if everything went OK, C{False} if not.
+        """
+        # Concurrent cross-correlations between volatility shocks and inverse
+        # of concurrent cross-correlations between innovation shocks
+        try:
+            rz = linalg.cholesky(self.covarMatrix(
+                paramContainer.fs, paramContainer.fr), lower=True)
+            ri = linalg.inv(linalg.cholesky(self.covarMatrix(
+                s.ones(self.p), paramContainer.cr), lower=True))
+        except linalg.LinAlgError:
+            return False
+        paramContainer.rz = rz
+        paramContainer.ri = ri
+        # Scaling constants for multivariate normal pdf
+        sc = s.empty((self.p, 2))
+        sc[:,0] = 1.0 / s.sqrt(1.0 - paramContainer.g**2)
+        sc[:,1] = s.log(s.sqrt(2*s.pi) / sc[:,0])
+        paramContainer.sc = sc
+        return True
+
-    def new(self):
-        """
-        Returns a new parameter container with values intialized from a random
-        draw of my priors.
-
-
-
-
-
-
-
-
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
-        paramContainer.Lp = Lp
-        paramContainer.Lj = 0 # There was no jump
-        return paramContainer
-
-Index, vWeights
+
-        """
-        Returns a new instance of L{ParameterContainer} with a random-walk
-
-
-
-
+
+
-
-        The returned parameter container object will have new values of I{Lp}
-        and I{Lj} corresponding to the proposal.
-
-
-
-
-
-
+
+
+
+
+
+
+
+
-        """
-        if self.typeChecking and \
@@ -458 +514 @@
-                "You must supply a ParameterContainer as the first "+\
-                "argument, not a '%s'" % type(paramContainer))
-
-
-
-
-
-
-
-
-
-
-
-
-
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
-        newVersion.Lp = Lp
-        newVersion.Lj = Lj

projects/AsynCluster/trunk/svpmc/pmc.py

@@ -157 +157 @@
-        self.pm = projectManager
-        self.mm = projectManager.mgr
+        self.V = projectManager.V
-        self.queue = projectManager.queue
-        self.resampler = Resampler(logWeights=True)            
@@ -171 +172 @@
-    def proposals(self, *args):
-        """
-n
-integer index of a particular jump deviation, and it will generate
-roposals from the existing p
+
+jump deviation, and it will generate proposals from the existing
+
-
-        Call it with just an integer number of proposals, and it will generate
-        new ones directly from the priors.
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
+
+
+
-        if len(args) == 2:
-            new = False
-Index
+
-            N = len(X)
-        else:
-            new = True
-            N = args[0]
-            vIndex = "+"
-        XP = params.FlexArray(N)
-        for k in xrange(N):
@@ -212 +194 @@
-                proposalPC = self.pm.priors.new()
-            else:
-
-
+
-            XP[k] = proposalPC
-
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
-            d.addCallback(simulationDone)
-            d.addErrback(self._oops, k=k)
-            dList.append(d)
-
-
-
+
+
+
+
+
+
-        """
-        Returns a deferred that fires with a 1-D array of log importance
-roposals supplied
+arameter proposals
-        
-        The importance weights are to be combined with those from other calls
@@ -241 +245 @@
-                    for x in ('Lx', 'Lp', 'Lj')])
-            else:
-                info = ""
-                L = -s.inf
-
+
+
+
+
-            return L
-        
@@ -250 +256 @@
-        dList = [self.mm.likelihood(XPk).addCallback(weight) for XPk in XP]
-        return defer.gatherResults(dList).addCallback(lambda W: s.array(W))
-
+
-    @defer.deferredGenerator
-    def run(self, N_iter):
@@ -266 +272 @@
-
-        """
-        def gotTheseWeights(Wv, I, Iv):
-            """
-            Returns a 1-D weights array of the same length as this allocation's
-            length, with infinitely negative values except for the valid
-            weights supplied in I{Wv}.
-            """
-            W = -s.inf * s.ones(len(I))
-            W[Iv] = Wv
-            return W
-
-        X = None
-        if self.socket:
@@ -286 +282 @@
-        for i in xrange(N_iter):
-            print "\n%03d\n%s" % (i+1, "-"*100)
-            self.mm.nextIteration()
-            if X is None:
-                print "Generating %d proposals from priors:" % N_members
-                # No existing members, generate proposals directly from the
-
-
+
+
+
+
+
-                yield wfd
-
-
-
-
-
-
-
+
+
+
+
+
+
+
+
+
+
+
+
+
-                for vIndex, I, d in self.allocator.assembler(resultList):
-
-
+
+
+
+
+
+
+
+
+
-                wfd = defer.waitForDeferred(defer.DeferredList(dList))
-                yield wfd
-                wfd.getResult()
-, validXP = resultList
+ = resultList[0]
-            print "\n"
-            # Evaluate the proposals
-            dList, resultList = [], []
-            for vIndex, I, d in self.allocator.assembler(resultList):
-
-
-
+
-                d.addErrback(self._oops, vIndex=vIndex)
-                dList.append(d)
-            # Record the previous iteration's members while the nodes work on
-            # evaluating the valid proposals for the current iteration...
-            if X is not None:
-                self.pm.writeParams(X)
-            # "Wait" here for all the proposals to be evaluated
-            wfd = defer.waitForDeferred(defer.DeferredList(dList))
-            yield wfd
@@ -332 +335 @@
-                self.allocator.updateAllocations(I0[I1])
-            else:
+                # No valid parameters, start over
+                X = None
-                self.allocator.updateAllocations()
-        wfd = defer.waitForDeferred(self.pm.done())

projects/AsynCluster/trunk/svpmc/project.py

@@ -70 +70 @@
-      iteration.
-
-parameter arrays in the parameter
-sequences used in the project
+arrays in the parameter containers
+used in the project for both primary and secondary parameters
-
-    @ivar D: The number of jump deviations used by the D-kernel PMC
@@ -91 +91 @@
-        'm':        'member',    
-        }
-    
-    def __init__(self, specFile, ncFile="svpmc.nc", m=1000, readOnly=False):
-        self.m = m
@@ -224 +223 @@
-                        pa[j,k] = params.Prior(**kw)
-        
-
+rimaryP
-        titles, dimensions = [], []
-        self.paramInfo = {}
@@ -236 +235 @@
-                constructPriors(priorSpec, shape)
-            titles.append(title)
-
-
+
+
+
+
+
+
+
-        return titles, dimensions
-
@@ -254 +258 @@
-        obs.long_name = "Observations, %d time series " % self.p +\
-                        "with %d samples each" % self.n
-        # Last-Sample Log-Volatility
-        logVol = cdf.createVariable(
-            "log_volatility", 'f', ('iteration', 'member', 'series'))
-        logVol.long_name = "Last Sample of Simulated Log-Volatilities"
-        # Parameters
-
+rimaryP
-            pDim = ['iteration', 'member'] + pDims[k]
-            param = cdf.createVariable(name, 'f', tuple(pDim))
-            param.long_name = pTitles[k]
-
-log_likelihood
+, log-density of priors and jumps
+Lx
-        L.long_name = "Log-Likelihood of Observations Given Parameters"
+        L = cdf.createVariable("Lp", 'f', ('iteration', 'member'))
+        L.long_name = "Log-Density of Parameters from Prior Distributions"
+        L = cdf.createVariable("Lj", 'f', ('iteration', 'member'))
+        L.long_name = "Log-Density of Proposal Jumps"
-        # Done setting up, save what we've got thus far
-        cdf.sync()
@@ -279 +283 @@
-                "Each iteration record must have %d parameters" % self.m)
-        # Write the parameters to their variables
-
+rimaryP
-            var = self.cdf.variables[name]
-            # Iterate over a 1-D FlexArray of the population's parameter arrays
-            for member, paramArray in enumerate(getattr(parameters, name)):
-                var[self.iteration, member, :] = paramArray.astype('f')
-ast sample's multivariate log-volatility value
+og-volatility, prior and jump log-density values
-        # population member
-
-
-
-
-
-
+
+
+
-        # Done writing this iteration of the CDF record
-        self.cdf.sync()

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

@@ -346 +346 @@
-
-class Test_Model_misc(BaseCase):
-    def test_covarMatrix(self):
-        def tryCoeffs(cs, cr, x):
-            x = s.array(x)
-            modelObj.y = s.empty((x.shape[0], 100))
-            y = modelObj.covarMatrix(cs, cr)
-            self.failUnlessEqual(y.shape, x.shape)
-            self.failUnless(
-                s.absolute(y-x).max() < 1E-8,
-                "Generated array \n%s\n != \n%s" % (y, x))
-
-        modelObj = model.Model()
-        tryCoeffs([1.0], [], [[1.0]])
-        tryCoeffs([1.0, 1.0], [0.2], [[1.0, 0.2], [0.2, 1.0]])
-        tryCoeffs([2, 3], [0], [[4.0, 0.0], [0.0, 9.0]])
-        tryCoeffs([3, 4], [0.5], [[9.0, 6.0], [6.0, 16.0]])
-        tryCoeffs(
-            [3, 4, 5], [0.1, 0.2, 0.3],
-            [[9.0, 1.2, 3.0], [1.2, 16.0, 6.0], [3.0, 6.0, 25.0]])
-
-    def test_support_uniform(self):
-        code = """
@@ -384 +365 @@
-          X1(k) = truncnorm(a, b);
-        """
-
+00
-        ranges = ((-1, +1), (-0.5, +1.5), (1.0, 2.0), (3.0, 5.0))
-        fig = self.fig
@@ -752 +733 @@
-            ('g', -0.9, +0.9)]
-        fig = self.fig
-1
+2
-        for j, paramInfo in enumerate(paramInfoList):
-            paramName = paramInfo[0]
@@ -761 +742 @@
-                Lx[k] = modelObj.inlineDirect(
-                    modelObj.code['likelihood'], **inlineVars)
-
+2*
-            sp.plot(values, Lx, '.')
-            sp.set_title("Parameter '%s'" % paramName)
+            sp = fig.add_subplot(spNumberBase+2*j+1)
+            Lx_max = Lx.max()
+            sp.plot(values, s.exp(Lx-Lx_max))
+            sp.set_title("Parameter '%s'" % paramName)
-            inlineVars[paramName] = originalParamValue
+
+    def test_Lx_vs_z(self):
+        wiggle = 1.0
+        eValue = 0.90
+        Ns, Ni, Nr = 1000, 40, 200
+
+        Lx = s.empty(Nr)
+        modelObj = model.Model()
+
+        inlineVars = self._inlineVars(wiggle, eValue, Ni)
+        for zReal, zSim, hReal, hSim in self._hybridGibbs(Ns, 1, inlineVars):
+            pass
+
+        fig = self.fig
+        zErrorList = [0.0005, 0.001, 0.01]
+        spNumberBase = 100*len(zErrorList) + 11
+        for j, zError in enumerate(zErrorList):
+            for k in xrange(Nr):
+                inlineVars['z'] = (
+                    zSim + zError*2*(s.rand(1, Ns) - 0.5)).reshape((1, 1, Ns))
+                Lx[k] = modelObj.inlineDirect(
+                    modelObj.code['likelihood'], **inlineVars)
+            sp = fig.add_subplot(spNumberBase+j)
+            sp.hist(Lx)
+            sp.set_title("Range (+/-) of zError: %f" % zError)
-    
-    def test_optimize_Ni(self):

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

@@ -222 +222 @@
-class Test_Prior(util.TestCase):
-    def setUp(self):
-
-
+
-        self.p.N_draw = 10
-        self.pdf = self.p.rdsObj.pdf
@@ -235 +234 @@
-        self.failUnless(s.equal(s.array([p.pdf(x) for x in X]), 0.5).all())
-
-v, pV
+Jump, pJump
-        pdf = self.p.rdsObj.pdf
-
-
+
-
-    def test_newDraw(self):
@@ -249 +247 @@
-    def test_rJumpDraw(self):
-        for k in xrange(103):
-v, pV
-v, pV
+Jump, pJump
+Jump, pJump
-
-    def test_rJump(self):
-        for startValue in s.linspace(-1, +3, 10):
-
-
-
+
+
-                jumpValue = jumpObj.x - startValue
-                rv = jumpValue / v
-
-
-
-
-
-
+
+
+
+
-    MPC = util.Mock_ParameterContainer
-    V = s.array([1.0, 0.1, 0.01])
-    
-    def setUp(self):
-        self.priorContainer = params.PriorContainer(
-, self.V
+
-        self.priorContainer.typeChecking = False
-
-
-
+
+