Changeset 139

Show
Ignore:
Timestamp:
04/05/08 17:11:44 (9 months ago)
Author:
edsuom
Message:

Got NormalWalk? working efficiently and tested

Files:

Legend:

Unmodified
Added
Removed
Modified
Copied
Moved

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

    r138 r139  
    8181    updated array returned. 
    8282    """ 
     83    m = 30 
     84     
    8385    def __init__(self, rows, cols): 
    8486        self.x = s.randn(rows, cols) 
    85         self.lp_x = 0 
     87        self.p = s.zeros_like(self.x) 
    8688     
    8789    def __call__(self, w): 
    …… 
    8991            raise ValueError( 
    9092                "Wiggle array must have same dimensions as walker array") 
    91         return self.c('x', 'w', w=w) 
     93        return self.c('x', 'p', 'm', 'w', w=w) 
    9294         
    9395             

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

    r138 r139  
    2222 
    2323// Supplied variables 
    24 // 'x', 'w' 
     24// 'x', 'p', 'm', 'w' 
    2525 
    2626 
    27 int i, j; 
    28 double xp, dlp, lp_xp; 
    29 double lp_x = 0; 
     27int i, j, k; 
     28double xp, dp, p_xp; 
    3029 
    3130for(i=0; i<Nx[0]; i++) { 
    3231  for(j=0; j<Nx[1]; j++) { 
    33     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; 
     32    for(k=0; k<m; k++) { 
     33      xp = X2(i,j) + W2(i,j) * (drand48() - 0.5); 
     34      p_xp = -0.5 * pow(xp, 2); 
     35      dp = p_xp - P2(i,j); 
     36      if(dp > 0  || log(drand48()) < dp) { 
     37        X2(i,j) = xp; 
     38        P2(i,j) = p_xp; 
     39      } 
    3940    } 
    4041  } 

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

    r138 r139  
    2525import sample 
    2626import util 
     27 
     28 
     29VERBOSE = False 
    2730 
    2831 
    …… 
    6871 
    6972class Test_NormalWalk(util.TestCase): 
    70     def test_univariate(self): 
     73    def check(self, x): 
     74        p = stats.normaltest(x)[1] 
     75        self.failUnless( 
     76            p > 0.05, "Deviation from normality with significance p=%f" % p) 
     77        if VERBOSE: 
     78            fig = self.fig 
     79            sp = fig.add_subplot(211) 
     80            sp.plot(x) 
     81            sp = fig.add_subplot(212) 
     82            sp.hist(x, bins=20) 
     83 
     84    def test_univariate_bigJumps(self): 
     85        walker = sample.NormalWalk(1, 1) 
     86        w = s.array([[1.0]]) 
     87        x = s.array([walker(w)[0,0] for k in xrange(200)]) 
     88        self.check(x) 
     89 
     90    def test_univariate_smallJumps(self): 
    7191        walker = sample.NormalWalk(1, 1) 
    7292        w = s.array([[0.1]]) 
    73         x = s.array([walker(w)[0,0] for k in xrange(100000)]) 
    74         fig = self.fig 
    75         sp = fig.add_subplot(211) 
    76         sp.plot(x) 
    77         sp = fig.add_subplot(212) 
    78         sp.hist(x, bins=20) 
    79  
     93        x = s.array([walker(w)[0,0] for k in xrange(30000)]) 
     94        self.check(x[::30]) 
    8095             
    8196