module Bayes.Continuous(
CNMonad(..)
, CV
, DN
, VariableName(..)
, BayesianVariable(..)
, Distri
, ContinuousNetwork(..)
, ContinuousSample(..)
, InstantiationValue(..)
, CVI
, uniform
, normal
, beta
, beta'
, exponential
, execCN
, runCN
, evalCN
, runSampling
, continuousMCMCSampler
, (=:)
, histogram
, samplingHistograms
) where
import Bayes.Sampling
import Bayes
import Data.Monoid
import Control.Applicative
import Bayes.PrivateTypes
import Bayes.Network(NetworkMonad(..),runGraph,execGraph,evalGraph)
import Control.Monad.State.Strict
import Data.Maybe(fromJust)
import Control.Monad
import Control.Monad.Reader
import Bayes.Sampling(runSampling, samplingHistograms,gibbsMCMCSampler, ContinuousNetwork(..), ContinuousSample(..))
import Math.Gamma(gamma)
type CNMonad a = GraphMonad DirectedSG () Distri a
instance Monoid (DistributionSupport (Double,Double)) where
mempty = BoundedSupport (0.0,0.0)
(BoundedSupport (xa,ya)) `mappend` (BoundedSupport (xb,yb)) = BoundedSupport (min xa xb, max ya yb)
(BoundedSupport _) `mappend` (Unbounded a) = Unbounded a
(Unbounded a) `mappend` (BoundedSupport _) = Unbounded a
(Unbounded a) `mappend` (Unbounded b) = Unbounded (a+b)
newtype RN = RN (Reader (Sample DirectedSG CVI) Double)
instance Num RN where
(RN a) + (RN b) = RN $ liftA2 (+) a b
(RN a) (RN b) = RN $ liftA2 () a b
(RN a) * (RN b) = RN $ liftA2 (*) a b
abs (RN a) = RN $ liftA abs a
signum (RN a) = RN $ liftA signum a
fromInteger a = RN (return (fromInteger a))
instance Fractional RN where
(RN a) / (RN b) = RN $ liftA2 (/) a b
fromRational a = RN (return (fromRational a))
data DN = DN !([CV]) !RN
instance BayesianVariable DN where
vertex (DN [] _) = error "No vertex for this expression"
vertex (DN (v:[]) _) = vertex v
vertex (DN l _) = error "Too many vertices for this expression"
dependencies :: DN -> [CV]
dependencies (DN l _) = l
instance Instantiable DN Double CVI where
(=:) (DN [] _) x = error "No variable is related to this expression"
(=:) (DN (a:[]) _) x = a =: x
(=:) (DN l _) x = error "Too many variables are related to this expression"
value :: DN -> Reader (Sample DirectedSG CVI) Double
value (DN _ (RN v)) = v
instance Num DN where
(DN da na) + (DN db nb) = DN (da `mappend`db) (na + nb)
(DN da na) (DN db nb) = DN (da `mappend`db) (na nb)
(DN da na) * (DN db nb) = DN (da `mappend`db) (na * nb)
abs (DN a b) = DN a (abs b)
signum (DN a b) = DN a (signum b)
fromInteger i = DN [] (fromInteger i)
instance Fractional DN where
(DN da na) / (DN db nb) = DN (da `mappend`db) (na / nb)
fromRational i = DN [] (fromRational i)
class VariableName m where
mkVariable :: m -> CNMonad CV
instance VariableName String where
mkVariable = variable
instance VariableName () where
mkVariable _ = unamedVariable
addVariableIfNotFound :: String -> CNMonad Vertex
addVariableIfNotFound vertexName = graphNode vertexName (D (CV (Vertex 0)) (Distri (return mempty) (\s -> return 0.0)))
unamedVariable :: CNMonad CV
unamedVariable = do
((namemap,count),g) <- get
let vertexName = "unamed" ++ show count
va <- addVariableIfNotFound vertexName
return (CV va)
whenJust Nothing _ = return ()
whenJust (Just i) f = f i >> return ()
setBayesianNode :: CV -> Distri -> CNMonad ()
setBayesianNode (CV v) newValue = do
(aux,oldGraph) <- get
let newGraph = changeVertexValue v newValue oldGraph
whenJust newGraph $ \nvm -> do
put $! (aux, nvm)
(<--) :: CV -> CV -> CNMonad ()
(CV va) <-- (CV vb) = newEdge vb va mempty
node :: CV -> DN
node v = DN [v] (RN $ do
r <- ask
return . instantiationValue . fromJust . vertexValue r $ (vertex v)
)
variable :: String
-> CNMonad CV
variable name = do
va <- addVariableIfNotFound name
return (CV va)
cpt :: CV -> [CV] -> CNMonad CV
cpt node conditions = do
mapM_ (node <--) (reverse conditions)
return node
exponential :: VariableName s
=> s
-> DN
-> CNMonad DN
exponential s na = do
v <- mkVariable s
let la = dependencies na
cpt v la
let bound = do
lambda <- value na
return (Unbounded (1.0 / lambda))
result = D v (Distri bound $ \inst -> do
lambda <- value na
let x = instantiationValue inst
if x < 0
then
return 0.0
else
return $ lambda * exp(lambda * x)
)
setBayesianNode v result
return (node v)
beta :: VariableName s
=> s
-> DN
-> DN
-> CNMonad DN
beta s na nb = do
v <- mkVariable s
let la = dependencies na
lb = dependencies nb
let l = la ++ lb
cpt v l
let bound = return (BoundedSupport (0.0,0.1) )
result = D v (Distri bound $ \inst -> do
a <- value na
b <- value nb
let x = instantiationValue inst
return $ gamma(a+b)/gamma(a)/gamma(b)*x**(a1)*(1x)**(b1)
)
setBayesianNode v result
return (node v)
beta' :: VariableName s
=> s
-> DN
-> DN
-> CNMonad DN
beta' s na nb = do
v <- mkVariable s
let la = dependencies na
lb = dependencies nb
let l = la ++ lb
cpt v l
let bound = do
a <- value na
b <- value nb
let std = a * (a + b 1)
return (Unbounded std)
result = D v (Distri bound $ \inst -> do
a <- value na
b <- value nb
let x = instantiationValue inst
if x <= 0
then
return 0.0
else
return $ gamma(a+b)/gamma(a)/gamma(b)*x**(a1)*(1+x)**(ab)
)
setBayesianNode v result
return (node v)
uniform :: VariableName s
=> s
-> DN
-> DN
-> CNMonad DN
uniform s na nb = do
v <- mkVariable s
let la = dependencies na
lb = dependencies nb
let l = la ++ lb
cpt v l
let bound = do
a <- value na
b <- value nb
return (BoundedSupport (a,b) )
let result = D v (Distri bound $ \inst -> do
a <- value na
b <- value nb
let d = 1.0 / (b a)
if instantiationValue inst < a || instantiationValue inst > b then return 0.0 else return d)
setBayesianNode v result
return (node v)
normal :: VariableName s
=> s
-> DN
-> DN
-> CNMonad DN
normal s mean std = do
v <- mkVariable s
let la = dependencies mean
lb = dependencies std
let l = la ++ lb
cpt v l
let bound = do
s <- value std
return (Unbounded s )
let result = D v (Distri bound $ \inst -> do
m <- value mean
s <- value std
let x = instantiationValue inst
return (1.0 / (sqrt (2 * pi)*s) * exp((xm)*(xm)/(2*s*s))
))
setBayesianNode v result
return (node v)
runCN :: CNMonad a -> (a,ContinuousNetwork)
runCN = runGraph
execCN :: CNMonad a -> ContinuousNetwork
execCN = execGraph
evalCN :: CNMonad a -> a
evalCN = evalGraph