commit 5851ec0d67f136a085a38e80c4c2aa202dd43a77
Author: Jared Tobin <jared@jtobin.ca>
Date: Mon, 6 Mar 2017 09:26:57 +1300
Initial commit.
Diffstat:
25 files changed, 749 insertions(+), 0 deletions(-)
diff --git a/.gitignore b/.gitignore
@@ -0,0 +1,4 @@
+data
+*swp
+*swo
+.stack-work
diff --git a/LICENSE b/LICENSE
@@ -0,0 +1,19 @@
+Copyright (c) 2017 Jared Tobin
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in
+all copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+THE SOFTWARE.
diff --git a/README.md b/README.md
@@ -0,0 +1,82 @@
+
+# deanie
+
+*deanie* is an embedded probabilistic programming language. It can be used to
+denote, sample from, and perform inference on probabilistic programs.
+
+## Usage
+
+Programs are written in a straightforward monadic style:
+
+``` haskell
+mixture :: Double -> Double -> Program Double
+mixture a b = do
+ p <- beta a b
+ accept <- bernoulli p
+ if accept
+ then gaussian (negate 2) 0.5
+ else gaussian 2 0.5
+```
+
+You can sample from them by first converting them into an *RVar* from
+[random-fu][rafu]:
+
+```
+> sample (rvar (mixture 1 3))
+```
+
+Sample more using standard monadic combinators like 'replicateM':
+
+```
+> replicateM 1000 (sample (rvar (mixture 1 3)))
+```
+
+![](assets/mixture.png)
+
+You can perform inference on models using rejection or importance sampling, or
+use a simple, stateful Metropolis backend. Here's a simple beta-bernoulli model
+conditioned on the provided observations:
+
+``` haskell
+betaBernoulli :: Double -> Double -> Program Bool
+betaBernoulli a b = do
+ p <- beta a b
+ bernoulli p
+
+observations :: [Bool]
+observations = [True, True, False, True, False, False, True, True, True]
+```
+
+Here's one way to encode a posterior via rejection sampling:
+
+``` haskell
+rposterior :: Double -> Double -> Program Double
+rposterior a b =
+ grejection
+ (\xs ys -> count xs == count ys)
+ observations (beta a b) bernoulli
+ where
+ count = length . filter id
+```
+
+![](assets/bb_rejection.png)
+
+Here's another, via importance sampling:
+
+``` haskell
+iposterior :: Double -> Double -> Program (Double, Double)
+iposterior a b =
+ importance observations (beta a b) logDensityBernoulli
+```
+
+There are also some Monte Carlo convenience functions provided, such as a
+weighted average for weighted samples returned via importance sampling:
+
+```
+> samples <- replicateM 1000 (sample (rvar (iposterior 1 1)))
+> print (mcw samples)
+0.6369246537796793
+```
+
+[rafu]: https://hackage.haskell.org/package/random-fu
+
diff --git a/Setup.hs b/Setup.hs
@@ -0,0 +1,2 @@
+import Distribution.Simple
+main = defaultMain
diff --git a/assets/bb_rejection.png b/assets/bb_rejection.png
Binary files differ.
diff --git a/assets/mixture.png b/assets/mixture.png
Binary files differ.
diff --git a/deanie.cabal b/deanie.cabal
@@ -0,0 +1,91 @@
+name: deanie
+version: 0.1.0
+synopsis: An embedded probabilistic programming language.
+license: MIT
+license-file: LICENSE
+author: Jared Tobin
+maintainer: jared@jtobin.io
+build-type: Simple
+cabal-version: >=1.10
+description:
+ An embedded probabilistic programming language.
+
+library
+ default-language: Haskell2010
+ hs-source-dirs: lib
+ default-extensions:
+ LambdaCase
+ , RecordWildCards
+ other-modules:
+ Control.Applicative.Extended
+ exposed-modules:
+ Deanie
+ , Deanie.Density
+ , Deanie.Expr
+ , Deanie.Inference
+ , Deanie.Inference.Importance
+ , Deanie.Inference.Metropolis
+ , Deanie.Inference.Rejection
+ , Deanie.Language
+ , Deanie.RVar
+ build-depends:
+ base
+ , foldl
+ , free
+ , math-functions
+ , monad-loops
+ , random-fu
+ , recursion-schemes
+
+Test-suite betabernoulli-rejection
+ type: exitcode-stdio-1.0
+ hs-source-dirs: test
+ Main-is: BetaBernoulliRejection.hs
+ default-language: Haskell2010
+ ghc-options:
+ -O2 -Wall
+ other-modules:
+ Model.BetaBernoulli
+ build-depends:
+ base
+ , deanie
+
+Test-suite betabernoulli-importance
+ type: exitcode-stdio-1.0
+ hs-source-dirs: test
+ Main-is: BetaBernoulliImportance.hs
+ default-language: Haskell2010
+ ghc-options:
+ -O2 -Wall
+ other-modules:
+ Model.BetaBernoulli
+ build-depends:
+ base
+ , deanie
+
+Test-suite betabernoulli-metropolis
+ type: exitcode-stdio-1.0
+ hs-source-dirs: test
+ Main-is: BetaBernoulliMetropolis.hs
+ default-language: Haskell2010
+ ghc-options:
+ -O2 -Wall
+ other-modules:
+ Model.BetaBernoulli
+ build-depends:
+ base
+ , deanie
+
+Test-suite mixture-importance
+ type: exitcode-stdio-1.0
+ hs-source-dirs: test
+ Main-is: MixtureImportance.hs
+ default-language: Haskell2010
+ ghc-options:
+ -O2 -Wall
+ other-modules:
+ Model.Mixture
+ build-depends:
+ base
+ , deanie
+
diff --git a/lib/Control/Applicative/Extended.hs b/lib/Control/Applicative/Extended.hs
@@ -0,0 +1,12 @@
+
+module Control.Applicative.Extended (
+ replicateA
+
+ , module Control.Applicative
+ ) where
+
+import Control.Applicative
+
+replicateA :: Applicative f => Int -> f a -> f [a]
+replicateA n = sequenceA . replicate n
+
diff --git a/lib/Deanie.hs b/lib/Deanie.hs
@@ -0,0 +1,10 @@
+
+module Deanie (
+ module Export
+ ) where
+
+import Deanie.Language as Export
+import Deanie.Expr as Export
+import Deanie.Density as Export
+import Deanie.Inference as Export
+import Deanie.RVar as Export
diff --git a/lib/Deanie/Density.hs b/lib/Deanie/Density.hs
@@ -0,0 +1,41 @@
+
+module Deanie.Density (
+ logDensityBernoulli
+ , logDensityBeta
+ , logDensityGamma
+ , logDensityGaussian
+ ) where
+
+import Numeric.SpecFunctions (logGamma)
+
+logDensityBernoulli :: Double -> Bool -> Double
+logDensityBernoulli p x
+ | p < 0 || p > 1 = log 0
+ | x = log p
+ | otherwise = log (1 - p)
+{-# INLINE logDensityBernoulli #-}
+
+logDensityBeta :: Double -> Double -> Double -> Double
+logDensityBeta a b x
+ | x <= 0 || x >= 1 = log 0
+ | a < 0 || b < 0 = log 0
+ | otherwise = (a - 1) * log x + (b - 1) * log (1 - x)
+{-# INLINE logDensityBeta #-}
+
+-- | The log-density for the gamma distribution with rate parameter b.
+--
+-- f(x) = b ^ a * x ^ (a - 1) * x ^ (-b * x) / Gamma(a)
+logDensityGamma :: Double -> Double -> Double -> Double
+logDensityGamma a b x
+ | a < 0 || b < 0 = log 0
+ | otherwise = a * log b + (a - 1) * log x - b * x - logGamma a
+{-# INLINE logDensityGamma #-}
+
+-- | The density for the normal distribution with specified mean and standard
+-- deviation.
+logDensityGaussian :: Double -> Double -> Double -> Double
+logDensityGaussian m sd x
+ | sd <= 0 = log 0
+ | otherwise = negate (log sd) - (x - m) * (x - m) / (2 * sd * sd)
+{-# INLINE logDensityGaussian #-}
+
diff --git a/lib/Deanie/Expr.hs b/lib/Deanie/Expr.hs
@@ -0,0 +1,65 @@
+
+module Deanie.Expr (
+ -- * utilities
+
+ product
+ , iid
+ , indep
+
+ -- * distributions
+
+ , binomial
+ , uniform
+ , exponential
+ , chisq
+ , coin
+ , lgaussian
+ , invgamma
+ , geometric
+ ) where
+
+import Control.Applicative.Extended
+import qualified Control.Foldl as L
+import Control.Monad
+import Data.Function (fix)
+import Deanie.Language
+import Prelude hiding (product)
+
+product :: Ap (Free ProgramF) a -> Program a
+product term = liftF (ProgramF (InR term))
+
+iid :: Int -> Program a -> Program [a]
+iid n term = product (replicateA n (liftAp term))
+
+indep :: f a -> Ap f a
+indep = liftAp
+
+binomial :: Int -> Double -> Program Int
+binomial n p = fmap count (replicateM n (bernoulli p)) where
+ count = L.fold (L.handles (L.filtered id) L.length)
+
+uniform :: Program Double
+uniform = beta 1 1
+
+exponential :: Double -> Program Double
+exponential = gamma 1
+
+chisq :: Integral a => a -> Program Double
+chisq k = gamma (fromIntegral k / 2) (1 / 2)
+
+coin :: Program Bool
+coin = bernoulli 0.5
+
+lgaussian :: Double -> Double -> Program Double
+lgaussian m sd = fmap exp (gaussian m sd)
+
+invgamma :: Double -> Double -> Program Double
+invgamma a b = fmap recip (gamma a b)
+
+geometric :: Double -> Program Int
+geometric p = fix $ \count -> do
+ accept <- bernoulli p
+ if accept
+ then return 1
+ else fmap succ count
+
diff --git a/lib/Deanie/Inference.hs b/lib/Deanie/Inference.hs
@@ -0,0 +1,43 @@
+
+module Deanie.Inference (
+
+ -- * functions
+
+ mc
+ , mcw
+
+ -- * modules
+
+ , module Export
+ ) where
+
+import qualified Control.Foldl as L
+import Data.Monoid
+import Deanie.Inference.Importance as Export
+import Deanie.Inference.Metropolis as Export
+import Deanie.Inference.Rejection as Export
+
+data Average a = Average !a !a
+
+instance Fractional a => Monoid (Average a) where
+ mempty = Average 0 0
+ mappend (Average nL xL) (Average nR xR) = Average n x where
+ n = nL + nR
+ x = (xL * nL + xR * nR) / (nL + nR)
+
+average :: Fractional a => L.Fold a a
+average = L.Fold tally mempty summarize where
+ tally x a = x <> Average 1 a
+ summarize (Average _ x) = x
+{-# INLINE average #-}
+
+mc :: (Foldable f, Fractional a) => f a -> a
+mc = L.fold average
+{-# INLINE mc #-}
+
+mcw :: (Foldable f, Fractional a) => f (a, a) -> a
+mcw pts = L.fold (L.premap weight L.sum) pts where
+ mass = L.fold (L.premap fst L.sum) pts
+ weight (w, v) = w / mass * v
+{-# INLINE mcw #-}
+
diff --git a/lib/Deanie/Inference/Comonadic.hs b/lib/Deanie/Inference/Comonadic.hs
@@ -0,0 +1,22 @@
+
+module Deanie.Inference.Comonadic where
+
+import Control.Comonad.Cofree
+import Deanie.Language
+
+type Execution a = Cofree Program a
+
+initialize :: ProgramF a -> ()
+initialize = \case
+ ProgramF branch -> case branch of
+ InL term -> foo term
+ InR term -> ()
+ where
+ foo = \case
+ BernoulliF p _ -> ()
+ BetaF a b _ -> ()
+ GammaF a b _ -> ()
+ GaussianF m s _ -> ()
+
+
+
diff --git a/lib/Deanie/Inference/Importance.hs b/lib/Deanie/Inference/Importance.hs
@@ -0,0 +1,16 @@
+
+module Deanie.Inference.Importance (
+ importance
+ ) where
+
+import qualified Control.Foldl as L
+
+importance
+ :: (Foldable f, Monad m, Floating w)
+ => f a -> m b -> (b -> a -> w)
+ -> m (w, b)
+importance obs prior model = do
+ parameter <- prior
+ let cost = L.fold (L.premap (model parameter) L.sum) obs
+ return (exp cost, parameter)
+
diff --git a/lib/Deanie/Inference/Metropolis.hs b/lib/Deanie/Inference/Metropolis.hs
@@ -0,0 +1,47 @@
+
+module Deanie.Inference.Metropolis (
+ metropolis
+ ) where
+
+import Deanie.Language
+import qualified Control.Foldl as L
+import Control.Monad.Loops
+
+data MHP a = MHP {
+ n :: {-# UNPACK #-} !Int
+ , current :: !a
+ , ccost :: {-# UNPACK #-} !Double
+ }
+
+metropolis
+ :: Foldable f
+ => Int -> f a -> Program b -> (b -> a -> Double)
+ -> Program [b]
+metropolis epochs obs prior model = do
+ current <- prior
+ unfoldrM mh MHP { n = epochs, ccost = cost current, .. }
+ where
+ cost param = L.fold (L.premap (model param) L.sum) obs
+
+ mh MHP {..}
+ | n <= 0 = return Nothing
+ | otherwise = do
+ proposal <- prior
+ let pcost = cost proposal
+ prob = moveProbability ccost pcost
+ accept <- bernoulli prob
+ let (nepochs, nlocation, ncost) =
+ if accept
+ then (pred n, proposal, pcost)
+ else (pred n, current, ccost)
+ return (Just (nlocation, MHP nepochs nlocation ncost))
+
+moveProbability :: Double -> Double -> Double
+moveProbability current proposal =
+ whenNaN 0 (exp (min 0 (proposal - current)))
+ where
+ whenNaN val x
+ | isNaN x = val
+ | otherwise = x
+{-# INLINE moveProbability #-}
+
diff --git a/lib/Deanie/Inference/Rejection.hs b/lib/Deanie/Inference/Rejection.hs
@@ -0,0 +1,25 @@
+
+module Deanie.Inference.Rejection (
+
+ rejection
+ , grejection
+ ) where
+
+import Control.Monad
+import qualified Data.Foldable as F
+
+grejection
+ :: (Foldable f, Monad m)
+ => ([a] -> [b] -> Bool) -> f b -> m c -> (c -> m a) -> m c
+grejection predicate observed proposal model = loop where
+ len = length observed
+ loop = do
+ parameters <- proposal
+ generated <- replicateM len (model parameters)
+ if predicate generated (F.toList observed)
+ then return parameters
+ else loop
+
+rejection :: (Foldable f, Monad m, Eq a) => f a -> m b -> (b -> m a) -> m b
+rejection = grejection (==)
+
diff --git a/lib/Deanie/Language.hs b/lib/Deanie/Language.hs
@@ -0,0 +1,57 @@
+{-# LANGUAGE DeriveFunctor #-}
+
+module Deanie.Language (
+ -- * types
+
+ ProbF(..)
+ , Prob
+ , ProgramF(..)
+ , Program
+
+ -- * terms
+
+ , bernoulli
+ , beta
+ , gaussian
+ , gamma
+
+ -- * re-export
+
+ , module Control.Applicative.Free
+ , module Control.Monad.Free
+ , module Data.Functor.Sum
+
+ ) where
+
+import Control.Applicative.Extended (replicateA)
+import Control.Applicative.Free hiding (Pure)
+import Control.Monad.Free
+import Data.Functor.Sum
+import Prelude hiding (product)
+
+data ProbF r =
+ BernoulliF {-# UNPACK #-} !Double (Bool -> r)
+ | BetaF {-# UNPACK #-} !Double {-# UNPACK #-} !Double (Double -> r)
+ | GaussianF {-# UNPACK #-} !Double {-# UNPACK #-} !Double (Double -> r)
+ | GammaF {-# UNPACK #-} !Double {-# UNPACK #-} !Double (Double -> r)
+ deriving Functor
+
+type Prob = Free ProbF
+
+newtype ProgramF a = ProgramF (Sum ProbF (Ap (Free ProgramF)) a)
+ deriving Functor
+
+type Program = Free ProgramF
+
+bernoulli :: Double -> Program Bool
+bernoulli p = liftF (ProgramF (InL (BernoulliF p id)))
+
+beta :: Double -> Double -> Program Double
+beta a b = liftF (ProgramF (InL (BetaF a b id)))
+
+gamma :: Double -> Double -> Program Double
+gamma a b = liftF (ProgramF (InL (GammaF a b id)))
+
+gaussian :: Double -> Double -> Program Double
+gaussian m s = liftF (ProgramF (InL (GaussianF m s id)))
+
diff --git a/lib/Deanie/RVar.hs b/lib/Deanie/RVar.hs
@@ -0,0 +1,28 @@
+
+module Deanie.RVar (
+ rvar
+
+ -- * re-exported
+
+ , sample
+ ) where
+
+import Control.Monad
+import Deanie.Language
+import Data.Random hiding (rvar)
+import qualified Data.Random.Distribution.Bernoulli as RF
+import qualified Data.Random.Distribution.Beta as RF
+import qualified Data.Random.Distribution.Gamma as RF
+import qualified Data.Random.Distribution.Normal as RF
+
+rvar :: Program a -> RVar a
+rvar = iterM $ \case
+ ProgramF (InL term) -> evalAlg term
+ ProgramF (InR term) -> join (runAp rvar term)
+ where
+ evalAlg = \case
+ BernoulliF p k -> RF.bernoulli p >>= k
+ BetaF a b k -> RF.beta a b >>= k
+ GammaF a b k -> RF.gamma a b >>= k
+ GaussianF m s k -> RF.normal m s >>= k
+
diff --git a/stack.yaml b/stack.yaml
@@ -0,0 +1,66 @@
+# This file was automatically generated by 'stack init'
+#
+# Some commonly used options have been documented as comments in this file.
+# For advanced use and comprehensive documentation of the format, please see:
+# http://docs.haskellstack.org/en/stable/yaml_configuration/
+
+# Resolver to choose a 'specific' stackage snapshot or a compiler version.
+# A snapshot resolver dictates the compiler version and the set of packages
+# to be used for project dependencies. For example:
+#
+# resolver: lts-3.5
+# resolver: nightly-2015-09-21
+# resolver: ghc-7.10.2
+# resolver: ghcjs-0.1.0_ghc-7.10.2
+# resolver:
+# name: custom-snapshot
+# location: "./custom-snapshot.yaml"
+resolver: lts-8.3
+
+# User packages to be built.
+# Various formats can be used as shown in the example below.
+#
+# packages:
+# - some-directory
+# - https://example.com/foo/bar/baz-0.0.2.tar.gz
+# - location:
+# git: https://github.com/commercialhaskell/stack.git
+# commit: e7b331f14bcffb8367cd58fbfc8b40ec7642100a
+# - location: https://github.com/commercialhaskell/stack/commit/e7b331f14bcffb8367cd58fbfc8b40ec7642100a
+# extra-dep: true
+# subdirs:
+# - auto-update
+# - wai
+#
+# A package marked 'extra-dep: true' will only be built if demanded by a
+# non-dependency (i.e. a user package), and its test suites and benchmarks
+# will not be run. This is useful for tweaking upstream packages.
+packages:
+- '.'
+# Dependency packages to be pulled from upstream that are not in the resolver
+# (e.g., acme-missiles-0.3)
+extra-deps: []
+
+# Override default flag values for local packages and extra-deps
+flags: {}
+
+# Extra package databases containing global packages
+extra-package-dbs: []
+
+# Control whether we use the GHC we find on the path
+# system-ghc: true
+#
+# Require a specific version of stack, using version ranges
+# require-stack-version: -any # Default
+# require-stack-version: ">=1.3"
+#
+# Override the architecture used by stack, especially useful on Windows
+# arch: i386
+# arch: x86_64
+#
+# Extra directories used by stack for building
+# extra-include-dirs: [/path/to/dir]
+# extra-lib-dirs: [/path/to/dir]
+#
+# Allow a newer minor version of GHC than the snapshot specifies
+# compiler-check: newer-minor
+\ No newline at end of file
diff --git a/test/BetaBernoulliImportance.hs b/test/BetaBernoulliImportance.hs
@@ -0,0 +1,19 @@
+
+module Main where
+
+import Control.Monad
+import Deanie
+import Model.BetaBernoulli as Model
+
+posterior :: Double -> Double -> Program (Double, Double)
+posterior a b =
+ importance Model.observations (beta a b) logDensityBernoulli
+
+main :: IO ()
+main = do
+ samples <- replicateM 1000 (sample (rvar (posterior 1 1)))
+ mapM_ print samples
+
+ putStrLn "\nestimate:"
+ print (mcw samples)
+
diff --git a/test/BetaBernoulliMetropolis.hs b/test/BetaBernoulliMetropolis.hs
@@ -0,0 +1,18 @@
+
+module Main where
+
+import Deanie
+import Model.BetaBernoulli as Model
+
+posterior :: Double -> Double -> Program [Double]
+posterior a b =
+ metropolis 1000 Model.observations (beta a b) logDensityBernoulli
+
+main :: IO ()
+main = do
+ samples <- sample (rvar (posterior 1 1))
+ mapM_ print samples
+
+ putStrLn "\nestimate:"
+ print (mc samples)
+
diff --git a/test/BetaBernoulliRejection.hs b/test/BetaBernoulliRejection.hs
@@ -0,0 +1,23 @@
+
+module Main where
+
+import Control.Monad
+import Deanie
+import Model.BetaBernoulli as Model
+
+posterior :: Double -> Double -> Program Double
+posterior a b =
+ grejection
+ (\xs ys -> count xs == count ys)
+ Model.observations (beta a b) bernoulli
+ where
+ count = length . filter id
+
+main :: IO ()
+main = do
+ samples <- replicateM 1000 (sample (rvar (posterior 1 1)))
+ mapM_ print samples
+
+ putStrLn "\nestimate:"
+ print (mc samples)
+
diff --git a/test/MixtureImportance.hs b/test/MixtureImportance.hs
@@ -0,0 +1,23 @@
+
+module Main where
+
+import Control.Monad
+import Deanie
+import Model.Mixture as Model
+
+logDensity :: Bool -> Double -> Double
+logDensity accept x
+ | accept = logDensityGaussian (negate 2) 0.5 x
+ | otherwise = logDensityGaussian 2 0.5 x
+
+posterior :: Program (Double, Bool)
+posterior = importance Model.observations (bernoulli 0.5) logDensity
+
+main :: IO ()
+main = do
+ samples <- replicateM 1000 (sample (rvar posterior))
+ mapM_ print samples
+
+ putStrLn "\nestimate:"
+ print (mcw (fmap (\(w, x) -> if x then (w, 1) else (w, 0)) samples))
+
diff --git a/test/Model/BetaBernoulli.hs b/test/Model/BetaBernoulli.hs
@@ -0,0 +1,16 @@
+
+module Model.BetaBernoulli (
+ betaBernoulli
+ , observations
+ ) where
+
+import Deanie
+
+betaBernoulli :: Double -> Double -> Program Bool
+betaBernoulli a b = do
+ p <- beta a b
+ bernoulli p
+
+observations :: [Bool]
+observations = [True, True, False, True, False, False, True, True, True]
+
diff --git a/test/Model/Mixture.hs b/test/Model/Mixture.hs
@@ -0,0 +1,19 @@
+
+module Model.Mixture (
+ mixture
+ , observations
+ ) where
+
+import Deanie
+
+mixture :: Double -> Double -> Program Double
+mixture a b = do
+ p <- beta a b
+ accept <- bernoulli p
+ if accept
+ then gaussian (negate 2) 0.5
+ else gaussian 2 0.5
+
+observations :: [Double]
+observations = [2.040498171768912,-1.4886800830664786,-2.2799168958984994,2.70438820866698,-2.019455618378096,2.542849699098371,1.3229918636104645,2.706715209209304,1.8811654040228547,2.350066155719298,-1.6047105477552586,1.615059226984037,1.8186851411137601,2.795297640244097,1.1168580729983073,-2.5920126754577066,-1.7885130341578905,2.6131218047808598,2.5812659415124797,2.3444242581367543,1.8645787746906526,2.032623371516321,1.4750402034595242,1.9556832053978666,1.48059807802758,2.068505992690237,1.6506110471451743,1.3006937877860025,-2.279825280993058,2.1501214015946157,2.514875162098206,2.2798785151277814,1.8027764410570382,-1.1167007133844695,2.800256264234736,-1.8198287602319874,-2.59931017506766,-1.6962609436891334,1.67663431336761,2.165842255227899,2.395793048292388,2.1615588717820886,2.049174438171109,1.947005287858322,1.6275033387698807,2.2939457358798756,1.5337043499515657,-2.1486603604944965,1.8624552341364018,1.9043877285273967,1.6670012786022568,2.3647792828000522,2.1650470885301236,1.3479390949501253,-1.7343052275581954,2.3179890069179723,-2.3162833991711236,2.4245220955809352,1.7494025814524208,-2.7111409574216996,-2.0024464199852217,-2.5636008775177284,2.4323929127660153,1.7101552901532557,-1.791164373611174,2.2427063921525283,2.465029060202423,1.743457323083459,1.4889768518700142,2.563421044085401,1.9106590571335695,-2.191307223089777,-3.3734326866988473,1.6404348623061318,2.070238335325899,1.8621420738658259,2.0054816238498425,2.659008825422621,2.4130497952997976,1.9428792631052152,-1.40472531559455,-2.821845794286177,2.3407130607846867,1.820387850984587,1.2950412106026628,-1.6421955632983278,1.617765736025464,-1.2251499331656295,-2.894589401411516,1.3249089780073369,-2.373796533503496,1.9469700677211814,2.1280250783510786,2.5246382238038785,1.695334875219744,1.6375722630131602,-1.918986263422856,1.7507547165317987,-3.3570645876201906,1.7967988598136233,1.9779700764086297,-1.8437070406163127,-1.9628147997404144,1.7021428722027625,1.6392796170905743,-2.2437768340122815,2.265129112943273,2.619915915473057,-2.457343228796371,2.096063867491764,1.164289627036515,2.2516480810969592,1.7010063089026983,2.22128576584823,1.2985689832254161,2.5387948661605595,2.7802167972272516,1.3914386060155521,1.6604899223036573,2.3201038277629533,2.8089974076957245,2.1249558435938494,1.7823961062480584,1.9379711257810657,1.7611594566002697,1.9570262296126966,2.957494454582769,2.0044354414321437,1.5407729229498404,-2.059194600981957,2.7552988100528,-1.960578535013273,2.4197252493237293,-2.5071214255060457,-2.3395615902640463,1.5723355725674484,2.5473477572458343,2.229775829378868,2.6556075869414792,2.635618524164154,1.959478273197359,2.0236817019605975,1.515259839656796,-2.0901323321012026,2.6296078287045357,1.731309314122821,2.350367671068628,1.1532998064912103,2.201152928057661,2.2426327780083954,2.624526663634715,2.043952189217419,-2.0589612195513958,1.9219343091509826,1.3240511031708453,1.878277565486727,2.5559712158418817,2.00825165671651,1.1144880526444845,-3.1435993237229445,2.0308742970614864,2.913723210397535,2.5033833802193923,1.1799244608290687,2.6572411422657325,1.4356635780893832,-2.4827727036772163,1.6377435979519928,-2.2793508977143913,2.3475760879462397,-1.2312566467674602,2.359934650546971,2.1118250757663946,-1.4062962387329447,2.08417690405576,1.3182807052533105,1.9864667409137695,-1.4182234363298254,-2.301669389320769,1.6843101822502946,2.0662291718538532,-2.0780960483321556,2.5691557652702124,2.1306764425119615,2.022654684177645,2.4628595864500915,2.0771417723163212,2.257807110553126,-2.1055314415608266,3.1495267960397575,1.9704459546686752,1.0351336893344747,1.5251858120897697,2.1148005853801157,-1.5491722313898224,1.6529036406181608,-2.047979384798056,1.4165978362194855,2.329631293476734,1.8688347278108233,1.7454270787629276,0.9153711664143733,2.182260377729185,-2.6622344676326937,2.2436976150656083,1.9940763091883458,-2.879724503644798,2.25217437212981,-1.8555922009428503,1.7405495892553904,2.248159180861862,-2.6792731076136533,3.0813081392448938,3.0432839813003394,1.8117896193546552,1.5981087947435981,1.440997265251592,2.1397232745974617,2.0973084071635792,-2.4806999516061126,-2.5831296822699628,1.8900637398158693,2.5330083157454983,1.8704560730328534,-1.848536371471426,1.3922044026905112,-1.9046141545487645,2.01787392140428,2.2000743556688542,2.700673982648665,1.6065343056442205,1.7376792667212215,2.2921140594391165,2.108079059243912,-2.1863108506646594,-3.1143160041828724,1.6744511148451917,2.299610738814311,-2.541820192337976,1.4852713757594227,1.0494011722510466,1.714472376406449,-2.0803360840922056,-1.2681960616841286,-1.4596858051649297,2.341520181250461,2.8722405346186632,2.389609346304171,2.070056803219863,1.5993790675799187,1.6231590739912507,2.9158526317692584,3.2703330696762962,-1.6083804308184102,2.1987177590300973,2.423016394418599,1.344882706182093,1.9543760819301816,1.88777265556394,-2.299033849221305,-1.914108941823389,1.8022748205025887,2.3179708797056247,-1.9717809517036262,2.5151545212695274,2.4959787693010163,2.672082864566547,2.132159635034685,1.9532667525922272,1.6510175535597846,2.1922149142815672,2.058673116400376,-2.3487757479170837,1.8192517017362597,1.8028076988577066,1.7833763401123128,1.863854364608229,1.9099584445545947,-1.5678722755798093,-1.558831216854349,1.441001770375314,-1.7965490909055246,1.1418721494170565,-2.0806295902242518,1.9602265359137412,2.6086560740403533,1.5128934988640932,2.31268099405124,2.3597162547164006,2.179418498263955,1.6864696771405991,1.934784522340504,1.482256499834551,2.651838425205662,1.7436448562818507,1.0075253605279166,2.670564162525278,2.2815719870675757,-3.909905695967328,1.3933498180879733,-1.9247794469694923,2.693515132286159,-1.6235089275375247,1.8834714469805123,2.0501127949290754,-1.5092968277247518,1.3533668840636506,-2.0365463610942114,2.641373963309193,-0.858807604930317,1.7626385996672458,2.1831963424890475,2.7484638980741525,-2.029583733457625,1.6532331920844667,1.482950186935497,2.432491636140653,0.9350457507196992,-2.3175443288422803,2.8841414010210107,1.6375900059832742,2.582745171770648,-1.6726351899449747,1.9031781325851966,-1.0268294866582255,-2.0721877688480927,-1.6592914910666052,1.9478946321814945,1.787711659317469,1.8305175333087327,2.355944609422427,2.1058872209084196,2.2379746183326406,2.1328007355432237,1.9587889843127044,-1.9897577860635147,1.5605872033839367,1.6236165478887863,1.2166343228438667,2.0005894065714194,1.6286842963769932,2.5008547597723854,-1.8601612104577685,-1.7634916858774967,1.9650365863587564,2.402134091605931,-2.2241722279043428,-2.7014322188932596,2.1283993710708717,-1.5940913068078677,-2.0227719851557846,2.695330678799133,2.0436146198470544,1.4585989193788431,1.281114565098774,2.3869426264495286,2.219061256348053,1.9076427573792785,-2.0415753607743623,2.4262363309930626,1.4753908300688812,-2.106963802772042,2.44933680607737,1.8980387085908987,2.226559005040951,2.1320406706486876,-1.678847182408624,2.066891735798001,2.562157490975004,1.454326217682123,1.9663958969619404,-1.0304412611828822,2.324276819209477,-1.3158809329511243,-1.2430336708977592,-1.954157028342249,1.761391132092033,-2.0667102344277892,3.0635648021542643,1.692985569020018,2.4770396462140956,-1.6084468052368541,2.174972273048207,2.2186631064029654,2.8718249053749765,1.6021949731605176,1.9673044237705235,1.773505663528724,-1.4527705820240195,1.9219452691476364,1.793026189230678,1.4131303157413986,-2.316955569212823,-3.5906442170025668,2.536798006181261,1.9988125056542612,2.0813528641534353,-2.4853456439408443,1.7084427187079603,1.8668415248634624,1.6592025390364724,2.291871258649108,-2.1197473368275976,2.9557928366486204,-2.048742720084168,-2.592053557079969,-1.046522390647616,-1.721403993368058,-2.429309873580671,1.1624847519097554,1.9800200060180877,2.1087971859167336,2.066365824590253,1.9121828364368474,2.0498662566050623,-1.1752572866411113,1.6009769136988559,-2.339616717274173,1.5237884887127335,2.9491527600644565,1.979167635273394,2.3432902490676804,2.034120331356947,-1.5419136869565637,2.682939153854116,1.9574425828071598,2.363264361740107,2.432871757210849,1.3726679680910974,-2.0937821014355174,2.161056431473673,1.399858107575219,1.3841693134410535,-2.6146400806155583,2.2903337013596516,2.231330018575711,2.3787021043106105,2.3919169159759344,2.3113155069093403,2.202146103148283,2.6412106548446577,2.078811157394137,1.18996641481856,2.8633028083643763,2.253163797834183,0.5189811495418195,2.3059768510490573,2.6149590712851336,1.6973767621783038,2.090938707828445,2.3357575211943953,1.4695762906747831,2.3036398241380134,-1.5754652217639977,2.0993273992195864,0.9466641565527283,-1.1263350091263185,2.5300256944162034,1.4116235839805276,3.0870736818446494,-1.5870325064096151,1.9123411324869086,-2.7271190721513183,1.5556107823809753,1.2572609458313397,2.6460775998218113,1.5109794411395285,2.0672231709904394,1.3717759584916447,1.711874513985119,-1.85670220997989,2.5470990733693526,2.189279118001529,1.4841172833653484,-2.523002321434689,-1.0367964919150356,2.168309256748491,1.2023089712551687,2.113479420836238,2.598308507073872,1.762605112970926,2.3461147031974936,2.778851545319778,2.126126336541477,-3.339397200949767,2.4881168695506855,2.6235132494429867,1.6576013392841396,1.4193311873485923,1.118021235216334,1.7631115702439668,1.3419696406410755,3.6079833994711485,1.8650256418343358,2.095667694484698,2.393871570672834,2.2723708375645506,1.526896650065234,2.260682770614744,-1.2945186786379965,1.3427113641844148,-2.357812228905531,1.6690855511761395,1.3850000102386424,1.756740361638026,2.3652285949760667,-1.3986035981882337,2.251363990429474,1.005001183423495,-2.115708201111327,1.0462396523181856,2.322707615587651,-1.6027262198666221,1.8186341608993195,1.9885741484289874,2.1830660281953502,-1.3433056793979041,-1.9442854487670562,1.300641333686504,2.508491523149976,1.8172312299689686,2.04405377217793,2.4570229864548154,1.3315595656721255,1.873487723269466,1.5225633925746536,-1.6825623796410034,1.6828935575717345,1.2772072473651068,1.875268704006719,3.1969042780195345,1.7613090711997532,2.1162258535030967,2.285948323239134,1.8262374852238366,2.208852987358055,-1.0364979573429385,3.026041213968913,1.716197037722503,1.7437900521962342,-1.967533085603372,2.4594060565558404,2.5367590895736214,2.2330133079174135,-1.2911589232829472,2.3037073706184477,1.7326885633829776,1.8740922464691778,1.6894049290137167,2.1127950576185306,-1.9103041727723964,1.4217797823476483,1.7281973533411792,-2.31101636477957,-1.4810989483569528,1.7735345614749654,2.4547733664644498,1.4382437155904253,-2.2926175218317213,2.1308609253866178,-2.585378602055923,-1.9837903014689091,-1.7420511755911894,1.6623543176239164,1.9616253950979516,2.185027850593901,2.1927943940059826,1.7376529008138482,2.860988373206612,1.4619444786561955,2.450699240222887,2.00662418205578,1.686138188650883,-1.7346787112656454,-2.1999868420601847,0.9132883953387188,1.9088820291893436,-1.9788466846069328,-1.4965895386559098,1.9073957174036509,2.029255162121264,2.66312410021487,1.4982798169301876,2.0469123058767074,2.6567761753901866,1.8932137778465519,-1.8556950124659686,-1.2989460440956984,2.333149864933306,-2.345094901729416,-2.4576900154940775,1.5726335673669063,-1.5927025483589716,2.229276962078651,0.9296083665915349,1.932895627794748,-2.5787147023831243,-1.829468854853467,2.1653053235258906,2.358404401653845,-2.011408977348583,3.015070210961185,2.3061022194405307,2.001749872096849,-1.652147091235575,3.8136198782163877,2.5713816273850614,1.4269332627813631,1.7603959898151311,-2.136655060700856,0.8548124473931684,1.9994057795510394,2.087806591176813,2.87221611584285,-2.527538581779394,2.777451852464437,3.0343702918564164,2.590240695447608,-2.0402993581854103,2.204850369181398,2.4015964685791107,1.8454226579586788,2.1129432422629018,1.0314838006302294,2.209147299497288,2.5438572232647223,2.2766937756937,3.1341053844184783,-1.235962113798965,2.4251691678829506,1.7265176647953555,2.5691116575484347,0.8158685705101785,-2.1866124782261567,1.2228130728722553,0.7326368186187586,1.490023572939612,1.7128969376145018,-1.6023910752405328,-2.082680059155702,-1.5703245247387845,0.8402728649077125,1.5411093734470631,1.2405125994706352,3.1614264263759058,1.259893895706659,2.096437899447434,1.4390812696346407,-2.1983546950788106,-1.7199046053078082,2.061118777134465,1.8570297003474736,-2.5446098462720563,2.2578660179948726,2.4744993725513376,2.4152815322125556,2.445755282692375,2.374213876589887,2.30740416465824,-1.3028401431459171,1.8360286708504276,1.7514883640947967,1.5783844488672922,1.8227002477771062,-1.478947591511421,1.9762563370858397,1.5705311745714037,1.923949720752053,2.19369805311981,1.1745193952707669,2.128097155475538,-1.98968010000829,1.9997220976147707,2.009501428269911,2.039548991068272,1.8694964348030332,1.9983979266002834,2.5037384196163246,2.0859981771330216,-2.532874776003714,-1.0334222467456966,1.4282799647922042,-1.4692828421676314,1.9066061528937914,2.643638708494721,1.5165453018749102,2.578453202105292,2.1409169379018325,1.7632217182960521,2.1758905710387197,-1.9351632754637869,1.8166105456695383,1.3340087345893707,1.126775070973284,1.2893943156118923,1.0090672046353701,2.067323985889357,2.2505033938132,1.307937547901393,-2.2842900456346373,2.5094891625607323,-1.5641593797382474,2.075874415235858,1.5344035355910859,2.595100953072377,-1.8728522004360657,1.4996156189373466,-1.7928861431459473,2.099139271753124,-2.194780505060433,2.7245783510067767,1.2217251349862672,2.5249052882797955,2.1508323157846445,1.4294734859052076,2.1731267982334104,1.007409019719604,-1.9443935935043615,2.0784681419094193,-2.419454321205,2.2634446100873076,1.4196631000830844,2.3463244367982554,1.9804154708053106,2.443045524433445,-2.492744089595733,1.1811801214678432,-1.416246955500415,1.8321290140170443,0.9669220527801639,1.7774649656725119,-1.6836846791280462,2.1621817852579035,1.7069605519240716,1.9698221692340858,1.0942333234736288,2.9129926941648407,1.8848897488113767,1.8048371813925534,2.364320917871216,2.7756931628567507,1.3401141601285858,-1.1537217558110893,1.9514627232145316,1.365089424060705,2.4486492707233394,-2.5768970802215745,2.0262544410232133,2.4733028483724935,1.7944646886772233,2.102556210962225,2.4982135646794874,-2.203014695288726,2.069908490669035,-1.6955181881196144,2.5727187035393517,2.275034768753278,-0.9446904847715745,2.581271196475106,1.7131757935387626,2.8931238437083193,-1.495498048630974,1.468522551217743,-2.8259916410092667,1.4571362201027258,3.0528368337526723,-2.4601841085468834,2.1528609085798935,-1.968827555114081,1.8407856180017623,2.1223956645615156,1.9201873339364806,2.9265867344905683,1.9672268432657738,1.9609139926294126,1.8919985592100006,-1.526258943618376,-2.188569365809944,2.067601679977788,2.1319678720789157,2.0290575719627846,2.1334130377073173,2.0691508006219284,2.1310579099102225,-2.227227957972257,2.0705795347560376,1.5243705477146778,1.4053544567146976,1.7492091509903902,-2.178111376257619,-2.001589181152252,-2.4285318417771222,1.5399644493751738,1.637916308378858,1.9863718598079345,-0.9153730840511975,1.6491606333977904,2.2445137332081635,1.1431832068177736,2.408461851871155,1.6551808758439521,1.74318866926646,-1.389387251811931,2.0423831835585937,-2.0492730929749094,-2.4048723408798462,-2.1766990151666064,2.2343184231816187,2.8468354358557866,1.6883265913969585,3.0727856956077337,2.4905735544516165,-1.4383709484589944,1.4354830759656845,-2.4431082115345046,-1.3642260439848533,-3.0871603949004385,1.3249669459231785,-2.776549245017973,-1.3801650509905306,2.531254655404986,2.204873962682287,1.3499027216356332,2.3513242633771263,1.9123963297216513,1.116341799875221,-2.9099693595467646,-1.6782690820741593,-2.36732152808579,1.782269011658509,2.5812203281209927,2.161892220370811,1.7822024647789407,-1.3461452916889427,1.2297738746441769,2.22681597807914,1.5808129580184467,2.2642524601368414,1.9819655286543434,-2.464734312981304,1.518599689647127,-2.4253063524430645,2.050286117643502,2.3156603144194245,1.9957550308541947,-1.9660259317117732,2.184403343522855,-2.5290337947244743,1.1931268691222776,-2.723119966757685,1.646996513728255,-2.2313925714156295,2.131695560477424,3.3989086630710608,2.7224780984467634,1.985477907829886,2.281228237687391,1.9415784065364088,2.0818210026106034,-1.2154025087170544,2.8160585568496463,-1.225223210594841,2.8078699500064035,-2.336244886373845,1.4251537493897675,3.0397538263505828,1.650303365371757,2.126927978739912,2.291599344312093,1.5051584786039152,-2.301762509172051,-1.8542613403661976,-2.1442539278633794,1.2884282760298258,1.548408060785982,2.1990106955473094,-2.435108586883474,-2.24110767240865,2.559503390765469,-2.126706690854664,1.4342064935797612,2.935149572783921,2.6703734191151756,2.042629948886332,-1.6569209603344839,2.631903759058977,1.5992532470473462,2.442349644015136,1.8981614211662798,2.728610354449893,-2.3885121892674963,2.221014071455413,2.3926745596774825,2.120005946219882,2.346175154991238,2.432412924754965,2.0383481499505574,-2.053485966850169,2.910707366530109,-3.1451800667915397,2.2997526429410367,2.2349354962856456,1.2259580818490494,3.297205264220768,2.193676807449321,1.8502197713835435,1.2813421496272577,1.9417259494012067,1.3166843975519882,1.2489535444642623,2.298333484804063,-1.558659896726138,-1.5447559659888432,1.5367383576456526,-1.4986587745135551,1.8375845134748738,-2.8403664056290383,2.2539723839423895,-2.579981795048657,-2.0361635508732836,2.3362402316715984,-1.549467453209816,1.6762537073744213,2.43134827814065,1.960221633685193,1.4161489168129244,-2.192863989552737,1.3675671285082402,2.194624480710983,1.9947933115199687,1.7207005912450635,1.5362410225434522,1.9378469881067169,0.8078920633651796,2.0605019652306864,1.8217101122967445,1.3158059101581514,2.4222250082330654,1.7139001185747207,-1.5672659458084428,0.9586752047146816,-2.1000823351718276,-1.9367670855405619,1.7663711706335692,2.6667816643058257,1.8724787726235608,1.8046221246730965,1.9976428507013653,1.5774629548550214,-2.1844031277397113,-2.37681286014546,1.8253248447171244,-1.5788810768815014,2.763614968711961,1.6681254090840125,-1.7888085025949607,2.3528244840292976,2.108851562509698,-2.123697969429917,2.256149605774371,1.683618624322591,1.7270267908891983,-1.9919179603858945,2.5292524650525894,-2.0123543155284263,2.2964102006786278,2.3238588922438432,-1.1039576402783546,-1.6468348544564142,2.079125306570687,2.5986107739265316,3.136814509289019,1.5415636566306685,2.7792702256080624,-1.2772677494202722,1.8851749570909602,2.195995012457669,2.5626678045398177,2.191554964023757,1.4811772660466813,2.344792664548995,-1.6613090531393202,-1.2668545152633315,1.7232909690875893,1.9383809994394323,0.8188351209426508,-2.0759834897304574,2.252866167883927,-2.156333825892606,-1.7088827259272346,-1.9786464505743182,-2.0733588217616914,-1.9619698757998387,2.5760746251915223,1.791771719172268,1.4835843181072768,1.357101887367243,2.6201130358305034,1.240096643255935,-2.3583108261315373]
+