Changeset 203

Show
Ignore:
Timestamp:
05/30/08 08:27:51 (6 months ago)
Author:
edsuom
Message:

Switching from attempting estimate of L with integration to Rao-Blackwellized evaluation over z sim for all param proposals

Files:

Legend:

Unmodified
Added
Removed
Modified
Copied
Moved

  • projects/AsynCluster/branches/GM_WWT_Normal-svpmc_model.c

    r201 r203  
    172172  static double logC = -log(2 * M_PI); 
    173173  return -0.5*x*x + logC; 
     174} 
     175 
     176 
     177// Normal random variate generation 
     178// 
     179// From http://people.scs.fsu.edu/~burkardt/cpp_src/ziggurat/ziggurat.C and 
     180// G. Marsaglia, et al, "The Ziggurat Method for Generating Random Variables," 
     181// Journal of Statistical Software, Volume 5, Number 8, October 2000. 
     182#define SHR3 (jz=jsr, jsr^=(jsr<<13), jsr^=(jsr>>17), jsr^=(jsr<<5), jz+jsr) 
     183#define UNI (0.5 + (signed) SHR3 * 0.2328306e-09) 
     184float normal(void) 
     185{ 
     186  const float r = 3.442620; 
     187  static float x, y; 
     188  static long int hz; 
     189  static unsigned long int iz, jz; 
     190  static unsigned long int jsr = 123456789; 
     191  static unsigned int ke[256], kn[128]; 
     192  static float fe[256], fn[128], we[256], wn[128]; 
     193 
     194  // Set things up on the first run only 
     195  static int ready = 0; 
     196  if(!ready) { 
     197    const double m1 = 2147483648.0; 
     198    int k; 
     199    double q; 
     200    double dn = 3.442619855899; 
     201    double tn = 3.442619855899; 
     202    double vn = 9.91256303526217e-03; 
     203    q = vn / exp(-0.5 * dn*dn); 
     204    kn[0] = (dn / q) * m1; 
     205    kn[1] = 0; 
     206 
     207    wn[0] = q / m1; 
     208    wn[127] = dn / m1; 
     209 
     210    fn[0] = 1.0; 
     211    fn[127] = exp(-0.5 * dn*dn); 
     212 
     213    for (k=126; 1 <= k; k--) { 
     214      dn = sqrt( -2.0 * log(vn/dn + exp(-0.5 * dn*dn))); 
     215      printf("%f %f %f %f\n", dn, tn, m1, (dn/tn)*m1); 
     216      kn[k+1] = (dn / tn) * m1; 
     217      tn = dn; 
     218      fn[k] = exp(-0.5 * dn*dn); 
     219      wn[k] = dn / m1; 
     220    } 
     221    ready = 1; 
     222  } 
     223 
     224  // Return the quick result if possible... 
     225  hz = SHR3; 
     226  iz = hz & 127; 
     227  x = hz * wn[iz]; 
     228  printf("%f %f\n", fabs(hz), kn[iz]); 
     229  if(fabs(hz) < kn[iz]) { 
     230    printf("A: %f\n",x); 
     231    return x; 
     232  } 
     233 
     234  // ...and compute things otherwise if not 
     235  for (;;) { 
     236    if(iz == 0) { 
     237      do { 
     238        x = -log(UNI) * 0.2904764;  
     239        y = -log(UNI); 
     240      } while (y+y < x*x); 
     241      printf("B: %f\n",(0<hz) ? r+x : -r-x); 
     242      return (0<hz) ? r+x : -r-x; 
     243    } 
     244    if(fn[iz] + UNI*(fn[iz-1] - fn[iz]) < exp(-0.5 * x*x)) { 
     245      printf("C: %f\n",x); 
     246      return x; 
     247    } 
     248    hz = SHR3; 
     249    iz = hz & 127; 
     250    if(fabs(hz) < kn[iz]) { 
     251      printf("D: %f\n",x); 
     252      return hz * wn[iz]; 
     253    } 
     254  } 
    174255} 
    175256 

  • projects/AsynCluster/branches/svpmc_model_Lx-Lz/model.c

    r202 r203  
    343343// The number of multivariate samples 
    344344const int Ns = Nx[1]; 
     345// The last iteration index 
     346const int ki_n = ni-1; 
    345347 
    346348int ki, ke; 
    …… 
    355357// first time-series sample of the evaluation interval 
    356358double he[2][Nt]; 
     359// Accumulator for log-probability of the simulated log-volatilities, relative 
     360// to that of the baseline zero-valued log-volatilities. 
     361double Lz; 
    357362// Log-likelihood of the innovation for the first time-series sample in current 
    358 // and proposal evalutions, working values and accumulated for each iteration 
    359 double Lxk[2], Lx[ni]; 
     363// and proposal evalutions, working values and accumulated for the last iteration 
     364double Lxk[2]; 
     365double Lx = 0; 
    360366 
    361367// A set of evaluation sample structs 
    …… 
    378384} 
    379385 
     386// Start with a baseline, all-zero set of log-volatilities... 
     387for(k0=0; k0<Ns; k0++) { 
     388  for(k1=0; k1<Nt; k1++) 
     389    Z2(k1, k0) = 0; 
     390} 
     391// ...and the negative of their log-densitities, for the reference value of 
     392// Log(P(z|w)) 
     393Lz = -Ns*Nt*normLogDensity(0); 
     394 
    380395// For each iteration... 
    381396for(ki=0; ki<ni; ki++) { 
    382  
    383   // Initialize this iteration's accumulator 
    384   Lx[ki] = 0; 
    385397 
    386398  // The "previous" log-volatility of the first time-series sample is set to 
    …… 
    486498    // 
    487499    // P(x,z|w) = P(x|z,w) * P(z|w) 
    488     // 
    489     // P(z|w) doesn't need to be accounted for, however, because the z's are 
    490     // simulated from it. 
    491  
    492     // Accumulate Log(P(xk,zk|w) = Lxk to eventually get Log(P(x,z|w) 
    493     Lx[ki] += Lxk[ke]; 
     500 
     501    if(ki==ki_n) { 
     502      // Last iteration, accumulate positive Log(P(zk|w)) for log-density... 
     503      for(k2=0; k2<Nt; k2++) 
     504        Lz += normLogDensity(Z2(k2, k0)); 
     505      // ...and Log(P(xk|zk,w) = Lxk for Log(P(x|z,w) 
     506      Lx += Lxk[ke]; 
     507    } 
    494508   
    495509    // Done with log-volatility draw for this time-series sample 
    …… 
    513527} 
    514528 
    515 // Compute the result, the integrated P(x|w) over the last 3/4 of the 
    516 // simulations on z 
    517 k0 = (int)(100*ni/400); 
    518  
    519 // Compute the maximum log density and save it to subtract from everybody so 
    520 // that the maximum log density is normalized to zero 
    521 maxVal = -HUGE_VAL; 
    522 for(ki=k0; ki<ni; ki++) { 
    523   if(Lx[ki] > maxVal) 
    524     maxVal = Lx[ki]; 
    525 } 
    526 // Accumulate the linearized normalized densities 
    527 val = 0; 
    528 for(ki=k0; ki<ni; ki++) 
    529   val += exp(Lx[ki] - maxVal); 
    530 // ...and return the denormalized, log mean of the linear densities... 
    531 return_val = log(val) + maxVal - log(ni-k0); 
     529// Return the log-likelihood of the innovations given the log-volatilities at 
     530// the end of the simulation, Log(P(x|z,w))... 
     531py::tuple result(2); 
     532result[0] = Lx; 
     533// ...along with the log-probability of the simulated log-volatilities, 
     534// Log(P(z|w)) 
     535result[1] = Lz; 
     536return_val = result; 
     537 

  • projects/AsynCluster/branches/svpmc_model_Lx-Lz/model.py

    r202 r203  
    298298 
    299299        # Run the hybrid Gibbs rounds 
    300         Lx = self.inline( 
     300        Lx, Lz = self.inline( 
    301301            'z', 'h', 'x', 'w', 'ni', 
    302302            'd', 'e', 'g', 'rz', 'ri', 'sc', 
    303303            w=self.wiggle, ni=self.Ni) 
    304         return Lx, z, h 
     304        return Lx, Lz, z, h 

  • projects/AsynCluster/branches/svpmc_model_Lx-Lz/test_model.py

    r202 r203  
    637637            if j==0 or not sameZ: 
    638638                inlineVars['z'] = z = s.zeros((1, Ns)) 
    639             Lx = modelObj.inlineDirect( 
     639            Lx, Lz = modelObj.inlineDirect( 
    640640                modelObj.code['likelihood'], **inlineVars) 
    641641            hSim = self._z_to_h(z[0,:], eValue, 1.0) 
    642             yield hReal, hSim, Lx 
     642            yield hReal, hSim, Lx, Lz 
    643643     
    644644    def test_basic(self): 
    …… 
    672672            sp.set_title("%s=%d" % (name, value)) 
    673673 
    674     def test_Lx(self): 
     674    def test_Lx_vs_Lz(self): 
    675675        # Wiggle seems fairly critical, best at 0.5-1.0. 
    676676        # Tried independent proposals +/- 4.0, didn't work well. 
    677677        wiggle = 1.0 
    678678        eValue = 0.97 
    679         Ns, Ni, Nr = 200, 100, 100 
     679        Ns, Ni, Nr = 200, 20, 100 
    680680        simResults = {} 
    681         L = s.empty(Nr) 
     681        Lx = s.empty(Nr) 
     682        Lz = s.empty(Nr) 
    682683        k = 0 
    683684        inlineVars = self._inlineVars(wiggle, eValue, Ni) 
    684         for hReal, hSim, Lx in self._likelihood( 
     685        for hReal, hSim, Lxk, Lzk in self._likelihood( 
    685686            Ns, Nr, inlineVars, sameX=True, sameZ=False): 
    686             L[k] = Lx 
    687             simResults[k] = (hReal, hSim, "Lx=%f" % (Lx,)) 
     687            Lx[k] = Lxk 
     688            Lz[k] = Lzk 
     689            simResults[k] = (hReal, hSim, "Lx=%f, Lz=%f" % (Lxk, Lzk)) 
    688690            if VERBOSE: 
    689                 print "%3d: %4.2f" % (k, Lx) 
     691                print "%3d: %4.2f\t%4.2f" % (k, Lxk, Lzk) 
    690692            k += 1 
    691693        print "Lx: mean=%4.2f, sdev=%4.2f, range=%4.2f" % ( 
    692             L.mean(), L.std(), L.max()-L.min()) 
     694            Lx.mean(), Lx.std(), Lx.max()-Lx.min()) 
     695        print "Lz: mean=%4.2f, sdev=%4.2f, range=%4.2f" % ( 
     696            Lz.mean(), Lz.std(), Lz.max()-Lz.min()) 
    693697        if VERBOSE: 
    694             self._plot(L) 
     698            self._plot(Lx, Lz) 
     699            self.fig.add_subplot(111).plot(Lx, Lz, '.') 
    695700            # Worst 
    696             plotArgs = [simResults[x] for x in L.argsort()[:4]] 
     701            I = Lx.argsort() 
     702            plotArgs = [simResults[x] for x in [I[-1], I[0]]] 
    697703            fig = self._plot(*plotArgs) 
    698             fig.text(0.05, 0.95, "Worst", size='x-large') 
    699             # Best 
    700             plotArgs = [simResults[x] for x in L.argsort()[::-1][:4]] 
    701             fig = self._plot(*plotArgs) 
    702             fig.text(0.05, 0.95, "Best", size='x-large') 
     704            fig.text(0.05, 0.95, "Best & Worst", size='x-large') 
    703705     
    704706    def test_optimize_Ni(self): 

  • projects/AsynCluster/tags/20080530-Integ_L_Estimate/svpmc/params.py

    r196 r203  
    347347        Returns probability density of the supplied parameter value I{x}. 
    348348        """ 
     349        if x.shape: 
     350            return self.distObj.pdf(x) 
    349351        return float(self.distObj.pdf(x)) 
    350352     

  • projects/AsynCluster/tags/20080530-Integ_L_Estimate/svpmc/project.py

    r201 r203  
    9292        self.m = m 
    9393        self.iteration = 0 
    94         #--- Use the NullQueue for debugging, thread queue normally ----------- 
    95         #self.queue = NullQueue() 
    96         self.queue = asynqueue.ThreadQueue(1) 
    97         #---------------------------------------------------------------------- 
    9894        specDir = os.path.dirname(specFile) 
    9995        self.tables = self._parseSpec(specFile) 
    …… 
    114110        # Parameters 
    115111        paramTitles, dimensions = self._setupParams() 
     112        # Open the NetCDF file 
     113        self.cdf = NetCDF.NetCDFFile( 
     114            os.path.join(specDir, ncFile), ('w', 'r')[int(readOnly)]) 
     115        # If in read-only mode, we're done now 
     116        if readOnly: 
     117            return 
     118        # Write mode, write some initial stuff to the NetCDF file... 
     119        self._setupCDF(tsData, seriesTitles, paramTitles, dimensions) 
     120        for k, title in enumerate(seriesTitles): 
     121            obs = self.cdf.variables['observations'] 
     122            obs[k,:] = tsData[k,:].astype('f') 
     123            setattr(obs, "series_%02d" % (k,), title) 
     124        # ...and construct a queue and model manager for a simulation 
     125        self.queue = asynqueue.ThreadQueue(1) 
    116126        self.mgr = model.ModelManager( 
    117127            self, tsData, wiggle=self.wiggle, Ni=self.Ni) 
    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) 
    128128     
    129129    def _parseSpec(self, filePath): 
    …… 
    236236        return titles, dimensions 
    237237 
    238     def _setupCDF(self, filePath, tsData, sTitles, pTitles, pDims): 
    239         """ 
    240         """ 
    241         cdf = NetCDF.NetCDFFile(filePath, 'w') 
    242         self.cdf = cdf 
    243         #self.cdf = CDF_Wrapper(cdf) 
     238    def _setupCDF(self, tsData, sTitles, pTitles, pDims): 
     239        """ 
     240        """ 
     241        cdf = self.cdf 
    244242        cdf.title = self.tables['project'][0][0] 
    245243        # Dimensions 
    …… 
    298296        Call this when the project is done. 
    299297        """ 
    300         d = self.mgr.shutdown() 
    301         d.addCallback(lambda _: self.queue.shutdown()) 
    302         d.addCallback(lambda _: self.cdf.close()) 
    303         return d 
    304  
     298        if hasattr(self, 'mgr'): 
     299            d = self.mgr.shutdown() 
     300            d.addCallback(lambda _: self.queue.shutdown()) 
     301            d.addCallback(lambda _: self.cdf.close()) 
     302            return d 
     303        self.cdf.close() 
     304 

  • projects/AsynCluster/tags/20080530-Integ_L_Estimate/svpmc/test/test_model.py

    r202 r203  
    679679        Ns, Ni, Nr = 200, 100, 100 
    680680        simResults = {} 
    681         L = s.empty(Nr) 
     681        Lx = s.empty(Nr) 
    682682        k = 0 
    683683        inlineVars = self._inlineVars(wiggle, eValue, Ni) 
    684         for hReal, hSim, Lx in self._likelihood( 
     684        for hReal, hSim, Lxk in self._likelihood( 
    685685            Ns, Nr, inlineVars, sameX=True, sameZ=False): 
    686             L[k] = Lx 
    687             simResults[k] = (hReal, hSim, "Lx=%f" % (Lx,)) 
     686            Lx[k] = Lxk 
     687            simResults[k] = (hReal, hSim, "Lx=%f" % (Lxk,)) 
    688688            if VERBOSE: 
    689                 print "%3d: %4.2f" % (k, Lx) 
     689                print "%3d: %4.2f" % (k, Lxk) 
    690690            k += 1 
    691691        print "Lx: mean=%4.2f, sdev=%4.2f, range=%4.2f" % ( 
    692             L.mean(), L.std(), L.max()-L.min()) 
     692            Lx.mean(), Lx.std(), Lx.max()-Lx.min()) 
    693693        if VERBOSE: 
    694             self._plot(L) 
    695             # Worst 
    696             plotArgs = [simResults[x] for x in L.argsort()[:4]] 
     694            self._plot(Lx) 
     695            # Best & Worst 
     696            I = Lx.argsort() 
     697            plotArgs = [simResults[x] for x in [I[-1], I[0]]] 
    697698            fig = self._plot(*plotArgs) 
    698             fig.text(0.05, 0.95, "Worst", size='x-large') 
    699             # Best 
    700             plotArgs = [simResults[x] for x in L.argsort()[::-1][:4]] 
    701             fig = self._plot(*plotArgs) 
    702             fig.text(0.05, 0.95, "Best", size='x-large') 
     699            fig.text(0.05, 0.95, "Best & Worst", size='x-large') 
    703700     
    704701    def test_optimize_Ni(self): 
    …… 
    710707        Nr = 1000 
    711708        wiggle = 1.0 
    712         eValue = 0.95 
     709        eValue = 0.97 
    713710        Ns = 1000 
    714711        inlineVars = self._inlineVars(wiggle, eValue, 1) 

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

    r196 r203  
    347347        Returns probability density of the supplied parameter value I{x}. 
    348348        """ 
     349        if x.shape: 
     350            return self.distObj.pdf(x) 
    349351        return float(self.distObj.pdf(x)) 
    350352     

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

    r201 r203  
    9292        self.m = m 
    9393        self.iteration = 0 
    94         #--- Use the NullQueue for debugging, thread queue normally ----------- 
    95         #self.queue = NullQueue() 
    96         self.queue = asynqueue.ThreadQueue(1) 
    97         #---------------------------------------------------------------------- 
    9894        specDir = os.path.dirname(specFile) 
    9995        self.tables = self._parseSpec(specFile) 
    …… 
    114110        # Parameters 
    115111        paramTitles, dimensions = self._setupParams() 
     112        # Open the NetCDF file 
     113        self.cdf = NetCDF.NetCDFFile( 
     114            os.path.join(specDir, ncFile), ('w', 'r')[int(readOnly)]) 
     115        # If in read-only mode, we're done now 
     116        if readOnly: 
     117            return 
     118        # Write mode, write some initial stuff to the NetCDF file... 
     119        self._setupCDF(tsData, seriesTitles, paramTitles, dimensions) 
     120        for k, title in enumerate(seriesTitles): 
     121            obs = self.cdf.variables['observations'] 
     122            obs[k,:] = tsData[k,:].astype('f') 
     123            setattr(obs, "series_%02d" % (k,), title) 
     124        # ...and construct a queue and model manager for a simulation 
     125        self.queue = asynqueue.ThreadQueue(1) 
    116126        self.mgr = model.ModelManager( 
    117127            self, tsData, wiggle=self.wiggle, Ni=self.Ni) 
    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) 
    128128     
    129129    def _parseSpec(self, filePath): 
    …… 
    236236        return titles, dimensions 
    237237 
    238     def _setupCDF(self, filePath, tsData, sTitles, pTitles, pDims): 
    239         """ 
    240         """ 
    241         cdf = NetCDF.NetCDFFile(filePath, 'w') 
    242         self.cdf = cdf 
    243         #self.cdf = CDF_Wrapper(cdf) 
     238    def _setupCDF(self, tsData, sTitles, pTitles, pDims): 
     239        """ 
     240        """ 
     241        cdf = self.cdf 
    244242        cdf.title = self.tables['project'][0][0] 
    245243        # Dimensions 
    …… 
    298296        Call this when the project is done. 
    299297        """ 
    300         d = self.mgr.shutdown() 
    301         d.addCallback(lambda _: self.queue.shutdown()) 
    302         d.addCallback(lambda _: self.cdf.close()) 
    303         return d 
    304  
     298        if hasattr(self, 'mgr'): 
     299            d = self.mgr.shutdown() 
     300            d.addCallback(lambda _: self.queue.shutdown()) 
     301            d.addCallback(lambda _: self.cdf.close()) 
     302            return d 
     303        self.cdf.close() 
     304 

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

    r202 r203  
    679679        Ns, Ni, Nr = 200, 100, 100 
    680680        simResults = {} 
    681         L = s.empty(Nr) 
     681        Lx = s.empty(Nr) 
    682682        k = 0 
    683683        inlineVars = self._inlineVars(wiggle, eValue, Ni) 
    684         for hReal, hSim, Lx in self._likelihood( 
     684        for hReal, hSim, Lxk in self._likelihood( 
    685685            Ns, Nr, inlineVars, sameX=True, sameZ=False): 
    686             L[k] = Lx 
    687             simResults[k] = (hReal, hSim, "Lx=%f" % (Lx,)) 
     686            Lx[k] = Lxk 
     687            simResults[k] = (hReal, hSim, "Lx=%f" % (Lxk,)) 
    688688            if VERBOSE: 
    689                 print "%3d: %4.2f" % (k, Lx) 
     689                print "%3d: %4.2f" % (k, Lxk) 
    690690            k += 1 
    691691        print "Lx: mean=%4.2f, sdev=%4.2f, range=%4.2f" % ( 
    692             L.mean(), L.std(), L.max()-L.min()) 
     692            Lx.mean(), Lx.std(), Lx.max()-Lx.min()) 
    693693        if VERBOSE: 
    694             self._plot(L) 
    695             # Worst 
    696             plotArgs = [simResults[x] for x in L.argsort()[:4]] 
     694            self._plot(Lx) 
     695            # Best & Worst 
     696            I = Lx.argsort() 
     697            plotArgs = [simResults[x] for x in [I[-1], I[0]]] 
    697698            fig = self._plot(*plotArgs) 
    698             fig.text(0.05, 0.95, "Worst", size='x-large') 
    699             # Best 
    700             plotArgs = [simResults[x] for x in L.argsort()[::-1][:4]] 
    701             fig = self._plot(*plotArgs) 
    702             fig.text(0.05, 0.95, "Best", size='x-large') 
     699            fig.text(0.05, 0.95, "Best & Worst", size='x-large') 
    703700     
    704701    def test_optimize_Ni(self): 
    …… 
    710707        Nr = 1000 
    711708        wiggle = 1.0 
    712         eValue = 0.95 
     709        eValue = 0.97 
    713710        Ns = 1000 
    714711        inlineVars = self._inlineVars(wiggle, eValue, 1)