Changeset 196

Timestamp:

05/24/08 00:05:18 (8 months ago)

Author:

edsuom

Message:

Huge progress, got PMC running with noticably less variance in population likelihoods

Files:

• projects/AsynCluster/trunk/doc/svpmc/example/svpmc.conf (modified) (1 diff)
• projects/AsynCluster/trunk/svpmc/model.c (modified) (2 diffs)
• projects/AsynCluster/trunk/svpmc/model.py (modified) (2 diffs)
• projects/AsynCluster/trunk/svpmc/params.py (modified) (7 diffs)
• projects/AsynCluster/trunk/svpmc/project.py (modified) (4 diffs)
• projects/AsynCluster/trunk/svpmc/test/test_params.py (modified) (2 diffs)

projects/AsynCluster/trunk/doc/svpmc/example/svpmc.conf

@@ -25 +25 @@
-
-
-
+
--------------------------------------------------------------------------------
-
-
+
+
+
+
+
-
-

projects/AsynCluster/trunk/svpmc/model.c

@@ -323 +323 @@
-
-// The number of time series making up a single multivariate sample
-z
+x
-// The number of multivariate samples
-z
+x
-
-int ki, k0, k1, k2, k3;
@@ -371 +371 @@
-
-  // The "previous" log-volatility of the first time-series sample is set to
-
-
-
+
+
-  for(k0=0; k0<np; k0++) {
-    for(k1=0; k1<Nt; k1++)

projects/AsynCluster/trunk/svpmc/model.py

@@ -64 +64 @@
-    ###########################################################################
-    """
-
-
+
+
+
+
-        self.queue = projectManager.queue
-
@@ -307 +309 @@
-        sc[:,1] = s.log(s.sqrt(2*s.pi) * sc[:,0])
-
+        print self.wiggle, self.N1, self.N2, self.Ne
-        # Run the hybrid Gibbs rounds, returning the likelihood from the last
-        # one along with the state of the IID variate vector at that point.

projects/AsynCluster/trunk/svpmc/params.py

@@ -268 +268 @@
-    the only pertinent parameter is the constant, unvarying I{loc}.
-
-    @ivar V: A 1-D array  of all the jump deviations that are used in
-      computing probability densities of jumps
-
-    """
-    keyAttrs = {'V':None}
@@ -322 +319 @@
-
-    def _newDraw(self):
+        D = len(self.V)
-        rJumps = self.rdsObj.rvs(self.N_draw)
-
-
-
-
+
+
+
+
+
+
+
+
+
-        return {'k':-1, 'N':self.N_draw, 'R':rJumps, 'P':pJumps}
-    
+
-    def _rJumpDraw(self):
-        if 'draw' not in self.cache:
@@ -376 +379 @@
-        each possible jump deviation in the I{p} array of the result.
-        """
+        D = len(self.V)
-        for k in xrange(self.N_attempts):
-pV
+P
-            jumpValue *= self.V[vIndex]
+            pV = P[D-vIndex-1:2*D-vIndex-1] / self.V
-            x = startValue + jumpValue
-            px = self.pdf(x)
@@ -384 +389 @@
-                break
-        else:
+            print "Oops"
-            # Could't come up with a jump from this point, fall back to an
-            # independent draw from the prior distribution
-            x = self.rvs()
-            px = self.pdf(x)
-_like(self.V
+(D
-        return self.Jump(x, px, pV)
-
@@ -402 +408 @@
-
-    @ivar paramNames: A list of the parameter names in sequence.
-
-    @ivar V: A 1-D array  of all the jump deviations that are used in
-      computing probability densities of jumps made via my prior objects.
-    
-    """
-    typeChecking = True
-    
-, V
+
-        self.p = p
-        self.n = n
-        self.V = V
-    
-    def priorArray(self, paramName, *shape):
@@ -455 +457 @@
-        I{vWeights} array.
-        """
-        if vIndex < 0 or vIndex >= len(self.V):
-            raise ValueError("Invalid jump deviation index %d" % vIndex)
-        if getattr(vWeights, 'shape', 0) != self.V.shape:
-            raise ValueError(
-                "Must supply a mixture weight array with one element per "+\
-                "jump deviation")
-        if self.typeChecking and \
-               not isinstance(paramContainer, ParameterContainer):
@@ -479 +475 @@
-            jumpsFA = priorsFA.call('rJump', startValues, vIndex)
-            Lp += s.log(jumpsFA.px).sum()
-self.V
+vWeights
-                pJumps[k] *= jumpsFA.pV.call('item', k).prod()
-            setattr(newVersion, name, jumpsFA.x)

projects/AsynCluster/trunk/svpmc/project.py

@@ -67 +67 @@
-      sequences used in the project.
-
-V: A 1-D array of the
+D: The number of
-      inference engine that will be running the project.
+
+    @ivar Df: The value of a given jump deviation as a fraction of the previous
+      one (the first value is 1.0)
-
-    @ivar priors: An instance of L{params.PriorContainer} constructed in
@@ -87 +90 @@
-        self.iteration = 0
-        #--- Use the NullQueue for debugging, thread queue normally -----------
-
-#
+#
+
-        #----------------------------------------------------------------------
-        specDir = os.path.dirname(specFile)
-        self.tables = self._parseSpec(specFile)
-
-
-
+
+
+
+
+
+
+
+
+
+
+
-        # Observations
-        tsData, seriesTitles = self._setupTimeSeries(specDir)
@@ -100 +111 @@
-        # Parameters
-        paramTitles, dimensions = self._setupParams()
-
+
+
-        self._setupCDF(
-            os.path.join(specDir, ncFileName),
@@ -206 +218 @@
-        titles, dimensions = [], []
-        xcorrs = self.xcorrs = int(0.5*(self.p**2 + self.p))
-, self.V
+
-        for name, dims, title in self.tables['parameter']:
-            shape = parseDims(dims)

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

@@ -221 +221 @@
-
-class Test_Prior(util.TestCase):
+    def setUp(self):
+        self.V = s.array([1.0, 0.1, 0.01])
+        self.p = params.Prior(dname='norm', loc=1, scale=1, D=3, Df=0.1)
+        self.p.N_draw = 10
+        self.pdf = self.p.rdsObj.pdf
+
-    def test_uniformPrior(self):
-        p = params.Prior(dname='uniform', loc=1, scale=2)
@@ -229 +235 @@
-        self.failUnless(s.equal(s.array([p.pdf(x) for x in X]), 0.5).all())
-
+    def _checkDraw(self, rv, pV):
+        pdf = self.p.rdsObj.pdf
+        for k, v in enumerate(self.V):
+            self.failUnlessAlmostEqual(pdf(rv/v)/v, pV[k], 5)
+
+    def test_newDraw(self):
+        draw = self.p._newDraw()
+        k = draw['k']
+        while k < draw['N']:
+            self._checkDraw(draw['R'][k], draw['P'][k])
+            k += 1
+
+    def test_rJumpDraw(self):
+        for k in xrange(103):
+            rv, pV = self.p._rJumpDraw()
+            self._checkDraw(rv, pV)
+
+    def test_rJump(self):
+        for startValue in s.linspace(-1, +3, 10):
+            for vIndex, v in enumerate(self.V):
+                jumpObj = self.p.rJump(startValue, vIndex)
+                print "\n", vIndex, v
+                jumpValue = jumpObj.x - startValue
+                rv = jumpValue / v
+                print jumpValue, rv, jumpObj.pV[vIndex]
+                self.failUnlessAlmostEqual(self.pdf(rv)/v, jumpObj.pV[vIndex])
+                print "OK"
+                
-
-class Test_PriorContainer(util.TestCase):