Changeset 201

Show
Ignore:
Timestamp:
05/28/08 17:18:11 (6 months ago)
Author:
edsuom
Message:

Working on reducing variance in likelihood estimates

Files:

Legend:

Unmodified
Added
Removed
Modified
Copied
Moved

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

    r200 r201  
    2727Variable    Value           Description 
    2828------------------------------------------------------------------------------- 
    29 D           4               Number of jump deviations in the D-kernel sim 
     29D           5               Number of jump deviations in the D-kernel sim 
    3030Ds          0.1             Starting jump deviation value 
    31 Df          0.2             Fractional value of each dev from prev (or start) 
    32 wiggle      0.5             Range +/- of uniform random walks on z for h sim 
     31Df          0.5             Fractional value of each dev from prev (or start) 
     32wiggle      0.3             Range +/- of uniform random walks on z for h sim 
    3333Ni          50              Hybrid-Gibbs rounds for simulation of h 
    3434 

  • projects/AsynCluster/trunk/svpmc/model.c

    r200 r201  
    126126 
    127127 
     128// Fast approximate exponent function EXP 
     129// 
     130// From N. Schraudolph, "A Fast, Compact Approximation of the Exponential 
     131// Function," Neural Computation 11, 853-862 (1999) 
     132static union { 
     133  double d; 
     134  struct { 
     135#ifdef LITTLE_ENDIAN 
     136    int j, i; 
     137#else 
     138    int i, j; 
     139#endif 
     140  } n; 
     141} _eco; 
     142#define EXP_A (1048576/M_LN2) 
     143#define EXP_C 60801 
     144 
     145#define EXP(y) (_eco.n.i = EXP_A*(y) + (1072693248 - EXP_C), _eco.d) 
     146  
     147 
    128148// Uniformly distributed random variate in [a,b) 
    129149double uniform(double a, double b) 
    …… 
    132152  static int seeded = 0; 
    133153  if (!seeded) { 
    134     int i; 
    135     int j = 0; 
     154    int k; 
    136155    unsigned short seed[3]; 
    137     for(i=0;i<3;i++) { 
    138       do { 
    139         j++; 
    140       } while (j < (i+1)*10000); 
    141       seed[i] = clock()*clock() % 65537; 
     156    long int timeval = time(0); 
     157    for(k=0; k<3; k++) { 
     158      seed[k] = (unsigned short) timeval; 
     159      timeval = timeval >> 16; 
    142160    } 
    143161    seed48(seed); 
    …… 
    202220{ 
    203221  register int k; 
    204   double *y; 
    205222  double val; 
    206223  const int N = (int) P->g->dimensions[0]; 
    …… 
    214231    // Intel QX9760 running at 3.8 GHz, though some of that is the unavoidable 
    215232    // array lookups and products in the equation. 
    216     SPA1(S, x0, k) = exp(-0.5 * SPA1(S, h, k)) * PA1(x, k); 
     233    // 
     234    //SPA1(S, x0, k) = exp(-0.5 * SPA1(S, h, k)) * PA1(x, k); 
     235    // 
     236    // This is very fast, but has error of about +/-2% 
     237    SPA1(S, x0, k) = EXP(-0.5 * SPA1(S, h, k)) * PA1(x, k); 
    217238 
    218239  // Now decorrelate the equi-variance innovations in preparation for the 
    …… 
    237258    // exponential being done in log space by subtraction of the log of the 
    238259    // scaling coefficient 
    239     val = SPA1(S, x0, k) / SPP2(P, sc, k, 0); 
     260    val = SPP2(P, sc, k, 0) * SPA1(S, x0, k); 
    240261    // Fast squaring instead of pow() 
    241262    S->Lxk -= 0.5 * val*val + SPP2(P, sc, k, 1); 
    …… 
    309330// A struct for the model parameters 
    310331struct parameters P; 
     332 
    311333 
    312334// Generate parameter struct 
    …… 
    412434    val = Lzk[1] + se[1].Lx  -  Lzk[0] - se[0].Lx; 
    413435    if(val > 0 || log(uniform(0, 1)) < val) { 
    414       // Proposal accepted... 
     436      // Proposal accepted 
    415437      ke = 1; 
    416       // ...overwrite z for this sample with the proposal on z 
    417       for(k1=0; k1<Nt; k1++) 
    418         Z2(k1, k0) = zp[k1]; 
     438      for(k2=0; k2<Nt; k2++) 
     439        Z2(k2, k0) = zp[k2]; 
    419440    } else { 
    420441      // Proposal rejected 
    421442      ke = 0; 
    422443    } 
    423  
     444     
    424445    // Set the initial log-volatility of the evaluations for the next 
    425446    // time-series sample to the log-volatility of the selected evaluation for 
    426447    // this time-series sample. 
    427     for(k1=0; k1<Nt; k1++) { 
    428       val = he[ke][k1]; 
    429       for(k2=0; k2<2; k2++) 
    430         PA1(se[k2].h, k1) = val; 
     448    for(k2=0; k2<Nt; k2++) { 
     449      val = he[ke][k2]; 
     450      for(k3=0; k3<2; k3++) 
     451        PA1(se[k3].h, k2) = val; 
    431452    } 
    432453 
    …… 
    439460    // P(x,z|w) = P(x|z,w) * P(z|w) 
    440461    // 
    441     // P(z|w) is conditional on w because of the correlation between sequential 
    442     // z values due to the VAR(1) in volatility shocks. Fortunately, the way 
    443     // that we've simulated z allows us to just evaluate P(z) using the joint 
    444     // normal distribution. 
    445  
    446     // Accumulate Log(P(xk,zk|w) = Lxk + Lzk) to eventually get Log(P(x,z|w) 
    447     Lx[ki] += Lxk[ke] + Lzk[ke]; 
     462    // P(z|w) doesn't need to be accounted for, however, because the z's are 
     463    // simulated from it. 
     464 
     465    // Accumulate Log(P(xk,zk|w) = Lxk to eventually get Log(P(x,z|w) 
     466    Lx[ki] += Lxk[ke]; 
    448467   
    449468    // Done with log-volatility draw for this time-series sample 
    …… 
    467486} 
    468487 
    469 // Compute the result, the integrated P(x|w) over the simulations on z 
     488// Compute the result, the integrated P(x|w) over the last half of the 
     489// simulations on z 
    470490 
    471491// Compute the maximum log density // and save it to subtract from everybody 
    472492// that the maximum log density is normalized to zero 
    473493val = -HUGE_VAL; 
    474 for(ki=0; ki<ni; ki++) { 
     494for(ki=ni/2; ki<ni; ki++) { 
    475495  if(Lx[ki] > val) 
    476496    val = Lx[ki]; 
    …… 
    478498// Accumulate the linearized normalized densities, borrowing Lzk[0] to do so 
    479499Lzk[0] = 0; 
    480 for(ki=0; ki<ni; ki++) 
     500for(ki=ni/2; ki<ni; ki++) 
    481501  Lzk[0] += exp(Lx[ki] - val); 
    482 // ...and return the denormalized, log mean of the linear densities 
    483 return_val = log(Lzk[0]) + val - log(ni); 
     502// ...and return the denormalized, log mean of the linear densities... 
     503return_val = log(Lzk[0]) + val - log(ni/2); 

  • projects/AsynCluster/trunk/svpmc/model.py

    r200 r201  
    258258               an integration of the join probability of the observations and 
    259259               simulated log-volatilities over the simulation rounds. 
    260                 
     260 
    261261            2. The IID normal variates underlying the simulated 
    262262               log-volatilities after the last simulation round. 
    …… 
    294294        # ...scaling constants for multivariate normal pdf 
    295295        sc = s.empty((self.p, 2)) 
    296         sc[:,0] = s.sqrt(1.0 - self.g**2) 
    297         sc[:,1] = s.log(s.sqrt(2*s.pi) * sc[:,0]) 
     296        sc[:,0] = 1.0 / s.sqrt(1.0 - self.g**2) 
     297        sc[:,1] = s.log(s.sqrt(2*s.pi) / sc[:,0]) 
    298298 
    299299        # Run the hybrid Gibbs rounds 

  • projects/AsynCluster/trunk/svpmc/pmc.py

    r200 r201  
    140140 
    141141    """ 
    142     chunkSize = 100 
     142    chunkSize = 200 
    143143     
    144144    def __init__(self, projectManager, socket=None): 

  • projects/AsynCluster/trunk/svpmc/project.py

    r200 r201  
    5858    @ivar p: The number of time series. 
    5959 
     60    @ivar tables: A dict whose entries contain rows of tables parsed from my 
     61      project specification. 
     62 
    6063    @ivar xcorrs: The number of cross-correlations between series, computed as 
    6164      C{(p**2 + p)/2} 
    …… 
    8689        } 
    8790     
    88     def __init__(self, specFile, ncFileName="svpmc.nc", m=10000): 
     91    def __init__(self, specFile, ncFile="svpmc.nc", m=1000, readOnly=False): 
    8992        self.m = m 
    9093        self.iteration = 0 
    …… 
    113116        self.mgr = model.ModelManager( 
    114117            self, tsData, wiggle=self.wiggle, Ni=self.Ni) 
    115         self._setupCDF( 
    116             os.path.join(specDir, ncFileName), 
    117             tsData, seriesTitles, paramTitles, dimensions) 
    118         # Write the observation data 
    119         for k, title in enumerate(seriesTitles): 
    120             obs = self.cdf.variables['observations'] 
    121             obs[k,:] = tsData[k,:].astype('f') 
    122             setattr(obs, "series_%02d" % (k,), title) 
    123  
     118        # If not in read-only mode, setup the NetCDF file for writing 
     119        if not readOnly: 
     120            self._setupCDF( 
     121                os.path.join(specDir, ncFile), 
     122                tsData, seriesTitles, paramTitles, dimensions) 
     123            # Write the observation data 
     124            for k, title in enumerate(seriesTitles): 
     125                obs = self.cdf.variables['observations'] 
     126                obs[k,:] = tsData[k,:].astype('f') 
     127                setattr(obs, "series_%02d" % (k,), title) 
     128     
    124129    def _parseSpec(self, filePath): 
    125130        tables = {} 
    …… 
    217222        self.paramNames = [] 
    218223        titles, dimensions = [], [] 
     224        self.paramInfo = {} 
    219225        xcorrs = self.xcorrs = int(0.5*(self.p**2 + self.p)) 
    220226        self.priors = params.PriorContainer(self.p, self.n) 
    221227        for name, dims, title in self.tables['parameter']: 
    222228            shape = parseDims(dims) 
     229            self.paramInfo[name] = shape, title 
    223230            pa = self.priors.priorArray(name, *shape) 
    224231            for priorSpec in self.tables[name]: 
    …… 
    259266        # Done setting up, save what we've got thus far 
    260267        cdf.sync() 
    261              
     268 
    262269    def writeParams(self, parameters): 
    263270        """ 

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

    r200 r201  
    337337        self.failUnlessElementsEqual(r, s.matrix(a)*b) 
    338338 
     339    def test_support_fastExp(self): 
     340        code = """ 
     341        int k; 
     342        for(k=0; k<Nx[0]; k++) 
     343          Y1(k) = EXP(X1(k)); 
     344        """ 
     345        x = s.linspace(-500, +500, 3000) 
     346        y = s.empty_like(x) 
     347        self.inline(code, x=x, y=y) 
     348        self.fig.add_subplot(111).plot(x, y/s.exp(x)) 
     349         
    339350 
    340351class Test_Model_logVol(BaseCase): 
    …… 
    583594                vector = [vector] 
    584595            for v in vector: 
    585                 sp.plot(v) 
     596                if isinstance(v, str): 
     597                    sp.set_title(v) 
     598                else: 
     599                    sp.plot(v) 
    586600            sp.grid(True) 
    587      
    588     def _likelihood(self, Ns, Nr, inlineVars): 
     601        return fig 
     602     
     603    def _likelihood(self, Ns, Nr, inlineVars, sameX=False, sameZ=True): 
    589604        """ 
    590605        Runs C code for L{model.Model.likelihood} independent of the method 
    591606        itself, repeating the run I{Nr} times with different simulated data. 
    592607        """ 
    593         z = s.zeros((1, Ns)) 
    594         inlineVars['z'] = z 
    595608        modelObj = model.Model() 
    596609        eValue = inlineVars['e'][0][0] 
    597610        for j in xrange(Nr): 
    598             iv, hReal, x = self._simData(Ns, eValue) 
    599             inlineVars['x'] = s.row_stack([x]) 
    600             modelObj.inlineDirect(modelObj.code['likelihood'], **inlineVars) 
     611            if j==0 or not sameX: 
     612                iv, hReal, x = self._simData(Ns, eValue) 
     613                inlineVars['x'] = s.row_stack([x]) 
     614            if j==0 or not sameZ: 
     615                inlineVars['z'] = z = s.zeros((1, Ns)) 
     616            Lx = modelObj.inlineDirect( 
     617                modelObj.code['likelihood'], **inlineVars) 
    601618            hSim = self._z_to_h(z[0,:], eValue, 1.0) 
    602             yield hReal, hSim 
     619            yield hReal, hSim, Lx 
    603620     
    604621    def test_basic(self): 
    605622        wiggle = 0.1 
    606         eValue = 0.90 
     623        eValue = 0.80 
    607624        Ns, Nr = 6, 50 
    608625        inlineVars = self._inlineVars(wiggle, eValue, 1) 
    609626        inlineVars.update({'x':2*s.ones((1,Ns)), 'z':s.zeros((1, Ns))}) 
    610627        modelObj = model.Model() 
    611         print "\n  k      Lx       Simulated Volatility Shocks" 
     628        print "\n  k      Lx     Simulated Volatility Shocks" 
    612629        print "-"*100 
    613630        for k in xrange(Nr): 
    …… 
    615632                modelObj.code['likelihood'], **inlineVars) 
    616633            zLine = "   ".join(["%+4.2f" % xx for xx in inlineVars['z'][0,:]]) 
    617             print "%3d: %+6.1f %s%s" % (k, Lx, " "*5, zLine) 
     634            print "%3d: %+6.1f %s%s" % (k, Lx, " "*5, zLine) 
    618635     
    619636    def _simRounds(self, Nr, Ns, inlineVars, name, values): 
    …… 
    623640            inlineVars[name] = value 
    624641            correlations = [] 
    625             for hReal, hSim in self._likelihood(Ns, Nr, inlineVars): 
     642            for hReal, hSim, Lx in self._likelihood(Ns, Nr, inlineVars): 
    626643                corrInfo = stats.spearmanr(hReal, hSim) 
    627644                if VERBOSE: 
    …… 
    632649            sp.set_title("%s=%d" % (name, value)) 
    633650 
    634  
     651    def test_Lx(self): 
     652        wiggle = 0.5 
     653        eValue = 0.97 
     654        Ns, Ni, Nr = 200, 100, 100 
     655        simResults = {} 
     656        L = s.empty(Nr) 
     657        k = 0 
     658        inlineVars = self._inlineVars(wiggle, eValue, Ni) 
     659        for hReal, hSim, Lx in self._likelihood( 
     660            Ns, Nr, inlineVars, sameX=True, sameZ=False): 
     661            L[k] = Lx 
     662            simResults[k] = (hReal, hSim, "Lx=%f" % (Lx,)) 
     663            if VERBOSE: 
     664                print "%3d: %4.2f" % (k, Lx) 
     665            k += 1 
     666        print "Lx: mean=%4.2f, sdev=%4.2f, range=%4.2f" % ( 
     667            L.mean(), L.std(), L.max()-L.min()) 
     668        if VERBOSE: 
     669            self._plot(L) 
     670            # Worst 
     671            plotArgs = [simResults[x] for x in L.argsort()[:4]] 
     672            fig = self._plot(*plotArgs) 
     673            fig.text(0.05, 0.95, "Worst", size='x-large') 
     674            # Best 
     675            plotArgs = [simResults[x] for x in L.argsort()[::-1][:4]] 
     676            fig = self._plot(*plotArgs) 
     677            fig.text(0.05, 0.95, "Best", size='x-large') 
     678         
     679     
    635680    def test_optimize_Ni(self): 
    636681        """ 
    …… 
    651696        Ns, Ni = 200, 20 
    652697        inlineVars = self._inlineVars(wiggle, eValue, Ni) 
    653         for hReal, hSim in self._likelihood(Ns, 1, inlineVars): 
    654             self._plot((hReal, hSim)) 
     698        for hReal, hSim, Lx, in self._likelihood(Ns, 1, inlineVars): 
     699            self._plot((hReal, hSim, "Lx=%f" % (Lx,))) 
    655700     
    656701    def test_highCorrelation_model(self):