Changeset 200

Timestamp:

05/27/08 16:46:38 (6 months ago)

Author:

edsuom

Message:

Integration over log-volatility simulations using Metropolis-within-Gibbs is looking most promising

Files:

• projects/AsynCluster/branches/metropolis-within-gibbs (deleted)
• projects/AsynCluster/trunk/doc/svpmc/example/svpmc.conf (modified) (2 diffs)
• projects/AsynCluster/trunk/svpmc/model.c (modified) (16 diffs)
• projects/AsynCluster/trunk/svpmc/model.py (modified) (4 diffs)
• projects/AsynCluster/trunk/svpmc/pmc.py (modified) (1 diff)
• projects/AsynCluster/trunk/svpmc/project.py (modified) (1 diff)
• projects/AsynCluster/trunk/svpmc/test/test_model.py (modified) (9 diffs)

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

@@ -27 +27 @@
-Variable    Value           Description
--------------------------------------------------------------------------------
-5
+4
-Ds          0.1             Starting jump deviation value
-Df          0.2             Fractional value of each dev from prev (or start)
-wiggle      0.5             Range +/- of uniform random walks on z for h sim
-1          10              Proposals on z for each importance-sampling
-
+i          50              Hybrid-Gibbs rounds for simulation
+
-
-
@@ -68 +68 @@
-a           Distribution    Loc         Scale       a       b
--------------------------------------------------------------------------------
-5       2       5
+2       2       3
-
-

projects/AsynCluster/trunk/svpmc/model.c

