Changeset 202

Timestamp:

05/28/08 23:43:21 (7 months ago)

Author:

edsuom

Message:

Using normal proposal on z instead of uniform; same result at convergence, but seems to have less Lx variance

Files:

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

projects/AsynCluster/trunk/svpmc/model.c

@@ -175 +175 @@
-
-
+// Random variate from a truncated normal distribution
+double truncnorm(double a, double b)
+{
+  int k;
+  register double x;
+  double c, cp;
+  const double diff = b-a;
+  // Compute the scale value c as the maximum possible pdf (regular normal
+  // distribution) within the truncated range
+  if(a < 0) {
+    if(b > 0) {
+      // Center is within the range, so use its pdf
+      c = 1;
+    } else {
+      // Range is below center, so use pdf at its right edge
+      c = EXP(-0.5*b*b);
+    }
+  } else {
+    // Range is above center, so use pdf at its left edge
+    c = EXP(-0.5*a*a);
+  }
+  // Accept-reject sampling for the truncated normal distribution
+  do {
+    // Generate a proposal within the truncated range
+    x = a + diff*drand48();
+    // Compute the pdf of the proposal, using the regular normal distribution
+    cp = EXP(-0.5*x*x);
+    // Accept (and exit the loop) when the scaled uniform variate is less than
+    // the proposal's pdf
+  } while(c*drand48() > cp);
+  return x;
+}
+
+
-// Matrix multiplication x*y -> z, where the square matrix 'x' is a
-// PyArrayObject and 'y' and 'z' are 1-D arrays.
@@ -234 +268 @@
-    //SPA1(S, x0, k) = exp(-0.5 * SPA1(S, h, k)) * PA1(x, k);
-    //
-
+
+
+
-    SPA1(S, x0, k) = EXP(-0.5 * SPA1(S, h, k)) * PA1(x, k);
-
@@ -310 +346 @@
-int ki, ke;
-int k0, k1, k2, k3;
-
+, maxVal
-
-// Multivariate innovation for one current time-series sample
@@ -319 +355 @@
-// first time-series sample of the evaluation interval
-double he[2][Nt];
-// Log probability density of the value of z for current and proposal
-// evaluations, working values
-double Lzk[2];
-// Log-likelihood of the innovation for the first time-series sample in current
-// and proposal evalutions, working values and accumulated for each iteration
@@ -372 +405 @@
-    // Clear the log-likelihood and normal log-density accumulators of each
-    // evaluation struct for the next pair of evaluations.
- {
+
-      se[k2].Lx = 0;
-      Lzk[k2] = 0;
-      // Lxk[.] will be set to a value in se[.].Lx that has already been
-      // multivariate-accumulated, so it doesn't need clearing here
-    }
-    
-    do {
@@ -395 +426 @@
-              // First time-series sample for the second evaluation, use a
-              // proposal on z instead of z itself
-+= uniform(-w, 
+= truncnorm(val-w, val
-              zp[k3] = val;
-            }
-            // Store the log-density of z for the first time-series sample to
-            // compare current vs. proposal
-            Lzk[k2] += normLogDensity(val);
-          }
-          PA1(se[k2].z, k3) = val;
@@ -430 +458 @@
-    } while (k1<Ns && (val > MIN_DIFF || val < -MIN_DIFF));
-
-
-
-
+
+
-    if(val > 0 || log(uniform(0, 1)) < val) {
-      // Proposal accepted
@@ -486 +513 @@
-}
-
-half
+3/4
-// simulations on z
-
-
+
+
+
-// that the maximum log density is normalized to zero
-v
-ni/2
-v
-v
-
-, borrowing Lzk[0] to do so
-Lzk[0]
-ni/2
-Lzk[0] += exp(Lx[ki] - v
+maxV
+k0
+maxV
+maxV
+
+
+val
+k0
+val += exp(Lx[ki] - maxV
-// ...and return the denormalized, log mean of the linear densities...
-Lzk[0]) + val - log(ni/2
+val) + maxVal - log(ni-k0

projects/AsynCluster/trunk/svpmc/model.py

@@ -187 +187 @@
-    
-    """
-3
+10
-
-    #--- Properties -----------------------------------------------------------

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

@@ -318 +318 @@
-            self.failUnlessEqual(len(s.unique(z)), N)
-
+    def test_support_truncnorm(self):
+        pass
+        code = """
+        int k;
+        for(k=0; k<Nx[0]; k++)
+          X1(k) = truncnorm(a, b);
+        """
+        Ns = 100
+        ranges = ((-1, +1), (-0.5, +1.5), (1.0, 2.0), (3.0, 5.0))
+        fig = self.fig
+        spNumBase = 100*len(ranges) + 11
+        k = 0
+        for a, b in ranges:
+            title = "%+3.1f to %+3.1f" % (a, b)
+            print "\n%s\n" % title
+            x = s.empty(Ns)
+            self.inline(code, x=x, a=a, b=b)
+            sp = fig.add_subplot(spNumBase+k)
+            sp.hist(x, bins=20)
+            sp.set_title(title)
+            k += 1
+        print "\n"
+
-    def test_support_matrixMultiply(self):
-        code = """
@@ -650 +673 @@
-
-    def test_Lx(self):
-
+
+
+
-        eValue = 0.97
-        Ns, Ni, Nr = 200, 100, 100
@@ -676 +701 @@
-            fig = self._plot(*plotArgs)
-            fig.text(0.05, 0.95, "Best", size='x-large')
-        
-    
-    def test_optimize_Ni(self):
@@ -685 +709 @@
-        """
-        Nr = 1000
-0.5
+1.0
-        eValue = 0.95
-        Ns = 1000