Changeset 138

Show
Ignore:
Timestamp:
04/05/08 15:45:20 (9 months ago)
Author:
edsuom
Message:

Working on MCMC for normal random walk

Files:

Legend:

Unmodified
Added
Removed
Modified
Copied
Moved

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

    r137 r138  
    8282    """ 
    8383    def __init__(self, rows, cols): 
    84         self.x = s.rand(rows, cols) 
     84        self.x = s.randn(rows, cols) 
     85        self.lp_x = 0 
    8586     
    8687    def __call__(self, w): 
    87         if s.shape(w) != self.w.shape: 
     88        if s.shape(w) != self.x.shape: 
    8889            raise ValueError( 
    8990                "Wiggle array must have same dimensions as walker array") 

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

    r137 r138  
    2626 
    2727int i, j; 
    28 double xp, pr_x, pr_xp; 
     28double xp, dlp, lp_xp; 
     29double lp_x = 0; 
    2930 
    3031for(i=0; i<Nx[0]; i++) { 
    3132  for(j=0; j<Nx[1]; j++) { 
    3233    xp = X2(i,j) + W2(i,j)*normal(); 
     34    lp_xp = -0.5 * pow(xp, 2); 
     35    dlp = lp_xp - lp_x; 
     36    if(dlp > 0  || log(drand48()) < dlp) { 
     37      X2(i,j) = xp; 
     38      lp_x = lp_xp; 
     39    } 
    3340  } 
    3441} 

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

    r137 r138  
    7070    def test_univariate(self): 
    7171        walker = sample.NormalWalk(1, 1) 
    72         w = s.array([[0.2]]) 
    73         x = s.array([walker(w)[0][0] for k in xrange(50000)]) 
     72        w = s.array([[0.1]]) 
     73        x = s.array([walker(w)[0,0] for k in xrange(100000)]) 
    7474        fig = self.fig 
    7575        sp = fig.add_subplot(211)