mighty-metropolis

The classic Metropolis algorithm.
Log | Files | Refs | README | LICENSE

commit 918a241c6a13aa8e1f49f2b25cfe07c0d7ad7aec
parent 5499a94f481089dd11744f5cbad65e20a1cad69d
Author: Jared Tobin <jared@jtobin.io>
Date:   Thu, 21 May 2020 20:23:09 +0400

Put 'Spec' in with the other test modules.

Diffstat:
Mmighty-metropolis.cabal | 5++---
Atest/Spec.hs | 144+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Dtest/test/Spec.hs | 143-------------------------------------------------------------------------------
3 files changed, 146 insertions(+), 146 deletions(-)

diff --git a/mighty-metropolis.cabal b/mighty-metropolis.cabal @@ -73,7 +73,7 @@ Test-suite bnn Test-suite tests type: exitcode-stdio-1.0 - hs-source-dirs: test/test + hs-source-dirs: test main-is: Spec.hs default-language: Haskell2010 ghc-options: @@ -85,4 +85,4 @@ Test-suite tests , mwc-probability >= 1.0.1 , hspec , mwc-random - , mcmc-types -\ No newline at end of file + , mcmc-types diff --git a/test/Spec.hs b/test/Spec.hs @@ -0,0 +1,144 @@ +{-# OPTIONS_GHC -Wall #-} +{-# LANGUAGE RecordWildCards #-} + +import Data.Functor.Identity +import Data.Maybe (mapMaybe) +import Data.Sampling.Types +import Numeric.MCMC.Metropolis (chain, chain') +import System.Random.MWC +import Test.Hspec + +withinPercent :: Double -> Double -> Double -> Bool +withinPercent b n a + | b == 0 = a == 0 + | otherwise = d / b < n / 100 + where + d = abs (a - b) + +mean :: [Double] -> Double +mean xs = sum xs / n where + n = fromIntegral (length xs) + +variance :: [Double] -> Double +variance xs = sum [(x - m) ** 2.0 | x <- xs] / (n - 1) where + m = mean xs + n = fromIntegral (length xs) + +stdDev :: [Double] -> Double +stdDev = sqrt . variance + +stdErr :: [Double] -> Double +stdErr xs = stdDev xs / sqrt n where + n = fromIntegral (length xs) + +thin :: Int -> [a] -> [a] +thin n xs = case xs of + (h:t) -> h : thin n (drop (pred n) t) + _ -> mempty + +data Params = Params { + pepochs :: Int + , pradial :: Double + , porigin :: Identity Double + , ptunable :: Maybe (Identity Double -> Double) + , pltarget :: Identity Double -> Double + , pthin :: Int + } + +testParams :: Params +testParams = Params { + pepochs = 1000000 + , pradial = 0.2 + , porigin = Identity 1.0 + , ptunable = Just (\(Identity x) -> x ** 3.0) + , pltarget = \(Identity x) -> if x > 0 then negate x else negate 1 / 0 + , pthin = 1000 + } + +getChainResults :: IO [Double] +getChainResults = do + let Params {..} = testParams + + boxed <- withSystemRandom . asGenIO $ + chain pepochs pradial porigin pltarget + + let positions = fmap (runIdentity . chainPosition) boxed + pure (thin pthin positions) + +getTunableResults :: IO [Double] +getTunableResults = do + let Params {..} = testParams + + boxed <- withSystemRandom . asGenIO $ + chain' pepochs pradial porigin pltarget ptunable + + let positions = mapMaybe chainTunables boxed + pure (thin pthin positions) + +testWithinPercent :: SpecWith () +testWithinPercent = describe "withinPercent" $ + it "works as expected" $ do + 106 `shouldNotSatisfy` withinPercent 100 5 + 105 `shouldNotSatisfy` withinPercent 100 5 + 104 `shouldSatisfy` withinPercent 100 5 + 96 `shouldSatisfy` withinPercent 100 5 + 95 `shouldNotSatisfy` withinPercent 100 5 + 94 `shouldNotSatisfy` withinPercent 100 5 + +testMean :: SpecWith () +testMean = describe "mean" $ + it "works as expected" $ do + mean [1, 2, 3] `shouldSatisfy` withinPercent 2 1e-3 + mean [1..100] `shouldSatisfy` withinPercent 50.5 1e-3 + mean [1..1000000] `shouldSatisfy` withinPercent 500000.5 1e-3 + +testVariance :: SpecWith () +testVariance = describe "variance" $ + it "works as expected" $ do + variance [0, 1] `shouldSatisfy` withinPercent 0.5 1e-3 + variance [1, 1, 1] `shouldSatisfy` withinPercent 0 1e-3 + variance [1..100] `shouldSatisfy` withinPercent 841.66666666 1e-3 + +testStdErr :: SpecWith () +testStdErr = describe "stdErr" $ + it "works as expected" $ do + stdErr [1..100] `shouldSatisfy` withinPercent 2.901149 1e-3 + stdErr [1..1000] `shouldSatisfy` withinPercent 9.133273 1e-3 + +testHelperFunctions :: SpecWith () +testHelperFunctions = describe "helper functions" $ do + testWithinPercent + testMean + testVariance + testStdErr + +testSamples :: [Double] -> SpecWith () +testSamples xs = describe "sampled trace over exp(1)" $ do + let meanStdErr = stdErr xs + varStdErr = stdErr (fmap (\x -> pred x ** 2.0) xs) + + it "has the expected mean" $ do + mean xs `shouldSatisfy` (< 1 + 2 * meanStdErr) + mean xs `shouldSatisfy` (> 1 - 2 * meanStdErr) + + it "has the expected variance" $ do + variance xs `shouldSatisfy` (< 1 + 2 * varStdErr) + variance xs `shouldSatisfy` (> 1 - 2 * varStdErr) + +testTunables :: [Double] -> SpecWith () +testTunables ts = describe "sampled tunables over exp(1)" $ do + let meanStdErr = stdErr ts + + it "has the expected third moment (i.e. 6)" $ do + mean ts `shouldSatisfy` (< 6 + 2 * meanStdErr) + mean ts `shouldSatisfy` (> 6 - 2 * meanStdErr) + +main :: IO () +main = do + xs <- getChainResults + ts <- getTunableResults + + hspec $ do + testHelperFunctions + testSamples xs + testTunables ts diff --git a/test/test/Spec.hs b/test/test/Spec.hs @@ -1,143 +0,0 @@ -{-# OPTIONS_GHC -Wall #-} -{-# LANGUAGE RecordWildCards #-} - -import Data.Functor.Identity -import Data.Maybe (mapMaybe) -import Data.Sampling.Types -import Numeric.MCMC.Metropolis (chain, chain') -import System.Random.MWC -import Test.Hspec - -withinPercent :: Double -> Double -> Double -> Bool -withinPercent b n a - | b == 0 = a == 0 - | otherwise = d / b < n / 100 - where - d = abs (a - b) - -mean :: [Double] -> Double -mean xs = sum xs / n where - n = fromIntegral (length xs) - -variance :: [Double] -> Double -variance xs = sum [(x - m) ** 2.0 | x <- xs] / (n - 1) where - m = mean xs - n = fromIntegral (length xs) - -stdDev :: [Double] -> Double -stdDev = sqrt . variance - -stdErr :: [Double] -> Double -stdErr xs = stdDev xs / sqrt n where - n = fromIntegral (length xs) - -thin :: Int -> [a] -> [a] -thin n xs = case xs of - (h:t) -> h : thin n (drop (pred n) t) - _ -> mempty - -data Params = Params { - pepochs :: Int - , pradial :: Double - , porigin :: Identity Double - , ptunable :: Maybe (Identity Double -> Double) - , pltarget :: Identity Double -> Double - , pthin :: Int - } - -testParams :: Params -testParams = Params { - pepochs = 1000000 - , pradial = 0.2 - , porigin = Identity 1.0 - , ptunable = Just (\(Identity x) -> x ** 3.0) - , pltarget = \(Identity x) -> if x > 0 then negate x else negate 1 / 0 - , pthin = 1000 - } - -getChainResults :: IO [Double] -getChainResults = do - let Params {..} = testParams - - boxed <- withSystemRandom . asGenIO $ - chain pepochs pradial porigin pltarget - - let positions = fmap (runIdentity . chainPosition) boxed - pure (thin pthin positions) - -getTunableResults :: IO [Double] -getTunableResults = do - let Params {..} = testParams - - boxed <- withSystemRandom . asGenIO $ - chain' pepochs pradial porigin pltarget ptunable - - let positions = mapMaybe chainTunables boxed - pure (thin pthin positions) - -testWithinPercent :: SpecWith () -testWithinPercent = describe "withinPercent" $ - it "works as expected" $ do - 106 `shouldNotSatisfy` withinPercent 100 5 - 105 `shouldNotSatisfy` withinPercent 100 5 - 104 `shouldSatisfy` withinPercent 100 5 - 96 `shouldSatisfy` withinPercent 100 5 - 95 `shouldNotSatisfy` withinPercent 100 5 - 94 `shouldNotSatisfy` withinPercent 100 5 - -testMean :: SpecWith () -testMean = describe "mean" $ - it "works as expected" $ do - mean [1, 2, 3] `shouldSatisfy` withinPercent 2 1e-3 - mean [1..100] `shouldSatisfy` withinPercent 50.5 1e-3 - mean [1..1000000] `shouldSatisfy` withinPercent 500000.5 1e-3 - -testVariance :: SpecWith () -testVariance = describe "variance" $ - it "works as expected" $ do - variance [0, 1] `shouldSatisfy` withinPercent 0.5 1e-3 - variance [1, 1, 1] `shouldSatisfy` withinPercent 0 1e-3 - variance [1..100] `shouldSatisfy` withinPercent 841.66666666 1e-3 - -testStdErr :: SpecWith () -testStdErr = describe "stdErr" $ - it "works as expected" $ do - stdErr [1..100] `shouldSatisfy` withinPercent 2.901149 1e-3 - stdErr [1..1000] `shouldSatisfy` withinPercent 9.133273 1e-3 - -testHelperFunctions :: SpecWith () -testHelperFunctions = describe "helper functions" $ do - testWithinPercent - testMean - testStdErr - -testSamples :: [Double] -> SpecWith () -testSamples xs = describe "sampled trace over exp(1)" $ do - let meanStdErr = stdErr xs - varStdErr = stdErr (fmap (\x -> pred x ** 2.0) xs) - - it "has the expected mean" $ do - mean xs `shouldSatisfy` (< 1 + 2 * meanStdErr) - mean xs `shouldSatisfy` (> 1 - 2 * meanStdErr) - - it "has the expected variance" $ do - variance xs `shouldSatisfy` (< 1 + 2 * varStdErr) - variance xs `shouldSatisfy` (> 1 - 2 * varStdErr) - -testTunables :: [Double] -> SpecWith () -testTunables ts = describe "sampled tunables over exp(1)" $ do - let meanStdErr = stdErr ts - - it "has the expected third moment (i.e. 6)" $ do - mean ts `shouldSatisfy` (< 6 + 2 * meanStdErr) - mean ts `shouldSatisfy` (> 6 - 2 * meanStdErr) - -main :: IO () -main = do - xs <- getChainResults - ts <- getTunableResults - - hspec $ do - testHelperFunctions - testSamples xs - testTunables ts