Changeset 184

Timestamp:

05/17/08 01:31:43 (8 months ago)

Author:

edsuom

Message:

Grinding along with C code for new model

Files:

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

projects/AsynCluster/trunk/svpmc/model.c

@@ -109 +109 @@
-  // Log-likelihood of the innovations up to the current one, given the history
-  // of h to this point
- = 0
+
-  // Multivariate normal proposal, which is just z for time-series samples
-  // after the first in each log-volatility computation
@@ -127 +127 @@
-  int kp;
-  double Lp;
-
+;
-
-
@@ -271 +271 @@
-// the log-probability of the proposal's having been selected.
-struct selection \
-s
+importanceS
-{
-  int k0, kpMax;
@@ -287 +287 @@
-      pMax = val;
-      kpMax = k0;
+    }
-  }
-
@@ -304 +305 @@
-
-
-hybridGibbs
+likelihood
-// 
-// Supplied variables
@@ -311 +312 @@
-//--- Supplied variables ---
-// z    Independent normal variates, [pn] array (updated)
+// x    Innovation values, [pn] array
-// h    Simulated last-sample multivariate log-volatility, [p] (overwritten)
-// x    Innovation values, [pn] array
-// w    Wiggle value for +/- proposals (float)
-// n    Number of proposals to try per time-series sample (int)
@@ -366 +367 @@
-P.rz = rz_array; P.ri = ri_array; P.sc = sc_array;
-
-arrays
+stuff
-xk = (double *) calloc(sizeof(double), Nt);
-for(k0=0; k0<n; k0++) {
-  for(k1=0; k1<2; k1++) {
+    sp[k0][k1].Lx = 0;
-    sp[k0][k1].h = (double *) calloc(sizeof(double), Nt);
-    sp[k0][k1].x0 = (double *) calloc(sizeof(double), Nt);
@@ -430 +432 @@
-
-  // Importance sample proposal using pProb=exp(Lz+Lx) for each
-s
+importanceS
-  Lp += log(spSelection.Lp);
-  for(k1=0; k1<Nt; k1++)
@@ -453 +455 @@
-
-// Free array memory
-0
+k
-for(k0=0; k0<n; k0++) {
-  for(k1=0; k1<2; k1++) {

projects/AsynCluster/trunk/svpmc/model.py

@@ -56 +56 @@
-        M = modelObj
-
-, sigma
-, sigma
+
+
-        if result is None:
-            return
@@ -106 +106 @@
-        return d
-    
-sigma, 
+
-        """
-        Computes the log-likelihood of the parameters in the supplied
@@ -123 +123 @@
-                # An error from the remote likelihood method call is treated as
-                # infinitely low likelihood
-[], [], 0, 0
+0, [], []
-            elif isinstance(result, str):
-                # Unpack string result from remote call
-                result = list(pack.Unpacker(result))
-h', 'z', 'Lh', 'acceptanceRate
+Lp', 'z', 'h
-                setattr(paramContainer, name, result[k])
-            return paramContainer
-
+
-        if (remote or self.remoteMode) and not local:
-            d = self.client.run(
-sigma, 
+
-        else:
-            d = self.queue.call(
-, sigma
+
-        return d.addCallbacks(gotLikelihood, self._oops)
-
@@ -177 +177 @@
-    @ivar y: A 2-D array containing my observation data.
-
-
-
-
-
+
+
+
+
+
+
+
+
+
+
-
-    #--- Properties -----------------------------------------------------------
@@ -221 +227 @@
-                i += 1
-        return cv
-
-    def logVolatilities(self, z, v, x, h0):
-        """
-        Call this method with a 2-D array I{z} of multivariate IID normal
-        random variates, a 2-D array I{v} of volatility shocks correlated from
-        the IID variates, a 2-D vector I{x} containing one multivariate
-        innovation for each random variate in I{z}, and a 1-D vector I{h0} that
-        defines an initial multivariate log-volatility value.
-
-        The method will compute a 2-D array I{h} of log-volatilities for each
-        volatility shock by running it through the stochastic-volatility VAR(1)
-        process, with the value provided in I{h0} as the VAR component for the
-        first sample. Then it will compute the total log-likelihood of the
-        innovations in I{x} given each respective log-volatility value in I{h}.
-        
-        The method returns a tuple with the total log-likelihood value and the
-        computed log-volatility array I{h}.
-
-        CALLS TO THIS METHOD IS WHERE MOST OF THE CPU TIME IS SPENT!
-        """
-        # Mostly empty log-volatility array
-        h = s.column_stack([h0, s.zeros_like(v)])
-        # Working arrays with intialized values
-        if 'lv_work' in self.cache:
-            y, ri = self.cache['lv_work']
-        else:
-            # Scratchpad with some constants
-            y = s.empty((self.p, 4))
-            y[:,0] = s.sqrt(1 - self.g**2)
-            y[:,1] = s.log(s.sqrt(2*s.pi) * y[:,0])
-            # Inverse of concurrent cross-correlations between innovation
-            # shocks, will raise exception if invalid
-            ri = linalg.inv(linalg.cholesky(
-                self.covarMatrix(s.ones(self.p), self.cr), lower=True))
-            self.cache['lv_work'] = y, ri
-        # Run the C code where the heavy lifting is done
-        L0, L = self.inline('z', 'v', 'x', 'h', 'ri', 'y', 'd', 'e', 'g')
-        return L0, L, h[:,1:]
-
-    def hybridGibbs(self, z, x, rv, sigma, tol=1E-3):
-        """
-        Does one iteration of a hybrid Gibbs sampler for the log-volatilities,
-        where each conditional draw is done with a Metropolis-Hastings step.
-        
-        Updates the IID normal variates in place. Returns the likelihood of the
-        innovations in I{x}, given the simulated log-volatilities, along with
-        the last sample's multivariate log-volatility value, and the
-        log-probability of the simulated result.
-        
-        The method will account for the effect of log-volatility proposals on
-        downstream samples. You can specify the maximum fractional effect to
-        account for with I{tol}. Smaller tolerances require more post-sample
-        log-volatility recomputations and hence more CPU time.
-        """
-        L_total = 0
-        N = self.decaySamples(tol)
-        # Initialize volatility shocks from cross-correlation of the IID
-        # variates, if positive definite...
-        v = s.array(rv*z)
-        # ...and initial log-volatilites, from the offset alone...
-        h0 = self.d
-        # ...and a set of random-walk proposals for the IID variates, along
-        # with their volatility shocks
-        zp, L_result = self.zProposal(z, sigma)
-        vp = s.array(rv*zp)
-        # Do conditional draws of each sample in turn
-        for k in xrange(self.n):
-            # We need to compare the current log-volatilities for
-            # this sample, given the previous log-volatility sample...
-            LV = self.logVolatilities(z[:,k:k+N], v[:,k:k+N], x[:,k:k+N], h0)
-            # ...to a log-volatility proposal
-            zpk = s.column_stack([zp[:,k], z[:,k+1:k+N]])
-            vpk = s.column_stack([vp[:,k], v[:,k+1:k+N]])
-            LVp = self.logVolatilities(zpk, vpk, x[:,k:k+N], h0)
-            # Metropolis-Hastings accept/reject of the proposal
-            dL = LVp[1] - LV[1]
-            if dL < -600 or not s.isfinite(dL):
-                # Bogus likelihood or infintesimally small probability of
-                # acceptance, just reject
-                accept = False
-            elif dL > 0:
-                # Always accept if the proposal's likelihood is better
-                accept = True
-            elif self.logUniform() < dL:
-                # Proposal likelihood worse, but accepted
-                L_result += dL
-                accept = True
-            else:
-                # Proposal likelihood worse and rejected
-                L_result += s.log(1 - s.exp(dL))
-                accept = False
-            if accept:
-                LV = LVp
-                z[:,k] = zp[:,k]
-                v[:,k] = vp[:,k]
-            # Update stuff for the next sample
-            h0 = LV[2][:,0]
-            L_total += LV[0]
-            # If the likelihood at this point is so bad that the parameter is
-            # guaranteed not to be resampled, just bail out now and save time
-            if L_total < -1E6:
-                return -s.inf, s.zeros_like(h0), 0.0, 0.0
-        return L_total, h0, L_result
-    
-
-    #--- The API --------------------------------------------------------------
-    
-, sigma
+
-        """
-        Call this method with a parameter object that defines:
@@ -338 +241 @@
-               draw of log volatilities.
-
-            3. The log-likelihood I{L} of my observations given the parameters
-               and after this method has done a simulation draw of
-               log-volatilities conditional on those parameter values and
-               starting with the IID variates in I{z}.
-
-        A new 2-D array of log-volatilities will be drawn from a hybrid Gibbs
-Metropolis-Hastings
-random walk simulation of the IID variates, followed by
-cross-correlation to form volatility shocks, and VAR(1) to form t
-log-volatilities. The Metropolis-Hastings proposals are scaled based on
-the supplied fractional I{sigma}. The simulation is governed by the
-log-likelihood of my observation data given the supplied parameters an
-the simulated 
+importance-sampling
+simulation of the IID variates, followed by cross-correlation to form
+volatility shocks, and VAR(1) to form the log-volatilities. T
+importance-sampling proposals are scaled based on the supplied
+fractional I{sigma}. The simulation is governed by the log-likelihood
+of my observation data given the supplied parameters and the simulate
+
-
-        If a L{linalg.LinAlgError} is raised due to an invalid correlation
@@ -357 +255 @@
-        the log-volatility simulation:
-        
-value of the log-likelihood.
-
-
-
-
-
-4. The log-probability of the simulated result, normalized to be
-
-
-5. The mean acceptance rate of the Metropolis-Hastings proposals
+log-likelihood of my observations given the simulated
+
+
+
+
+
+   log-volatilities.
+
+
+   sample, useful for extrapolating from the fitted model
-        
-        The returned likelihood does not consider the prior probability of the
@@ -378 +276 @@
-        # ...including the latent parameter of IID normal variates
-        z = paramContainer.z.copy()
-
-
+
+
+
+
+
+
+
+
+
-        try:
-
-
+
+
+
+
-        except linalg.LinAlgError:
-            return
-
-
+
+
+
+
-        # Run the hybrid Gibbs rounds, returning the likelihood from the last
-
-
-
-
-
-
-
-
-
-
-
-
-
+
+
+
+
+
+
+
+
+
+

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

@@ -461 +461 @@
-    struct sample S;
-    struct parameters P;
+    double xk[Nx[0]];
-
-    P.g = g_array; P.ri = ri_array; P.sc = sc_array;
-
-
-
-
-
-
-
+
+
+
+
+
+
+
+
+
+
-        PA1(S.z, k1) = Z2(k1, k0);
+        PA1(S.h, k1) = H2(k1, k0);
+      }
-      // This is what's being tested
-
-
-
+
-    }
+    return_val = S.Lx;
-    """
-    
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
-
-class Test_Model_likelihood(BaseCase):