@@ -86 +86 @@
-// ----------------------------------------------------------------------------
-
-// Run evaluations until odds of error in selection between proposals null
-// proposal on z will be less than 1/100.
-#define MIN_DIFF 0.01
-
-
-// Model parameters
@@ -108 +104 @@
-
-
- proposal
+n evaluation
-struct sample
-{
@@ -117 +113 @@
-  // of h to this point
-  double Lx;
-
-
+
+
+
-  double *z;
-proposal
+evaluation
-  // previous proposal to use this struct if the current proposal hasn't been
-  // crunched yet
@@ -126 +123 @@
-  // Scratchpad 1-D arrays of the same length as z, h
-  double *x0, *x1;
-};
-
-
-// A selection of a proposed log-volatility simulation
-struct selection
-{
-  int kp;
-  double Lp;
-};
-
@@ -219 +208 @@
-  // Inverse-scale the innovation by the square root of the volatilities in
-  // preparation for decorrelation
+
-  for(k=0; k<N; k++)
-
+
+
+
+
-    SPA1(S, x0, k) = exp(-0.5 * SPA1(S, h, k)) * PA1(x, k);
-
+
-  // Now decorrelate the equi-variance innovations in preparation for the
-  // log-likelihood computation for this multivariate log-volatility proposal
@@ -232 +225 @@
-
-  S->Lxk = 0;
-
+
-  // For each time series...
-  for(k=0; k<N; k++) {
@@ -252 +245 @@
-    S->Lxk -= 0.5 * SPA1(S, h, k);
-  }
+  // Add the current log-likelihood to the accumulated log-likelihoods
-  S->Lx += S->Lxk;
-}
-
-
-// Returns with zero only if the non-null proposals have converged close enough
-// to the baseline value (the null proposal) to finish likelihood computation
-int notCloseEnough(struct sample sp[], int n)
-{
-  int k;
-  double L;
-  for(k=1; k<n; k++) {
-    L = sp[k].Lxk - sp[0].Lxk;
-    if(L > MIN_DIFF || L < -MIN_DIFF)
-      return 1;
-  }
-  return 0;
-}
-
-
-// Importance-samples one of the non-null proposals based on a weight computed
-// from the log-density of z and the log-likelihood of the innovation given the
-// proposal's log-volatility, conditional on the null proposal.
-//
-// Returns a struct with an integer index to the sampled proposal's struct and
-// the log-probability of the proposal's having been selected.
-struct selection \
-importanceSample(struct sample sp[], double Lz[], int n)
-{
-  int k;
-  double val;
-  double Dp_sum = 0;
-  double Lp[n], Dp[n];
-  double L_max = -HUGE_VAL; // Anything will be more positive
-  struct selection result;
-
-  // Log-densities densities are relative to the null-proposal, and with a
-  // maximum for normalization in selection
-  for(k=1; k<n; k++) {
-    val = sp[k].Lx + Lz[k] - sp[0].Lx - Lz[0];
-    if(val > L_max)
-      L_max = val;
-    Lp[k] = val;
-  }
-
-  // Compute linear probability density for each proposal, relative to the
-  // maximum one
-  for(k=1; k<n; k++) {
-    val = exp(Lp[k] - L_max);
-    Dp[k] = val;
-    Dp_sum += val;
-  }
-  
-  // Accept-reject sampling with uniform random proposal selection
-  do {
-    // Randomly select one of the proposals...
-    result.kp = (int) uniform(1, n);
-    // ...until the selected proposal is accepted, with probability
-    // proportional to its relative weight
-  } while(uniform(0, 1.0) > Dp[result.kp]);
-
-  // The probability density of the selected log-volatility proposal
-  result.Lp = log(Dp[result.kp] / Dp_sum);
-  return result;
-}
-
@@ -330 +262 @@
-// w    Wiggle value for +/- proposals (float)
-// ni   Number of hybrid-Gibbs iterations (int)
-// np   Number of proposals to try per log-volatility simulation (int)
-
-//--- Model parameters, direct ---
@@ -342 +273 @@
-// sc   Constants for multivariate normal PDF
-
-
-
-
-
+
+
+
+
+
+
+
+
-
-// The number of time series making up a single multivariate sample
@@ -352 +287 @@
-const int Ns = Nx[1];
-
-
+
+
-double val;
-double Lx = 0;
-double Lh = 0;
-
-// Multivariate innovation for one current time-series sample
-double *xk;
-
-
-
-
-
-
-
-
-
-
-
-
-
+
+
+
+
+
+
+
+
+
+
+
+
+
+
-// A struct for the model parameters
-struct parameters P;
-
-// The scalar results go here
-py::tuple result(2);
-
-// Generate parameter struct
@@ -384 +316 @@
-// Initialize various arrays
-xk = (double *) malloc(sizeof(double) * Nt);
-
-
-
-
-
-
-
+
+
+
+
+
+
-
-// For each iteration...
-for(ki=0; ki<ni; ki++) {
+
+  // Initialize this iteration's accumulator
+  Lx[ki] = 0;
-
-  // The "previous" log-volatility of the first time-series sample is set to
-  // the log-volatility offset plus a simulated, latent-parameter value.
-  // TODO
-np
+2
-    for(k1=0; k1<Nt; k1++)
-p
+e
-  }
-
@@ -406 +340 @@
-  for(k0=0; k0<Ns; k0++) {
-    
-
-
-
-
-
+
+
+
+
+
+
-    k1 = k0;
-  
-
-
-
-
-
+
+
+
+
+
+
+
+
+
-    do {
-  
@@ -424 +363 @@
-        PA1(xk, k2) = X2(k2, k1);
-      
-proposal
-np
+evaluation
+2
-  
-        // Set z...
-        for(k3=0; k3<Nt; k3++) {
-          val = Z2(k3, k1);
-
-
-
-
-
+
+
+
+
-              val += uniform(-w, +w);
-
-
+
+
+
+
+
-          }
-p
+e
-        }
-        // ...compute the multivariate log-volatility
-p
+e
-        // ...and the log-likelihood of the innovation for this time-series
-        // sample, given the log-volatility
-p
+e
-        
-        // Save the first log-volatility and the log-likelihood based on it (and
@@ -452 +393 @@
-        if(k1==k0) {
-          for(k3=0; k3<Nt; k3++)
-p[k2][k3] = PA1(sp
-p
+e[k2][k3] = PA1(se
+e
-        }
-      }
@@ -459 +400 @@
-      // Ready for the next time series sample, assuming there is one
-      k1++;
+
+      // Continue only if difference in log-likelihood of the current
+      // time-series sample of the evaluations between current and proposed z
+      // is great enough to warrant doing so.
+      val = se[0].Lxk - se[1].Lxk;
-    
-notCloseEnough(sp, np
-
-Importance sample proposal using pProb=exp(Lz+Lx) for each...
-spSelection = importanceSample(sp, Lz, np);
-// ...add the log-likelihood of the innovation given the selected
-// log-volatility proposal to what will be the total log-likelihood of the
-// innovations given the entire simulated h
-Lx += Lxk[spSelection.kp]
-// ...add the log-density of the selected log-volatility proposal to what
-// will be the total log-density of the entire simulated h...
-Lh += spSelection.Lp
-// ...and update the normal variate underlying the current time-series sample
-// to that on which the selected log-volatility proposal was based.
-for(k1=0; k1<Nt; k1++)
-  Z2(k1, k0) = zp[spSelection.kp];
-
-proposals for the next time-series
-sample to the log-volatility of the selected proposal for this time-series
-
+(val > MIN_DIFF || val < -MIN_DIFF
+
+Metropolis-hastings selection of current or proposal using Lz+Lx for
+// each evaluation
+val = Lzk[1] + se[1].Lx  -  Lzk[0] - se[0].Lx;
+if(val > 0 || log(uniform(0, 1)) < val) {
+  // Proposal accepted
+  ke = 1
+  // ...overwrite z for this sample with the proposal on z
+  for(k1=0; k1<Nt; k1++)
+    Z2(k1, k0) = zp[k1]
+} else {
+  // Proposal rejected
+  ke = 0;
+}
+
+evaluations for the next
+time-series sample to the log-volatility of the selected evaluation for
+this time-series 
-    for(k1=0; k1<Nt; k1++) {
-p[spSelection.kp
-np
-p
+e[ke
+2
+e
-    }
+
+    // The joint probability of the innovation and the variates underlying the
+    // log-volatility simulation, conditional on the parameters, is equal to
+    // the likelihood of the innovation, conditional on the simulated
+    // variates and the parameters, times the probability of the
+    // simulated variates, conditional on the parameters:
+    //
+    // P(x,z|w) = P(x|z,w) * P(z|w)
+    //
+    // P(z|w) is conditional on w because of the correlation between sequential
+    // z values due to the VAR(1) in volatility shocks. Fortunately, the way
+    // that we've simulated z allows us to just evaluate P(z) using the joint
+    // normal distribution.
+
+    // Accumulate Log(P(xk,zk|w) = Lxk + Lzk) to eventually get Log(P(x,z|w)
+    Lx[ki] += Lxk[ke] + Lzk[ke];
-  
-    // Done with log-volatility draw for this time-series sample
@@ -494 +456 @@
-// sample in the last iteration
-for(k0=0; k0<Nt; k0++)
-p[spSelection.kp
+e[ke
-
-// Free array memory
-free(xk);
-
-
-
-
-
-
-
-
-
-
-
-
-
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+

projects/AsynCluster/trunk/svpmc/model.py

@@ -129 +129 @@
-                # Unpack string result from remote call
-                result = list(pack.Unpacker(result))
-Lh', '
+
-                setattr(paramContainer, name, result[k])
-            return paramContainer
@@ -184 +184 @@
-      generator.
-
-
-
-
-
-
-
-
-
-
-
-
+
+
+
+
-
-    #--- Properties -----------------------------------------------------------
@@ -260 +253 @@
-        If a L{linalg.LinAlgError} is raised due to an invalid correlation
-        matrix parameter, the method returns with no result. Otherwise, it
-
-
+
-        
-simulated
-log-volatility values.
-
-2. The log-probability of the simulated log-volatilities.
-
-3. The IID normal variates underlying the simulated
-
-
-4. The multivariate log-volatility value for the last time-series
-sample, useful for extrapolating from the fitted model.
+parameters, from
+an integration of the join probability of the observations and
+
+   
+
+   log-volatilities after the last simulation round.
+
+
+   sample after the last simulation round. (Useful for
+extrapolating from the fitted model.)
-        
-        The returned likelihood does not consider the prior probability of the
@@ -305 +297 @@
-        sc[:,1] = s.log(s.sqrt(2*s.pi) * sc[:,0])
-
-
-
-
-
-
+
+
+
-            'd', 'e', 'g', 'rz', 'ri', 'sc',
-p=self.N1+1, ni=self.N2
-Lh, 
+i=self.Ni
+

projects/AsynCluster/trunk/svpmc/pmc.py

@@ -198 +198 @@
-        def weight(paramContainer):
-            if s.isfinite(paramContainer.Lx):
-
-
+
-                info = " ".join([
-                    "%+12.2f" % (getattr(paramContainer, x),)
-, 'Lh'
+
-            else:
-                info = ""

projects/AsynCluster/trunk/svpmc/project.py

@@ -112 +112 @@
-        paramTitles, dimensions = self._setupParams()
-        self.mgr = model.ModelManager(
-1=self.N1, N2=self.N2
+i=self.Ni
-        self._setupCDF(
-            os.path.join(specDir, ncFileName),

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

@@ -506 +506 @@
-
-
-class Test_Model_importanceSample(BaseCase):
-    code = """
-    int k;
-    const int n = NX[0];
-    struct sample sp[n];
-    double Lz[n];
-    py::tuple result(2);
-    struct selection spSelection;
-
-    for(k=0; k<n; k++) {
-      sp[k].Lx = X1(k);
-      Lz[k] = Z1(k);
-    }
-    spSelection = importanceSample(sp, Lz, n);
-    result[0] = spSelection.kp;
-    result[1] = spSelection.Lp;
-    return_val = result;
-    """
-
-    def _importanceSample(self, Lx_array, Lz_array=None):
-        if Lz_array is None:
-            Lz_array = s.zeros_like(Lx_array)
-        return self.inline(self.code, X=Lx_array, Z=Lz_array)
-
-    def test_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 _check(self, z0, wiggle):
-        N = 20
-        N_iter = 1000
-        distObj = stats.distributions.truncnorm(z0-wiggle, z0+wiggle)
-        X = 2*wiggle*(s.rand(N_iter, N) - 0.5) + z0
-        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], s.pi/N*distObj.pdf(Y[I]))
-        self.failUnlessNormal(Y)
-
-    def test_full(self):
-        return self._check(0, 4.0)
-
-    def test_closePortion(self):
-        return self._check(1.5, 1.0)
-
-    def test_farPortion(self):
-        return self._check(2.5, 0.5)
-
-
-class Test_Model_likelihood(BaseCase):
-    corr = 0.0
@@ -602 +542 @@
-        return iv, h, x
-
-1, N2=1
+i
-        return {
-            'h':s.empty(1),
-            'sc':s.array([1, s.log(s.sqrt(2*s.pi))]),
-            'w':wiggle,
-
-
+
-            'd':s.array([0.0]),
-            'e':s.array([[float(eValue)]]),
@@ -616 +555 @@
-            }
-    
-1, N2
+i
-        """
-        Returns a parameter container populated from my I{params} dict,
@@ -632 +571 @@
-                s.array(kw.get(name, self.params[name]), dtype='d'))
-        # Return the parameter container along with a model to be tested
-1=N1, N2=N2
+i=Ni
-        modelObj.paramNames = util.PARAM_NAMES
-        return pc, modelObj
@@ -650 +589 @@
-        """
-        Runs C code for L{model.Model.likelihood} independent of the method
-
+, repeating the run I{Nr} times with different simulated data
-        """
-        z = s.zeros((1, Ns))
@@ -666 +605 @@
-        wiggle = 0.1
-        eValue = 0.90
-1, N2 = 10, 5, 2
-N1, N2
- 
+r = 6, 5
+1
+
-        modelObj = model.Model()
-Lh     
+
-        print "-"*100
-2
-, Lh
+r
+
-                modelObj.code['likelihood'], **inlineVars)
-            zLine = "   ".join(["%+4.2f" % xx for xx in inlineVars['z'][0,:]])
-+7.1f%s%s" % (k, Lx, Lh
-
+s%s" % (k, Lx
+
-    def _simRounds(self, Nr, Ns, inlineVars, name, values):
-        fig = self.fig
-        bins = s.linspace(0.65, 0.95, 20)
-        for k, value in enumerate(values):
-            print k, value
-            inlineVars[name] = value
-            correlations = []
@@ -694 +632 @@
-            sp.set_title("%s=%d" % (name, value))
-
-
-
-
-
-
+
+
+
+
+
+
+
+
-        wiggle = 0.5
-        eValue = 0.95
-
-
-
-
-
-
-
-
-
-
-
-
-
+
+
-        self._simRounds(Nr, Ns, inlineVars, 'ni', (5,10,20,40))
-    
@@ -719 +649 @@
-        wiggle = 0.5
-        eValue = 0.97
-1, N2 = 200, 1
-1, N2
+i = 20
+i
-        for hReal, hSim in self._likelihood(Ns, 1, inlineVars):
-            self._plot((hReal, hSim))
@@ -727 +657 @@
-        wiggle = 0.5
-        eValue = 0.97
-1, N2 = 1000, 1
+i = 100
-        iv, h, x = self._simData(Ns, eValue)
-1, N2
+i
-        pc.z = s.zeros((1, Ns))
-        modelObj.y = self._x_to_y(x)
-Lh, 
+
-        self._plot(x, (iv[1,:], z[0,:]), (h, self._z_to_h(z[0,:], eValue)))