hnuts

commit 25c15b9f2eaaa13dc844d6ad85a25a17944cdfa3
parent 5aebd6f773fa2bbbff1378fd634c43d82a8ce0bd
Author: Jared Tobin <jared@jtobin.ca>
Date:   Tue, 15 Oct 2013 16:03:58 +1300

DA bug fix (looks reasonably correct).

Diffstat:
M
src/Numeric/MCMC/NUTS.hs
 | 
9
++++-----
M
tests/Test.hs
 | 
2
+-

2 files changed, 5 insertions(+), 6 deletions(-)
diff --git a/src/Numeric/MCMC/NUTS.hs b/src/Numeric/MCMC/NUTS.hs
@@ -273,7 +273,7 @@ buildTreeDualAvg lTarget glTarget g t r logu v 0 e t0 r0 = do
       joint    = log $ auxilliaryTarget lTarget t1 r1
       n        = indicate (logu <= joint)
       s        = indicate (logu - 1000 <  joint)
-      a        = min 1 (acceptanceRatio lTarget t1 r1 t0 r0)
+      a        = min 1 (acceptanceRatio lTarget t0 t1 r0 r1)
   return (t1, r1, t1, r1, t1, n, s, a, 1)
       
 buildTreeDualAvg lTarget glTarget g t r logu v j e t0 r0 = do
@@ -320,14 +320,13 @@ findReasonableEpsilon lTarget glTarget t0 g = do
   let (t1, r1) = leapfrog glTarget (t0, r0) 1.0
       a        = 2 * indicate (acceptanceRatio lTarget t0 t1 r0 r1 > 0.5) - 1
 
-      go j e t r 
-        | j <= 0 = e -- no infinite loops
+      go e t r 
         | (acceptanceRatio lTarget t0 t r0 r) ^^ a > 2 ^^ (-a) = 
             let (tn, rn) = leapfrog glTarget (t, r) e
-            in  go (pred j) (2 ^^ a * e) tn rn 
+            in  go (2 ^^ a * e) tn rn 
         | otherwise = e
 
-  return $ go 1000 1.0 t1 r1
+  return $ go 1.0 t1 r1
 
 -- | Simulate a single step of Hamiltonian dynamics.
 leapfrog :: Gradient -> Particle -> Double -> Particle
diff --git a/tests/Test.hs b/tests/Test.hs
@@ -73,7 +73,7 @@ genMoves = replicateM 1000 genMove
 
 main = do
   test <- withSystemRandom . asGenIO $ 
-    nutsDualAveraging lTarget glTarget 100 10 t0
+    nutsDualAveraging lTarget glTarget 5000 1500 t0
     -- nuts lTarget glTarget 5000 0.1 t0 
     -- genMovesDa
     -- genMoves