{-# LANGUAGE BangPatterns #-}
module Mcmc.Prior
( Prior,
PriorFunction,
PriorFunctionG,
noPrior,
greaterThan,
positive,
lessThan,
negative,
exponential,
gamma,
gammaMeanVariance,
gammaMeanOne,
gammaShapeScaleToMeanVariance,
gammaMeanVarianceToShapeScale,
logNormal,
normal,
uniform,
poisson,
product',
)
where
import Control.Monad
import Data.Maybe (fromMaybe)
import Data.Typeable
import Mcmc.Internal.SpecFunctions
import Mcmc.Statistics.Types
import Numeric.Log
type Prior = Log Double
type PriorFunction a = a -> Log Double
type PriorFunctionG a b = a -> Log b
noPrior :: RealFloat b => PriorFunctionG a b
noPrior :: forall b a. RealFloat b => PriorFunctionG a b
noPrior = forall a b. a -> b -> a
const Log b
1.0
{-# SPECIALIZE noPrior :: PriorFunction Double #-}
greaterThan :: RealFloat a => LowerBoundary a -> PriorFunctionG a a
greaterThan :: forall a. RealFloat a => a -> PriorFunctionG a a
greaterThan a
a a
x
| a
x forall a. Ord a => a -> a -> Bool
> a
a = Log a
1.0
| Bool
otherwise = Log a
0.0
{-# SPECIALIZE greaterThan :: Double -> PriorFunction Double #-}
positive :: RealFloat a => PriorFunctionG a a
positive :: forall a. RealFloat a => PriorFunctionG a a
positive = forall a. RealFloat a => a -> PriorFunctionG a a
greaterThan a
0
{-# SPECIALIZE positive :: PriorFunction Double #-}
lessThan :: RealFloat a => UpperBoundary a -> PriorFunctionG a a
lessThan :: forall a. RealFloat a => a -> PriorFunctionG a a
lessThan a
a a
x
| a
x forall a. Ord a => a -> a -> Bool
< a
a = Log a
1.0
| Bool
otherwise = Log a
0.0
{-# SPECIALIZE lessThan :: Double -> PriorFunction Double #-}
negative :: RealFloat a => PriorFunctionG a a
negative :: forall a. RealFloat a => PriorFunctionG a a
negative = forall a. RealFloat a => a -> PriorFunctionG a a
lessThan a
0.0
{-# SPECIALIZE negative :: PriorFunction Double #-}
exponential :: RealFloat a => Rate a -> PriorFunctionG a a
exponential :: forall a. RealFloat a => a -> PriorFunctionG a a
exponential a
l a
x
| a
l forall a. Ord a => a -> a -> Bool
<= a
0.0 = forall a. HasCallStack => [Char] -> a
error [Char]
"exponential: Rate is zero or negative."
| a
x forall a. Ord a => a -> a -> Bool
<= a
0.0 = Log a
0.0
| Bool
otherwise = Log a
ll forall a. Num a => a -> a -> a
* forall a. a -> Log a
Exp (forall a. Num a => a -> a
negate a
l forall a. Num a => a -> a -> a
* a
x)
where
ll :: Log a
ll = forall a. a -> Log a
Exp forall a b. (a -> b) -> a -> b
$ forall a. Floating a => a -> a
log a
l
{-# SPECIALIZE exponential :: Double -> PriorFunction Double #-}
gamma :: (Typeable a, RealFloat a) => Shape a -> Scale a -> PriorFunctionG a a
gamma :: forall a. (Typeable a, RealFloat a) => a -> a -> PriorFunctionG a a
gamma a
k a
t a
x
| a
k forall a. Ord a => a -> a -> Bool
<= a
0.0 = forall a. HasCallStack => [Char] -> a
error [Char]
"gamma: Shape is zero or negative."
| a
t forall a. Ord a => a -> a -> Bool
<= a
0.0 = forall a. HasCallStack => [Char] -> a
error [Char]
"gamma: Scale is zero or negative."
| a
x forall a. Ord a => a -> a -> Bool
<= a
0.0 = Log a
0.0
| Bool
otherwise = forall a. a -> Log a
Exp forall a b. (a -> b) -> a -> b
$ forall a. Floating a => a -> a
log a
x forall a. Num a => a -> a -> a
* (a
k forall a. Num a => a -> a -> a
- a
1.0) forall a. Num a => a -> a -> a
- (a
x forall a. Fractional a => a -> a -> a
/ a
t) forall a. Num a => a -> a -> a
- forall a. (Typeable a, RealFloat a) => a -> a
logGammaG a
k forall a. Num a => a -> a -> a
- forall a. Floating a => a -> a
log a
t forall a. Num a => a -> a -> a
* a
k
{-# SPECIALIZE gamma :: Double -> Double -> PriorFunction Double #-}
gammaMeanVariance :: (Typeable a, RealFloat a) => Mean a -> Variance a -> PriorFunctionG a a
gammaMeanVariance :: forall a. (Typeable a, RealFloat a) => a -> a -> PriorFunctionG a a
gammaMeanVariance a
m a
v = forall a. (Typeable a, RealFloat a) => a -> a -> PriorFunctionG a a
gamma a
k a
t
where
(a
k, a
t) = forall a. Fractional a => a -> a -> (a, a)
gammaMeanVarianceToShapeScale a
m a
v
{-# SPECIALIZE gammaMeanVariance :: Double -> Double -> PriorFunction Double #-}
gammaMeanOne :: (Typeable a, RealFloat a) => Shape a -> PriorFunctionG a a
gammaMeanOne :: forall a. (Typeable a, RealFloat a) => a -> PriorFunctionG a a
gammaMeanOne a
k = forall a. (Typeable a, RealFloat a) => a -> a -> PriorFunctionG a a
gamma a
k (forall a. Fractional a => a -> a
recip a
k)
{-# SPECIALIZE gammaMeanOne :: Double -> PriorFunction Double #-}
gammaShapeScaleToMeanVariance :: Num a => Shape a -> Scale a -> (Mean a, Variance a)
gammaShapeScaleToMeanVariance :: forall a. Num a => a -> a -> (a, a)
gammaShapeScaleToMeanVariance a
k a
t = let m :: a
m = a
k forall a. Num a => a -> a -> a
* a
t in (a
m, a
m forall a. Num a => a -> a -> a
* a
t)
{-# SPECIALIZE gammaShapeScaleToMeanVariance :: Double -> Double -> (Double, Double) #-}
gammaMeanVarianceToShapeScale :: Fractional a => Mean a -> Variance a -> (Shape a, Scale a)
gammaMeanVarianceToShapeScale :: forall a. Fractional a => a -> a -> (a, a)
gammaMeanVarianceToShapeScale a
m a
v = (a
m forall a. Num a => a -> a -> a
* a
m forall a. Fractional a => a -> a -> a
/ a
v, a
v forall a. Fractional a => a -> a -> a
/ a
m)
{-# SPECIALIZE gammaMeanVarianceToShapeScale :: Double -> Double -> (Double, Double) #-}
mLnSqrt2Pi :: RealFloat a => a
mLnSqrt2Pi :: forall a. RealFloat a => a
mLnSqrt2Pi = a
0.9189385332046727417803297364056176398613974736377834128171
{-# INLINE mLnSqrt2Pi #-}
logNormal :: RealFloat a => Mean a -> StandardDeviation a -> PriorFunctionG a a
logNormal :: forall a. RealFloat a => a -> a -> PriorFunctionG a a
logNormal a
m a
s a
x
| a
s forall a. Ord a => a -> a -> Bool
<= a
0.0 = forall a. HasCallStack => [Char] -> a
error [Char]
"logNormal: Standard deviation is zero or negative."
| a
x forall a. Ord a => a -> a -> Bool
<= a
0.0 = Log a
0.0
| Bool
otherwise = forall a. a -> Log a
Exp forall a b. (a -> b) -> a -> b
$ a
t forall a. Num a => a -> a -> a
+ a
e
where
t :: a
t = forall a. Num a => a -> a
negate forall a b. (a -> b) -> a -> b
$ forall a. RealFloat a => a
mLnSqrt2Pi forall a. Num a => a -> a -> a
+ forall a. Floating a => a -> a
log (a
x forall a. Num a => a -> a -> a
* a
s)
a :: a
a = forall a. Fractional a => a -> a
recip forall a b. (a -> b) -> a -> b
$ a
2.0 forall a. Num a => a -> a -> a
* a
s forall a. Num a => a -> a -> a
* a
s
b :: a
b = forall a. Floating a => a -> a
log a
x forall a. Num a => a -> a -> a
- a
m
e :: a
e = forall a. Num a => a -> a
negate forall a b. (a -> b) -> a -> b
$ a
a forall a. Num a => a -> a -> a
* a
b forall a. Num a => a -> a -> a
* a
b
normal :: RealFloat a => Mean a -> StandardDeviation a -> PriorFunctionG a a
normal :: forall a. RealFloat a => a -> a -> PriorFunctionG a a
normal a
m a
s a
x
| a
s forall a. Ord a => a -> a -> Bool
<= a
0 = forall a. HasCallStack => [Char] -> a
error [Char]
"normal: Standard deviation is zero or negative."
| Bool
otherwise = forall a. a -> Log a
Exp forall a b. (a -> b) -> a -> b
$ (-a
xm forall a. Num a => a -> a -> a
* a
xm forall a. Fractional a => a -> a -> a
/ (a
2 forall a. Num a => a -> a -> a
* a
s forall a. Num a => a -> a -> a
* a
s)) forall a. Num a => a -> a -> a
- a
denom
where
xm :: a
xm = a
x forall a. Num a => a -> a -> a
- a
m
denom :: a
denom = forall a. RealFloat a => a
mLnSqrt2Pi forall a. Num a => a -> a -> a
+ forall a. Floating a => a -> a
log a
s
{-# SPECIALIZE normal :: Double -> Double -> PriorFunction Double #-}
uniform :: RealFloat a => LowerBoundary a -> UpperBoundary a -> PriorFunctionG a a
uniform :: forall a. RealFloat a => a -> a -> PriorFunctionG a a
uniform a
a a
b a
x
| a
a forall a. Ord a => a -> a -> Bool
> a
b = forall a. HasCallStack => [Char] -> a
error [Char]
"uniform: Lower boundary is greater than upper boundary."
| a
x forall a. Ord a => a -> a -> Bool
< a
a = Log a
0.0
| a
x forall a. Ord a => a -> a -> Bool
> a
b = Log a
0.0
| Bool
otherwise = Log a
1.0
{-# SPECIALIZE uniform :: Double -> Double -> PriorFunction Double #-}
poisson :: (RealFloat a, Typeable a) => Rate a -> PriorFunctionG Int a
poisson :: forall a. (RealFloat a, Typeable a) => a -> PriorFunctionG Int a
poisson a
l Int
n
| a
l forall a. Ord a => a -> a -> Bool
< a
0.0 = forall a. HasCallStack => [Char] -> a
error [Char]
"poisson: Rate is zero or negative."
| Int
n forall a. Ord a => a -> a -> Bool
< Int
0 = Log a
0.0
| Bool
otherwise = forall a. a -> Log a
Exp forall a b. (a -> b) -> a -> b
$ forall a. Floating a => a -> a
log a
l forall a. Num a => a -> a -> a
* forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n forall a. Num a => a -> a -> a
- forall a b. (Integral a, RealFloat b, Typeable b) => a -> b
logFactorialG Int
n forall a. Num a => a -> a -> a
- a
l
product' :: RealFloat a => [Log a] -> Log a
product' :: forall a. RealFloat a => [Log a] -> Log a
product' = forall a. a -> Maybe a -> a
fromMaybe Log a
0 forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a. RealFloat a => [Log a] -> Maybe (Log a)
prodM
{-# SPECIALIZE product' :: [Log Double] -> Log Double #-}
prodM :: RealFloat a => [Log a] -> Maybe (Log a)
prodM :: forall a. RealFloat a => [Log a] -> Maybe (Log a)
prodM = forall (t :: * -> *) (m :: * -> *) b a.
(Foldable t, Monad m) =>
(b -> a -> m b) -> b -> t a -> m b
foldM (\ !Log a
acc Log a
x -> (Log a
acc forall a. Num a => a -> a -> a
* Log a
x) forall (f :: * -> *) a b. Functor f => a -> f b -> f a
<$ forall (f :: * -> *). Alternative f => Bool -> f ()
guard (Log a
acc forall a. Eq a => a -> a -> Bool
/= Log a
0.0)) Log a
1.0
{-# SPECIALIZE prodM :: [Log Double] -> Maybe (Log Double) #-}