module WienerDemo where
import LazyPPL
import LazyPPL.Distributions
import LazyPPL.Distributions.GP (wiener)
import Data.List
import Data.Map (empty,lookup,insert,size,keys)
import Control.Monad
import Data.IORef
import System.IO.Unsafe
import Graphics.Matplotlib hiding (density)
wiener :: Prob(Double -> Double)
which describes a Wiener 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, because it is no-where differentiable. This is
all dealt with using laziness. We can’t plot this unbounded,
undifferentiable function precisely, but when we plot it with a fixed
resolution and viewport, the necessary finite random choices are
triggered.
=
plotWienerPrior do
<- mh 1 $ sample wiener
fws let xys = map (\f -> map (\x -> (x,f x)) [0,0.1..6]) $ map fst $ take 100 $ fws
"images/wiener-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 wiener
example <- sample $ normal 0 3
a let f x = a + 2 * (g x)
-> score $ normalPdf (f x) 0.3 y)
forM_ dataset (\(x,y) return f
=
plotWienerRegression do
<- mh 0.1 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/wiener-reg.svg" dataset xys (-2) 10 0.1 plotCoords
Recall the splice function
splice :: Prob [Double] -> Prob (Double -> Double) -> Prob (Double -> Double)
from the linear
regression example, which pieces together draws from a random
function using a point process. We can immediately apply this to the
Wiener process, to get a jump diffusion process.
splice :: Prob [Double] -> Prob (Double -> Double) -> Prob (Double -> Double)
=
splice pointProcess randomFun do
<- pointProcess
xs <- mapM (const randomFun) xs
fs <- randomFun
default_f let h :: [(Double, Double -> Double)] -> Double -> Double
= default_f x
h [] x : xfs) x | x <= a = f x
h ((a, f) : xfs) x | x > a = h xfs x
h ((a, f) return (h (zip xs fs))
poissonPP :: Double -> Double -> Prob [Double]
=
poissonPP lower rate do
<- exponential rate
step let x = lower + step
<- poissonPP x rate
xs return (x : xs)
regress :: Double -> Prob (a -> Double) -> [(a, Double)] -> Meas (a -> Double)
=
regress sigma prior dataset do
<- sample prior
f -> score $ normalPdf (f x) sigma y)
forM_ dataset (\(x, y) return f
jump :: Prob (Double -> Double)
= let p = do f <- wiener
jump <- normal 0 3
a return $ \x -> a + f x
in splice (poissonPP 0 0.2) p
Here are six samples from this distribution.
=
plotJumpPrior do
<- mh 1 $ sample jump
fws let xys = map (\f -> map (\x -> (x,f x)) [0,0.02..6]) $ map fst $ take 6 $ fws
"images/wiener-jump-prior.svg" [] xys (-7) 7 1 plotCoords
datasetB :: [(Double, Double)]
= [(0,0.6), (1, 0.7), (2,8.2), (3,9.1), (4,3.2), (5,4.9), (6,2.9)]
datasetB
=
plotJumpRegression do
<- mhirreducible 0.2 0.1 (regress 0.3 jump datasetB)
fws let xys = map (\f -> map (\x -> (x,f x)) [0,0.02..6]) $ map fst $ take 100 $ every 1000 $ drop 300000 $ fws
"images/wiener-jump-reg.svg" datasetB xys (-2) 10 0.1 plotCoords
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 { plotWienerPrior ; plotWienerRegression ; plotJumpPrior ; plotJumpRegression } main