commit 2748f973e951ffb348aa2c95222a744228ef5d87
parent aab096ca92619239cdd75cc0144063c9d4275cb2
Author: Jared Tobin <jared@jtobin.ca>
Date: Wed, 16 Oct 2013 10:09:27 +1300
Add examples. Remove old testing code.
Diffstat:
3 files changed, 56 insertions(+), 100 deletions(-)
diff --git a/src/Numeric/MCMC/Examples/Examples.hs b/src/Numeric/MCMC/Examples/Examples.hs
@@ -0,0 +1,56 @@
+-- Various examples, using NUTS with dual-averaging. Insert whatever trace
+-- (rosenbrockTrace, bnnTrace, etc.) you want into 'main' in order to spit out
+-- some observations.
+--
+-- A convenient R script to display these traces:
+--
+-- require(ggplot2)
+-- system('runhaskell Examples.hs > trace.csv')
+-- d = read.csv('../tests/trace.dat', header = F)
+-- ggplot(d, aes(V1, V2)) + geom_point(alpha = 0.05, col = 'darkblue')
+--
+
+import Numeric.AD
+import Numeric.MCMC.NUTS
+import System.Random.MWC
+
+logRosenbrock :: RealFloat a => [a] -> a
+logRosenbrock [x0, x1] = negate (5 * (x1 - x0 ^ 2) ^ 2 + 0.05 * (1 - x0) ^ 2)
+
+rosenbrockTrace :: IO [Parameters]
+rosenbrockTrace = withSystemRandom . asGenST $
+ nutsDualAveraging logRosenbrock (grad logRosenbrock) 10000 1000 [0.0, 0.0]
+
+logHimmelblau :: RealFloat a => [a] -> a
+logHimmelblau [x0, x1] = negate ((x0 ^ 2 + x1 - 11) ^ 2 + (x0 + x1 ^ 2 - 7) ^ 2)
+
+himmelblauTrace :: IO [Parameters]
+himmelblauTrace = withSystemRandom . asGenST $
+ nutsDualAveraging logHimmelblau (grad logHimmelblau) 10000 1000 [0.0, 0.0]
+
+logBnn :: RealFloat a => [a] -> a
+logBnn [x0, x1] = -0.5 * (x0 ^ 2 * x1 ^ 2 + x0 ^ 2 + x1 ^ 2 - 8 * x0 - 8 * x1)
+
+bnnTrace :: IO [Parameters]
+bnnTrace = withSystemRandom . asGenST $
+ nutsDualAveraging logBnn (grad logBnn) 10000 1000 [0.0, 0.0]
+
+logBeale :: RealFloat a => [a] -> a
+logBeale [x0, x1]
+ | and [x0 >= -4.5, x0 <= 4.5, x1 >= -4.5, x1 <= 4.5]
+ = negate $
+ (1.5 - x0 + x0 * x1) ^ 2
+ + (2.25 - x0 + x0 * x1 ^ 2) ^ 2
+ + (2.625 - x0 + x0 * x1 ^ 3) ^ 2
+ | otherwise = - (1 / 0)
+
+bealeTrace :: IO [Parameters]
+bealeTrace = withSystemRandom . asGenST $
+ nutsDualAveraging logBeale (grad logBeale) 10000 1000 [0.0, 0.0]
+
+printTrace :: Show a => [a] -> IO ()
+printTrace = mapM_ (putStrLn . filter (`notElem` "[]") . show)
+
+main :: IO ()
+main = bnnTrace >>= printTrace
+
diff --git a/src/Numeric/MCMC/Examples/Rosenbrock.hs b/src/Numeric/MCMC/Examples/Rosenbrock.hs
@@ -1,17 +0,0 @@
-import Numeric.AD
-import Numeric.MCMC.NUTS
-import System.Random.MWC
-
-lTarget :: RealFloat a => [a] -> a
-lTarget [x0, x1] = (-1) * (5 * (x1 - x0 ^ 2) ^ 2 + 0.05 * (1 - x0) ^ 2)
-
-glTarget :: [Double] -> [Double]
-glTarget = grad lTarget
-
-inits :: [Double]
-inits = [0.0, 0.0]
-
-epochs :: Int
-epochs = 100
-
-
diff --git a/tests/Test.hs b/tests/Test.hs
@@ -1,83 +0,0 @@
-import Control.Lens
-import Control.Monad
-import Control.Monad.Primitive
-import Data.Vector (singleton)
-import Numeric.AD
-import Numeric.MCMC.NUTS
-import System.Random.MWC
-
-import Debug.Trace
-
-lTarget :: RealFloat a => [a] -> a
-lTarget [x0, x1] = (-1) * (5 * (x1 - x0 ^ 2) ^ 2 + 0.05 * (1 - x0) ^ 2)
-
--- glTarget :: [Double] -> [Double]
--- glTarget = grad lTarget
-
-glTarget [x, y] =
- let dx = 20 * x * (y - x ^ 2) + 0.1 * (1 - x)
- dy = -10 * (y - x ^ 2)
- in [dx, dy]
-
-t0 = [0.0, 0.0] :: [Double]
-r0 = [0.0, 0.0] :: [Double]
-logu = -0.12840 :: Double
-v = -1 :: Double
-n = 5 :: Int
-m = 1
-madapt = 0
-e = 0.1 :: Double
-eAvg = 1
-h = 0
-
-t0da = [0.0, 0.0] :: [Double]
-r0da = [0.0, 0.0] :: [Double]
-
-runBuildTree :: PrimMonad m => Gen (PrimState m) -> m BuildTree
-runBuildTree g = do
- liftM BuildTree $ buildTree lTarget glTarget g t0 r0 logu v n e
-
-getThetas :: IO [Parameters]
-getThetas = replicateM 1000 getTheta
-
-getTheta :: IO Parameters
-getTheta = do
- bt <- withSystemRandom . asGenIO $ \g ->
- buildTree lTarget glTarget g t0 r0 logu v n e
- return $ bt^._5
-
-getThetasDa :: IO [Parameters]
-getThetasDa = replicateM 1000 getThetaDa
-
-getThetaDa :: IO Parameters
-getThetaDa = do
- bt <- withSystemRandom . asGenIO $ \g ->
- buildTreeDualAvg lTarget glTarget g t0 r0 logu v n e t0da r0da
- return $ bt^._5
-
-genMoveDa :: IO Parameters
-genMoveDa = withSystemRandom . asGenIO $ \g -> do
- eps <- findReasonableEpsilon lTarget glTarget t0 g
- let daParams = basicDualAveragingParameters eps madapt
- blah <- nutsKernelDualAvg lTarget glTarget eps eAvg h m daParams t0 g
- return $ blah^._4
-
-genMovesDa :: IO [Parameters]
-genMovesDa = replicateM 1000 genMoveDa
-
-genMove :: IO Parameters
-genMove = withSystemRandom . asGenIO $ nutsKernel lTarget glTarget e t0
-
-genMoves :: IO [Parameters]
-genMoves = replicateM 1000 genMove
-
-main = do
- test <- withSystemRandom . asGenIO $
- nutsDualAveraging lTarget glTarget 5000 1000 t0
- -- nuts lTarget glTarget 5000 0.1 t0
- -- genMovesDa
- -- genMoves
- mapM_ (putStrLn . filter (`notElem` "[]") . show) test
- -- mapM_ print test
-
-
