Changeset 191

Timestamp:

05/22/08 17:21:43 (7 months ago)

Author:

edsuom

Message:

Figured out N1=10, N2=2, Ne=3 sim params

Files:

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

projects/AsynCluster/trunk/svpmc/model.c

@@ -110 +110 @@
-  // of h to this point
-  double Lx;
-  // Change in Lx from this time-series sample to the previous one
-  double L_diff;
-  // Multivariate normal proposal, which is just z for time-series samples
-  // after the first in each log-volatility computation
@@ -227 +225 @@
-  
-  // Now the log-likelihoood. For each time series...
-  S->L_diff = 0;
-  for(k=0; k<N; k++) {
-    // Offset the decorrelated innovation by the mean, which is the correlation
@@ -240 +237 @@
-    val = SPA1(S, x0, k) / SPP2(P, sc, k, 0);
-    // Fast squaring instead of pow()
-_diff
+x
-    // Since the innovation has been inverse-scaled by the square root of the
-    // proposed volatility, the probability density needs to be scaled as
-    // well. In log space, this is a subtraction of half the log-volatility.
-
-
-
+
+
-}
-
@@ -329 +325 @@
-
-int ki, k0, k1, k2, k3;
-, L_diff
+
-double Lx = 0;
-double Lh = 0;
@@ -450 +446 @@
-      Z2(k1, k0) = zp[spSelection.kp];
-  
-    //--- DEBUG -----------------------------------------------------------------
-    //printf("\n%d %d\n---------------------------\n",k0, spSelection.kp);
-    //for(k1=0; k1<n; k1++)
-    //  printf("%d %f %f\n", k1, Lxk[k1], sp[k1].Lx);
-    //---------------------------------------------------------------------------
-  
-    // 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

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

@@ -611 +611 @@
-        return iv, h, x
-
+    def _inlineVars(self, wiggle, eValue, N1, N2=1, Ne=3):
+        return {
+            'h':s.empty(1),
+            'sc':s.array([1, s.log(s.sqrt(2*s.pi))]),
+            'w':wiggle,
+            'np':N1,
+            'ni':N2,
+            'ne':Ne,
+            'd':s.array([0.0]),
+            'e':s.array([[float(eValue)]]),
+            'g':s.array([0.0]),
+            'rz':s.array([[1.0]]),
+            'ri':s.array([[1.0]])
+            }
+    
-    def _modelAndParamFactory(self, wiggle, N1, N2, **kw):
-        """
@@ -642 +657 @@
-            sp.grid(True)
-    
-eValue, wiggle, Ns, Nr, N1, N2, Ne, param
+Ns, Nr, inlineVar
-        """
-        Runs C code for L{model.Model.likelihood} independent of the method
@@ -648 +663 @@
-        """
-        z = s.zeros((1, Ns))
-param
+inlineVar
-        modelObj = model.Model()
+        eValue = inlineVars['e'][0][0]
-        for j in xrange(Nr):
-            iv, hReal, x = self._simData(Ns, eValue)
-param
-param
+inlineVar
+inlineVar
-            hSim = self._z_to_h(z[0,:], eValue, 1.0)
-            yield hReal, hSim
@@ -661 +677 @@
-        eValue = 0.90
-        Ns, N1, N2 = 10, 5, 20
-
-
-
+
+
+
+
-        print "-"*100
-
-
-
-
-
-
-
-
-
-
-
+
+
+
+
+
+
+
-        fig = self.fig
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
+
+
+
-            correlations = []
-
-
+
-                corrInfo = stats.spearmanr(hReal, hSim)
-                if VERBOSE:
-                    print corrInfo[0]
-                correlations.append(corrInfo[0])
-
-
-
+
+
+
+
+
+
+
+
+
+
+
+
+
-        #self._check(x, iv[1,:], zSim[-1,:], eValue)
+
+    def test_optimize_N1(self):
+        # Empirically, it looks like np=10 is best. Going to np=20 barely gives
+        # any improvement, and obviously doubles the amount of computation
+        # time. Dropping to np=5 is significantly worse.
+        Nr = 1000
+        wiggle = 3.0
+        eValue = 0.95
+        Ns, N2 = 1000, 2
+        inlineVars = self._inlineVars(wiggle, eValue, 1, N2)
+        self._simRounds(Nr, Ns, inlineVars, 'np', (2,5,10,20))
+        #self._check(x, iv[1,:], zSim[-1,:], eValue)
+
+    def test_optimize_N2(self):
+        # Empirically, it looks like ni=2 is best. It's noticeably better than
+        # ni=1 (though not dramatically so), but going to anything larger
+        # doesn't show any improvement.
+        Nr = 1000
+        wiggle = 2.0
+        eValue = 0.95
+        Ns, N1 = 1000, 10
+        inlineVars = self._inlineVars(wiggle, eValue, N1, 1)
+        self._simRounds(Nr, Ns, inlineVars, 'ni', (1,2,3,4))
+        #self._check(x, iv[1,:], zSim[-1,:], eValue)
-    
-    def test_highCorrelation_c(self):
-
-
-
-
-
+
+
+
+
+
+
-            self._plot((hReal, hSim))
-