Changeset 186

Timestamp:

05/20/08 15:42:21 (8 months ago)

Author:

edsuom

Message:

Working on putting the likelihood method together, finally

Files:

• projects/AsynCluster/trunk/svpmc/model.c (modified) (16 diffs)
• projects/AsynCluster/trunk/svpmc/model.py (modified) (3 diffs)
• projects/AsynCluster/trunk/svpmc/test/test_model.py (modified) (8 diffs)

projects/AsynCluster/trunk/svpmc/model.c

@@ -122 +122 @@
-
-
-al sample
+ed log-volatility simulation
-struct selection
-{
-  int kp;
-p
+h
-};
-
@@ -213 +213 @@
-  // preparation for decorrelation
-  for(k=0; k<N; k++)
+    // TODO: Should we use a stripped-down version of exp here for speed?
-    SPA1(S, x0, k) = exp(-0.5 * SPA1(S, h, k)) * PA1(x, k);
-    
@@ -247 +248 @@
-// Returns with zero only if the post-first proposals are close enough to
-// finish likelihood computation
-*
+
-{
-  int k;
@@ -254 +255 @@
-  double Lmax = -HUGE_VAL;   // Anything will be more positive
-  for(k=0; k<n; k++) {
-->
+.
-    if(L < Lmin)
-      Lmin = L;
@@ -271 +272 @@
-// the log-probability of the proposal's having been selected.
-struct selection \
-*
+
-{
-  int k, kpMax;
-  double val;
-pSum = 0
+Dp[n]
-  double pMax = -HUGE_VAL;  // Anything will be more positive
-  double Dp[n];
-  struct selection result;
-
-
+ and their maximum
-  for(k=0; k<n; k++) {
-->
+.
-    Dp[k] = val;
-    pSum += val;
-    if(val > pMax) {
-      pMax = val;
@@ -299 +298 @@
-  } while(result.kp != kpMax && uniform(0, pMax) > Dp[result.kp]);
-
-
-
+
+
+
-  return result;
-}
@@ -342 +342 @@
-double val;
-double Lx = 0;
-p
+h
-
-// Multivariate innovation for one current time-series sample
@@ -348 +348 @@
-// Log probability density of each proposal on z
-double Lz[n];
-// Linear probability density of each log-volatility proposal
-double Dp[n];
-
-// A struct for the result of the importance-sampling function
@@ -373 +371 @@
-  for(k1=0; k1<2; k1++) {
-    sp[k0][k1].Lx = 0;
-
-
-
-
+
+
+
+
+
+
+
+
-}
-
@@ -389 +391 @@
-  for(k2=0; k2<n; k2++)
-    Lz[k2] = 0;
-
+
-  do {
-
@@ -405 +407 @@
-          Lz[k2] += normLogDensity(val);
-        }
-SPA1(sp[k2][notFirst], 
+PA1(sp[k2][notFirst].
-      }
-    }
@@ -426 +428 @@
-      for(k2=0; k2<n; k2++) {
-        for(k3=0; k3<Nt; k3++)
-SPA1(sp[k2][1], h, k3) = SPA1(sp[k2][0], 
+PA1(sp[k2][1].h, k3) = PA1(sp[k2][0].
-      }
-    }
@@ -432 +434 @@
-  } while (k1==k0+1 || (k1<Ns && notCloseEnough(sp, n)));
-
-
-
-
+
+
+
+
+
+
+
-  for(k1=0; k1<Nt; k1++)
-SPA1(sp[spSelection.kp][1], 
-
+PA1(sp[spSelection.kp][0].
+
-  // Set the initial log-volatility of the proposals for the next time-series
-  // sample to the log-volatility of the selected proposal for this time-series
-  // sample.
-n
-SPA1(sp[spSelection.kp][1], 
-Nt
-SPA1(sp[k1][0], h, k2
+Nt
+PA1(sp[spSelection.kp][1].
+n
+PA1(sp[k2][0].h, k1
-  }
-
@@ -453 +459 @@
-// sample
-for(k0=0; k0<Nt; k0++)
-SPA1(sp[kp][1], 
+PA1(sp[kp][1].
-
-// Free array memory
@@ -459 +465 @@
-for(k0=0; k0<n; k0++) {
-  for(k1=0; k1<2; k1++) {
+    free(sp[k0][k1].z);
-    free(sp[k0][k1].h);
-    free(sp[k0][k1].x0);
@@ -465 +472 @@
-}
-
-
-
-
-
-
+
+
+
+
+
+
+
-return_val = result;

projects/AsynCluster/trunk/svpmc/model.py

@@ -276 +276 @@
-        # ...including the latent parameter of IID normal variates
-        z = paramContainer.z.copy()
+        print self.fs, self.fr
-
-        # Set up arary variables for the inline call, including innovations as
@@ -293 +294 @@
-            return
-        # ...scaling constants for multivariate normal pdf
+        sc = s.empty((self.p, 2))
-        sc[:,0] = s.sqrt(1 - self.g**2)
-        sc[:,1] = s.log(s.sqrt(2*s.pi) * sc[:,0])
@@ -299 +301 @@
-        # one along with the state of the IID variate vector at that point.
-        Lx = 0
-
-
+
+
+
-            result = self.inline(
-                'z', 'h', 'x', 'w', 'n', 'd', 'e', 'g', 'rz', 'ri', 'sc',
-                w=self.wiggle, n=self.N1)
+            print "\n", k, result
+            print z[0,:]
-            Lx += result[0]
-p
-p
+h
+h

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

@@ -27 +27 @@
-from twisted.spread import pb
-from asynqueue import ThreadQueue
-
-
-
+
+
+
+
-import util
-
@@ -40 +41 @@
-
-
-params.
+
-    drift = 1.5
-    p, n = 8, 937
@@ -109 +110 @@
-            self.nextCall(Mock_Client.calls)
-            self.nextCall(
-params.
+
-            self.nextCall(
-                'update', ('newModel', self.model))
@@ -144 +145 @@
-        }
-
-    def modelFactory(self, **kw):
-        """
-        Instantiate a model to be tested and set its parameters
-        """
-        kw['paramNames'] = util.PARAM_NAMES
-        for name, value in self.params.iteritems():
-            if name != 'y':
-                if name in kw:
-                    value = kw[name]
-                kw[name] = s.array(value)
-        if 'y' not in kw:
-            if hasattr(self, 'y'):
-                kw['y'] = self.y
-            else:
-                kw['y'] = s.zeros((len(kw['fs']), self.N))
-        return model.Model(**kw)
-
-    def inline(self, code, **kw):
-        """
-        Do an inline call with the model's Weaver API and support code.
-        """
-self.modelFactory
+model.Model
-    
-    def _covarMatrix(self, cs, cr):
@@ -313 +297 @@
-                "Generated array \n%s\n != \n%s" % (y, x))
-
-self.modelFactory
+model.Model
-        tryCoeffs([1.0], [], [[1.0]])
-        tryCoeffs([1.0, 1.0], [0.2], [[1.0, 0.2], [0.2, 1.0]])
@@ -523 +507 @@
-
-class Test_Model_spArray(BaseCase):
-
+Base
-    int k;
-
-
+
-    struct sample sp[n][2];
+    for(k=0; k<n; k++)
+      sp[k][1].Lx = X1(k);
+    """
+    
+    codeNCE = codeBase + """
+    return_val = notCloseEnough(sp, n);
+    """
+    
+    codeIS = codeBase + """
-    double Lz[n];
-    py::tuple result(2);
-
- {
+struct selection spSelection;
+
-      Lz[k] = Z1(k);
-  sp[k][1].Lx = L1(k
-}
-
-S.Lx
+spSelection = importanceSample(sp, Lz, n
+result[0] = spSelection.kp;
+result[1] = spSelection.Lp;
+result
-    """
-
+    def _notCloseEnough(self, Lx_array):
+        return self.inline(self.codeNCE, X=Lx_array)
+
+    def _importanceSample(self, Lx_array, Lz_array=None):
+        if Lz_array is None:
+            Lz_array = s.zeros_like(Lx_array)
+        return self.inline(self.codeIS, X=Lx_array, Z=Lz_array)
+
+    def test_notCloseEnough_basic(self):
+        N = 8
+        Lx_array = s.zeros(N)
+        self.failUnlessEqual(self._notCloseEnough(Lx_array), 0)
+        Lx_array[0] = 1
+        self.failUnlessEqual(self._notCloseEnough(Lx_array), 1)
+
+    def test_importanceSample_basic(self):
+        N = 8
+        N_iter = 1000
+        Lx_array = s.zeros(N)
+        I = [self._importanceSample(Lx_array)[0] for k in xrange(N_iter)]
+        for k in xrange(N):
+            self.failUnless(k in I)
+
+    def test_importanceSample_normal(self):
+        N = 10
+        N_iter = 1000
+        distObj = stats.distributions.norm()
+        X = 8.0*s.rand(N_iter, N) - 4.0
+        Y = s.empty(N_iter)
+        Lp = s.empty(N_iter)
+        for k, Xk in enumerate(X):
+            Lx_array = s.log(distObj.pdf(Xk))
+            kp, Lpk = self._importanceSample(Lx_array)
+            Y[k] = Xk[kp]
+            Lp[k] = Lpk
+        if VERBOSE:
+            I = Y.argsort()
+            sp = self.fig.add_subplot(111)
+            sp.plot(Y, s.exp(Lp), '.', Y[I], distObj.pdf(Y[I]))
+        self.failUnlessNormal(Y)
+
-
-class Test_Model_likelihood(BaseCase):
-1000
+6
-    corr = 0.0
-    params = {
-        'a': [0],
-        'b': [[0]],
-        'cs': [1.0],
-        'cr': [],
-        'd': [0],
@@ -555 +587 @@
-        }
-
-
-
-
+
+
+
+
+
+
+
+
+
+
+
+
+
+
-        self.iv = random.multivariate_normal(
-            [0, 0],
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
-        if VERBOSE:
-            fig = self.fig
-[0,:], (self.h, h[0,:]
+, (self.h, h
-            for k, vector in enumerate(vectors):
-(len(vectors)+1
+len(vectors
-                sp = fig.add_subplot(spNumber)
-                if not isinstance(vector, (tuple, list)):
@@ -590 +632 @@
-                for v in vector:
-                    sp.plot(v)
-
-
+
-        
-
-
-
-
+
+
-        
-    def test_hybridGibbs_lowCorrelation(self):