hamiltonian seems to work!

cspollard · Apr 7, 2019 · e6c822a · e6c822a
1 parent 3f6840d
commit e6c822a
Show file tree

Hide file tree

Showing 2 changed files with 11 additions and 10 deletions.
diff --git a/ext/plot.py b/ext/plot.py
@@ -16,7 +16,8 @@
 xs_sorted = np.flipud(xs_unsorted[xs_unsorted[:,0].argsort()])
 
 bests = xs_sorted[0]
-xs = xs_unsorted.transpose()
+xs = xs_sorted.transpose()
+xs_unsorted = xs_unsorted.transpose()
 
 print xs.shape
 
@@ -60,7 +61,7 @@
     width = 0.7 * (bins[1] - bins[0])
     center = (bins[:-1] + bins[1:]) / 2
 
-    best = bests[i]
+    best = param[0]
     med = np.median(param)
     (q16, q84) = np.percentile(param, [16, 84])
     globq16 = np.min(best68)
@@ -108,7 +109,7 @@
     fig, axis = plt.subplots()
     fig.suptitle(name)
 
-    plt.plot(param)
+    plt.plot(xs_unsorted[i])
     plt.savefig("%s_t.pdf" % name)
 
     plt.clf()

diff --git a/src/MarkovChain.hs b/src/MarkovChain.hs
@@ -57,7 +57,7 @@ metropolisStep proposal logLH (currloc, currllh) = do
 
 
 hamiltonian
-  :: (PrimMonad m, Traversable f, Additive f)
+  :: (PrimMonad m, Traversable f, Additive f, Show (f Double))
   => Int                        -- ^ the number of leapfrog steps per sample
   -> Double                     -- ^ epsilon in the leapfrog steps
   -> (f Double -> Double)       -- ^ the potential energy function
@@ -68,7 +68,7 @@ hamiltonian steps eps u du =
 
 
 hamiltonianStep
-  :: (PrimMonad m, Traversable f, Additive f)
+  :: (PrimMonad m, Traversable f, Additive f, Show (f Double))
   => Int
   -> Double
   -> (f Double -> Double)
@@ -101,17 +101,17 @@ leapfrog
   -> (f a, f a)
   -> (f a, f a)
 leapfrog steps eps du (q0s, p0s) =
-  let p1s = dienan "p1s" <$> lerp (-eps/2) (du q0s) p0s
+  let p1s = dienan "p1s" <$> ((-eps/2) *^ du q0s ^+^ p0s)
 
       its =
         flip iterate (q0s, p1s) $ \(qs, ps) ->
-          let qs' = dienan "qs'" <$> lerp eps ps qs
-              ps' = dienan "ps'" <$> lerp (-eps) (du qs') ps
+          let qs' = dienan "qs'" <$> (eps *^ ps ^+^ qs)
+              ps' = dienan "ps'" <$> ((-eps) *^ du qs' ^+^ ps)
           in (qs', ps')
 
       (q2s, p2s) = its !! (steps-1)
-      qfs = lerp eps p2s q2s
-      pfs = lerp (-eps/2) (du qfs) p2s
+      qfs = eps *^ p2s ^+^ q2s
+      pfs = (-eps/2) *^ du qfs ^+^ p2s
 
   in (qfs, pfs)