Changeset 190

Timestamp:

05/21/08 23:13:52 (7 months ago)

Author:

edsuom

Message:

Got C-based hybrid Gibbs working; Three-sample Lh evaluation appears optimal, to my surprise

Files:

• projects/AsynCluster/trunk/svpmc/model.c (modified) (4 diffs)
• projects/AsynCluster/trunk/svpmc/model.py (modified) (1 diff)
• projects/AsynCluster/trunk/svpmc/test/test_model.py (modified) (7 diffs)

projects/AsynCluster/trunk/svpmc/model.c

@@ -250 +250 @@
-
-
-// Returns with zero only if the proposals are close enough to finish
-// likelihood computation
-int notCloseEnough(struct sample sp[], int n)
-{
-  int k;
-  double L;
-  double Lmin = +HUGE_VAL;   // Anything will be less positive
-  double Lmax = -HUGE_VAL;   // Anything will be more positive
-  for(k=0; k<n; k++) {
-    L = sp[k].L_diff;
-    if(L < Lmin)
-      Lmin = L;
-    if(L > Lmax)
-      Lmax = L;
-  }
-  return (Lmax - Lmin > 0.001);
-}
-
-
-// Importance-samples one of the proposals based on a weight computed from the
-// log-density of z and the log-likelihood of the innovation given the
@@ -323 +304 @@
-// h    Simulated last-sample multivariate log-volatility, [p] (overwritten)
-// w    Wiggle value for +/- proposals (float)
+// ni   Number of hybrid-Gibbs iterations (int)
-// np   Number of proposals to try per log-volatility simulation (int)
-// ne   Max number of time-series samples to evaluate per proposal (int)
@@ -346 +328 @@
-const int Ns = Nz[1];
-
-
+i, k
-double val, L_diff;
-double Lx = 0;
@@ -383 +365 @@
-  sp[k0].x0 = (double *) malloc(sizeof(double) * Nt);
-  sp[k0].x1 = (double *) malloc(sizeof(double) * Nt);
+}
+
+
+// For each iteration...
+for(ki=0; ki<ni; ki++) {
+
-  // The "previous" log-volatility of the first time-series sample is set to
-  // the log-volatility offset.
-  //
-  // TODO: Treat it as another latent parameter and simulate its value
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
+
+
+
+
+
+
+
-    
-
-
-
-
-
-
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
-        if(k1==k0) {
-
-
-
-
+
+
+
-        }
-        PA1(sp[k2].z, k3) = val;
-      }
-
-
-
-
-
-
-
-
-
-
-
-
-
-
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
-    }
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
+
+
+
+
+
-}
-
-// Save the multivariate log-volatility simulated for the last time-series
-
+ in the last iteration
-for(k0=0; k0<Nt; k0++)
-  H1(k0) = PA1(sp[spSelection.kp].h, k0);

projects/AsynCluster/trunk/svpmc/model.py

@@ -299 +299 @@
-        # Run the hybrid Gibbs rounds, returning the likelihood from the last
-        # one along with the state of the IID variate vector at that point.
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
+
+
+
+
+
-        return Lx, Lh, z, h

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

@@ -595 +595 @@
-        return signal.lfilter([1], [1.0, -eValue], fsValue*z)
-
-1
+0
-        """
-        Runs the innovations through a VAR(1) to get simulated
@@ -611 +611 @@
-        return iv, h, x
-
-s, N
+
-        """
-        Returns a parameter container populated from my I{params} dict,
@@ -627 +627 @@
-                s.array(kw.get(name, self.params[name]), dtype='d'))
-        # Return the parameter container along with a model to be tested
-
+N1=N1, N2=N2, wiggle=wiggle
-        modelObj.paramNames = util.PARAM_NAMES
-        modelObj.N1=N1
-        modelObj.N2=N2
-        return pc, modelObj
-
@@ -644 +642 @@
-            sp.grid(True)
-    
-
+, params
-        """
-        Runs C code for L{model.Model.likelihood} independent of the method
@@ -650 +648 @@
-        """
-        z = s.zeros((1, Ns))
-
-
-
-
-
-
-
-
-
-
-
-
-
+
-        modelObj = model.Model()
-        for j in xrange(Nr):
-            iv, hReal, x = self._simData(Ns, eValue)
-
-
-
+
+
-            hSim = self._z_to_h(z[0,:], eValue, 1.0)
-            yield hReal, hSim
@@ -684 +669 @@
-            print "%3d: %+5.1f  %+6.1f%s%s" % (k, Lx[k], Lh[k], " "*5, zLine)
-
-
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
-        Ne = 2
-        wiggle = 3.0
@@ -692 +712 @@
-            self._plot((hReal, hSim))
-    
-    def test_highCorrelation_correlation(self):
-        # Empirically, it looks like Ne=2 is best (possibly Ne=3, but why?)
-        Ne = 2
-        Nr = 10000
-        wiggle = 3.0
-        eValue = 0.99
-        Ns, N1, N2 = 200, 20, 1
-        correlations = []
-        for hReal, hSim in self._likelihood(eValue, wiggle, Ns, Nr, N1, N2, Ne):
-            corrInfo = stats.spearmanr(hReal, hSim)
-            if VERBOSE:
-                print corrInfo[0]
-            correlations.append(corrInfo[0])
-            #self.failUnless(corrInfo[0] > 0.8)
-        sp = self.fig.add_subplot(111)
-        sp.hist(s.array(correlations))
-        sp.set_title("Ne=%d" % Ne)
-        #self._check(x, iv[1,:], zSim[-1,:], eValue)
-
-    def test_highCorrelation_model(self):
-        wiggle = 3.0
-        eValue = 0.95
-300, 20
+1000, 2
-        iv, h, x = self._simData(Ns, eValue)
-
-
+
-        pc.z = s.zeros((1, Ns))
-        modelObj.y = self._x_to_y(x)