module GaussianProcessDemo where
import LazyPPL
import Distr
import Prelude hiding ((<>))
import Data.List
import Data.Map (empty,lookup,insert,size,keys,elems,Map)
import Numeric.LinearAlgebra hiding (step,size)
import Numeric.LinearAlgebra.Data hiding (step,size)
import Control.Monad
import Data.IORef
import System.IO.Unsafe
import Debug.Trace
import Graphics.Matplotlib hiding (density)
We can define a random function gp :: (Double -> Double -> Double) -> Prob (Double -> Double)
which describes a Gaussian
process. This requires an infinite number of random choices, first
because it is defined for arbitrarily large x
, but also it
needs an infinite number of random choices in each finite interval. This
is all dealt with using laziness. We cannot plot the function precisely,
but when we plot it with a fixed resolution and viewport, it is fine
because only finitely many random choices are triggered.
The following defines a GP prior generated with the radial basis function, a.k.a. squared exponential covariance.
= var * exp (-(x1 - x2)^2 / (2 * l^2))
rbf var l x1 x2 =
plotGPPrior do
<- mh 1 $ sample $ gp $ rbf 0.25 0.5
fws let xys = map (\f -> map (\x -> (x,f x)) [0,0.1..6]) $ map fst $ take 100 $ fws
"images/gp-prior.svg" [] xys (-7) 7 0.2 plotCoords
dataset :: [(Double, Double)]
= [(0,0.6), (1, 0.7), (2,1.2), (3,3.2), (4,6.8), (5, 8.2), (6,8.4)] dataset
g
plus
a random start point a
. (Note that we are treating this
g
as a function like any other. And we could also have
built this model with the second-order regress
function
from the other regression examples.)
example :: Meas (Double -> Double)
= do g <- sample $ gp $ rbf 1 1
example <- sample $ normal 0 3
a <- sample $ normal 0 3
b let f x = a + b * x + (g x)
-> score $ normalPdf (f x) 0.3 y)
forM_ dataset (\(x,y) return f
=
plotGPRegression do
<- mh 0.2 example
fws let xys = map (\f -> map (\x -> (x,f x)) [0,0.1..6]) $ map fst $ take 100 $ every 1000 $ drop 10000 $ fws
"images/gp-reg.svg" dataset xys (-2) 10 0.1 plotCoords
step :: (Double -> Double -> Double) -> Data.Map.Map Double Double -> Double -> Double -> Double
=
step cov table x seed let sig11 = cov x x
= matrix (size table) [cov x b | b <-keys table]
sig12 = matrix 1 [cov a x | a <-keys table]
sig21 = matrix (size table) [cov a b | a <- keys table , b <-keys table]
sig22 = sig12 <> (pinvTol 1E8 sig22)
regCoeff = (regCoeff <> (matrix 1 $ elems table)) `atIndex` (0,0)
mu = sig11 - ((regCoeff <> sig21) `atIndex` (0,0)) in
var + seed * (sqrt $ if var > -0.01 then (abs var) else trace (show var) $ abs var)
mu
gp :: (Double -> Double -> Double) -> Prob (Double -> Double)
= do ns <- iid $ normal 0 1
gp cov return $ unsafePerformIO $ do
<- newIORef Data.Map.empty
ref 0 0)
modifyIORef' ref (Data.Map.insert return $ \x -> unsafePerformIO $ do
<- readIORef ref
table case Data.Map.lookup x table of
Just y -> do {return y}
Nothing -> do let y = step cov table x $ ns !! (1 + size table)
modifyIORef' ref (Data.Map.insert x y)return y
plotCoords :: String -> [(Double,Double)] -> [[(Double,Double)]] -> Double -> Double -> Double -> IO ()
=
plotCoords filename dataset xyss ymin ymax alpha do putStrLn $ "Plotting " ++ filename ++ "..."
$ foldl (\a xys -> a % plot (map fst xys) (map snd xys) @@ [o1 "go-", o2 "linewidth" (0.5 :: Double), o2 "alpha" alpha, o2 "ms" (0 :: Int)]) (scatter (map fst dataset) (map snd dataset) @@ [o2 "c" "black"] % xlim (0 :: Int) (6 :: Int) % ylim ymin ymax) xyss
file filename putStrLn "Done."
return ()
main :: IO ()
= do { plotGPPrior ; plotGPRegression } main