module Bayes.Sampling(
Sampler(..)
, Sample(..)
, runSampling
, topologicalOrder
, discreteAncestralSampler
, gibbsSampler
, gibbsMCMCSampler
, samplingHistograms
, histogram
, ContinuousNetwork(..)
, ContinuousSample(..)
, Distri(..)
, continuousMCMCSampler
) where
import Bayes
import Bayes.Factor
import Control.Monad.Reader
import Data.Maybe(fromJust)
import Bayes.PrivateTypes
import System.Random.MWC.CondensedTable
import System.Random.MWC.Distributions(normal)
import qualified Data.Vector as V
import qualified Data.Vector.Unboxed as UV
import qualified Data.Vector.Unboxed.Mutable as MV(write,read)
import System.Random.MWC(GenIO,create,withSystemRandom,asGenIO,uniformR,Variate)
type SamplerGraph g a = BayesianNetwork g (Sample g a -> IO a)
type SamplingScheme g b a = GenIO -> SamplerGraph g a -> b -> Sample g a -> IO (b,Sample g a)
data Sampler g a = forall b . Sampler !b !(GenIO -> IO (Sample g a)) !(GenIO -> SamplerGraph g a) !(SamplingScheme g b a)
data Distri = D !CV !(DistributionF DirectedSG (Double,Double) CVI)
type ContinuousNetwork = SBN Distri
type ContinuousSample = SBN CVI
distributionOnNetwork ::(DistributionF DirectedSG (Double,Double) CVI) -> Vertex -> Reader (Sample DirectedSG CVI) Double
distributionOnNetwork (Distri _ f) v = do
r <- ask
let inst = fromJust . vertexValue r $ v
f inst
instance Graph g => Show (Sample g [DVI]) where
show = show . allNodes
instance Graph g => Show (Sample g CVI) where
show = show . allNodes
instance Graph g => Show (Sample g [(Double,Double,Double)]) where
show = show . allNodes
topologicalOrder :: DirectedGraph g => g a b -> [Vertex]
topologicalOrder g = _topologicalOrder g []
where
_topologicalOrder g current = case rootNode g of
Nothing -> reverse current
Just r -> _topologicalOrder (removeVertex r g) (r:current)
generateNewSample :: Graph g => SamplerGraph g a -> Sample g a -> Vertex -> IO (Sample g a)
generateNewSample sampler sample theVertex = do
let vertexSampler = fromJust $ vertexValue sampler theVertex
newValue <- vertexSampler sample
return . fromJust . changeVertexValue theVertex newValue $ sample
ancestralSampling :: DirectedGraph g => SamplingScheme g [Vertex] a
ancestralSampling gen sampler nodes s = do
r <- foldM (generateNewSample sampler) s nodes
return (nodes,r)
gibbsSampling :: DirectedGraph g => SamplingScheme g (CondensedTable V.Vector Vertex) a
gibbsSampling gen samplerGraph table s = do
aVertex <- genFromTable table gen
r <- generateNewSample samplerGraph s aVertex
return (table,r)
sequenceG :: (FunctorWithVertex g, Graph g) => [g a b] -> g a [b]
sequenceG [] = error "Can't sequence an empty list of graphs"
sequenceG l@(h:_) =
let setVertexList vertex oldValue = map (\g -> fromJust . vertexValue g $ vertex) l
in
fmapWithVertex setVertexList h
_range :: Double -> Double -> [Double] -> (Double,Double)
_range mi ma [] = (mi,ma)
_range mi ma (h:t) | h < mi = _range h ma t
| h > ma = _range mi h t
| otherwise = _range mi ma t
range :: [Double] -> (Double,Double)
range = _range (1/0) (1/0)
histogram :: Int
-> [Double]
-> [(Double,Double,Double)]
histogram bins values | bins <= 0 = error "Bins must be >= 1"
| bins == 1 =
let (bmin,bmax) = range values
in
[(bmin,bmax,1.0)]
| otherwise =
let (bmin,bmax) = range values
in
if bmax == bmin
then
[(bmin,bmax,1.0)]
else
let delta = (bmax bmin) / fromIntegral bins
v = UV.replicate bins 0.0
orMAX a | a >= bins = bins 1
| otherwise = a
countSample value = UV.modify (\s -> do
let index = orMAX $ floor ((value bmin) / delta)
d <- MV.read s index
MV.write s index (d+1))
v' = foldr countSample v values
t = UV.sum v'
scaled = UV.map (/ t) v'
addBounds value i = (bmin + delta * i ,bmin + delta + delta*i,value)
in
zipWith addBounds (UV.toList scaled) [0..]
samplingHistograms :: (InstantiationValue i v,BayesianVariable i, FunctorWithVertex g, Graph g)
=> Int
-> [Sample g i]
-> Sample g [(Double,Double,Double)]
samplingHistograms bins g' =
let g = sequenceG g'
in
fmapWithVertex (createHistogram bins g) g
where
createHistogram bins g vertex l = histogram bins . map toDouble $ l
runSampling :: (DirectedGraph g, FunctorWithVertex g)
=> Int
-> Int
-> Sampler g a
-> IO [Sample g a]
runSampling n b (Sampler initState startSample sampler samplingScheme) = do
withSystemRandom . asGenIO $ \gen -> do
start <- startSample gen
l <- _runAncestralSampling gen (n+b) initState start []
return (drop b $ l)
where
_runAncestralSampling gen n istate i r | n <= 1 = return $ reverse (i:r)
| otherwise = do
(newState,i') <- samplingScheme gen (sampler gen) istate i
_runAncestralSampling gen (n1) newState i' (i':r)
gibbsSampler :: (Factor f, FunctorWithVertex g, DirectedGraph g)
=> BayesianNetwork g f
-> [DVI]
-> Sampler g DVI
gibbsSampler g evidence =
let isEvidence v = v `elem` (map vertex evidence)
v = filter (not . isEvidence) (allVertices g)
t = tableFromWeights $ V.fromList $ zip v (repeat 1.0)
in
Sampler t (gibbsInitValue evidence g) (gibbsGraphSampler g) gibbsSampling
gibbsMCMCSampler :: (Factor f, FunctorWithVertex g, DirectedGraph g)
=> BayesianNetwork g f
-> [DVI]
-> Sampler g DVI
gibbsMCMCSampler g evidence =
let isEvidence v = v `elem` (map vertex evidence)
v = filter (not . isEvidence) (allVertices g)
t = tableFromWeights $ V.fromList $ zip v (repeat 1.0)
in
Sampler t (gibbsInitValue evidence g) (gibbsMCMCGraphSampler g) gibbsSampling
continuousMCMCSampler :: ContinuousNetwork
-> [CVI]
-> Sampler DirectedSG CVI
continuousMCMCSampler g evidence =
let isEvidence v = v `elem` (map vertex evidence)
v = filter (not . isEvidence) (allVertices g)
t = tableFromWeights $ V.fromList $ zip v (repeat 1.0)
in
Sampler t (continuousInitValue evidence g) (gibbsContinuousMCMCGraphSampler g) gibbsSampling
discreteAncestralSampler :: (Factor f, FunctorWithVertex g, DirectedGraph g) => BayesianNetwork g f -> Sampler g DVI
discreteAncestralSampler g =
let nodes = topologicalOrder g
in
Sampler nodes (discreteInitValue [] g) (discreteSampler g) ancestralSampling
distributionFromFactor :: (DirectedGraph g,Factor f) => f -> DistributionF g DV DVI
distributionFromFactor f =
let vars@(main:_) = factorVariables f
in
Distri (return $ BoundedSupport main) $ \current -> do
sample <- ask
let currentV = vertex main
sample' = fromJust . changeVertexValue currentV current $ sample
return $ factorValue f (variableInstantiations sample' vars)
distributionFromMarkovBlanket :: (DirectedGraph g,Factor f) => BayesianNetwork g f -> f -> DistributionF g DV DVI
distributionFromMarkovBlanket g f =
let vars@(main:_) = factorVariables f
in
Distri (return $ BoundedSupport main) $ \current -> do
sample <- ask
let currentV = vertex main
sample' = fromJust . changeVertexValue currentV current $ sample
children = childrenNodes sample currentV
childrenF = map (fromJust . vertexValue g) children
values = map (\aFactor -> factorValue aFactor (variableInstantiations sample' (factorVariables aFactor))) (f:childrenF)
return $ product values
continuousDistributionFromMarkovBlanket :: ContinuousNetwork
-> Distri
-> DistributionF DirectedSG (Double,Double) CVI
continuousDistributionFromMarkovBlanket g (D cv d@(Distri b distri)) = Distri b $ \current -> do
sample <- ask
let currentV = vertex cv
children = childrenNodes sample currentV
dist (D _ d) = d
childrenF = map (dist . fromJust . vertexValue g) children
local (fromJust . changeVertexValue currentV current) $ do
s' <- ask
values <- mapM (uncurry distributionOnNetwork) ((d,currentV):zip childrenF children)
return $ product $ values
forAllRange :: Graph g => Sample g DVI -> DistributionF g DV DVI -> [Double]
forAllRange g (Distri bound r) =
let range = allInstantiationsForOneVariable theBound
in
map (flip runReader g . r) range
where
BoundedSupport theBound = runReader bound g
distributionValue :: Graph g => DistributionF g bounds inst -> Sample g inst -> inst -> Double
distributionValue (Distri bound r) g v = runReader (r v) g
distribution :: Graph g => Sample g DVI -> DistributionF g DV DVI -> CondensedTableV Int
distribution g f =
let notNullProba (i,p) = p /= 0.0
values = forAllRange g f
in
tableFromWeights $ V.fromList $ filter notNullProba $ zip [0..] values
variableInstantiations :: (DirectedGraph g, BayesianVariable v) => Sample g a -> [v] -> [a]
variableInstantiations g vars =
let vertices = map vertex vars
in
map (fromJust . vertexValue g) vertices
withEvidence :: (Monad m, BayesianVariable v, BayesianVariable vb) => [v] -> vb -> (vb -> m v) -> m v
withEvidence e v f =
if (vertex v) `elem` (map vertex e)
then
return . head . filter (\x -> vertex v == vertex x) $ e
else
f v
nullInstantiation :: Factor f => [DVI] -> Vertex -> f -> DVI
nullInstantiation evidence v f = setDVValue (factorMainVariable f) 0
randomInstantiation :: Factor f => GenIO -> [DVI] -> Vertex -> f -> IO DVI
randomInstantiation gen evidence vertex f = withEvidence evidence vertex $ \v -> do
let main = factorMainVariable f
t = tableFromWeights $ V.fromList $ zip [0..dimension main 1] (repeat 1.0)
r <- genFromTable t gen
return (setDVValue main r)
continuousRandomInstantiation :: GenIO -> [CVI] -> Vertex -> Distri -> IO CVI
continuousRandomInstantiation gen evidence vertex (D cv (Distri b _)) = withEvidence evidence vertex $ \v -> do
value <- uniformR (10.0,10.0) gen :: IO Double
return (cv =: value)
updateOneSample :: (DirectedGraph g,Factor f) => GenIO -> Vertex -> f -> (Sample g DVI -> IO DVI)
updateOneSample gen v f = \s -> do
let d = distributionFromFactor f
newVal <- genFromTable (distribution s d) gen
let v = factorMainVariable f
return (setDVValue v newVal)
gibbsUpdateOneSample :: (DirectedGraph g,Factor f) => BayesianNetwork g f -> GenIO -> Vertex -> f -> (Sample g DVI -> IO DVI)
gibbsUpdateOneSample g gen v f =
let d = distributionFromMarkovBlanket g f
v = factorMainVariable f
in
\s -> do
newVal <- genFromTable (distribution s d) gen
return (setDVValue v newVal)
class SamplingBounds bound inter | bound -> inter where
samplingBounds :: Num inter => bound -> (inter,inter)
class SamplingBounds bound inter => SampleGeneration var inst bound inter | var -> inst, var -> bound, var -> inter where
generateSample :: Instantiable var inter inst => Sample g inst -> DistributionF g bound inst -> GenIO -> var -> IO inst
instance SamplingBounds DV Int where
samplingBounds dv = (0, dimension dv 1)
instance SampleGeneration DV DVI DV Int where
generateSample s (Distri b _) gen var = do
let theBound = runReader b s
case theBound of
BoundedSupport s -> do
let (sa,sb) = samplingBounds s
u <- uniformR (sa,sb) gen
return $ setDVValue var u
Unbounded v -> error "DVI are bounded"
instance SamplingBounds (Double,Double) Double where
samplingBounds (sa,sb) = (sa,sb)
instance SampleGeneration CV CVI (Double,Double) Double where
generateSample s (Distri b _) gen var = do
let theBound = runReader b s
case theBound of
BoundedSupport s -> do
let (sa,sb) = samplingBounds s
u <- uniformR (sa,sb) gen
return $ var =: u
Unbounded sigma -> do
n <- normal 0.0 sigma gen
return $ var =: n
gibbsMCMCUpdateOneSample :: (DirectedGraph g,Factor f)
=> BayesianNetwork g f
-> GenIO
-> Vertex
-> f
-> (Sample g DVI -> IO DVI)
gibbsMCMCUpdateOneSample g gen v f =
let d = distributionFromMarkovBlanket g f
main = factorMainVariable f
v = vertex main
in
_gibbsMCMCUpdateOneSample d gen v main
gibbsContinuousMCMCUpdateOneSample :: ContinuousNetwork
-> GenIO
-> Vertex
-> Distri
-> (ContinuousSample -> IO CVI)
gibbsContinuousMCMCUpdateOneSample g gen v distri@(D cv _) =
let d = continuousDistributionFromMarkovBlanket g distri
v = vertex cv
in
_gibbsMCMCUpdateOneSample d gen v cv
_gibbsMCMCUpdateOneSample :: (DirectedGraph g,BayesianVariable var, Instantiable var inter inst, SampleGeneration var inst bound inter)
=> DistributionF g bound inst
-> GenIO
-> Vertex
-> var
-> (Sample g inst -> IO inst)
_gibbsMCMCUpdateOneSample d@(Distri _ _) gen v var = \s -> do
newSample <- generateSample s d gen var
let oldSample = fromJust . vertexValue s $ v
oldProbaValue = distributionValue d s oldSample
newProbaValue = distributionValue d s newSample
if oldProbaValue == 0.0
then
return newSample
else
let r = newProbaValue / oldProbaValue
in
if (r >= 1.0)
then
return newSample
else do
z <- uniformR (0.0, 1.0) gen
if z < r
then
return newSample
else
return oldSample
discreteSampler :: (FunctorWithVertex g, Factor f, DirectedGraph g) => BayesianNetwork g f -> GenIO -> SamplerGraph g DVI
discreteSampler g gen = fmapWithVertex (updateOneSample gen) g
gibbsGraphSampler :: (FunctorWithVertex g, Factor f, DirectedGraph g) => BayesianNetwork g f -> GenIO -> SamplerGraph g DVI
gibbsGraphSampler g gen = fmapWithVertex (gibbsUpdateOneSample g gen) g
gibbsMCMCGraphSampler :: (FunctorWithVertex g, Factor f, DirectedGraph g) => BayesianNetwork g f -> GenIO -> SamplerGraph g DVI
gibbsMCMCGraphSampler g gen = fmapWithVertex (gibbsMCMCUpdateOneSample g gen) g
gibbsContinuousMCMCGraphSampler :: ContinuousNetwork
-> GenIO
-> SamplerGraph DirectedSG CVI
gibbsContinuousMCMCGraphSampler g gen = fmapWithVertex (gibbsContinuousMCMCUpdateOneSample g gen) g
discreteInitValue :: (FunctorWithVertex g, Factor f, DirectedGraph g) => [DVI] -> BayesianNetwork g f -> GenIO -> IO (Sample g DVI)
discreteInitValue evidences g gen = fmapWithVertexM (randomInstantiation gen evidences) g
gibbsInitValue :: (FunctorWithVertex g, Factor f, DirectedGraph g) => [DVI] -> BayesianNetwork g f -> GenIO -> IO (Sample g DVI)
gibbsInitValue = discreteInitValue
continuousInitValue ::[CVI] -> ContinuousNetwork -> GenIO -> IO ContinuousSample
continuousInitValue evidences g gen = fmapWithVertexM (continuousRandomInstantiation gen evidences) g