{-# LANGUAGE DeriveDataTypeable #-}
{-# LANGUAGE DeriveGeneric #-}
module ELynx.Tree.Distribution.BirthDeath
( BirthDeathDistribution (..),
cumulative,
density,
quantile,
)
where
import Data.Data
( Data,
Typeable,
)
import ELynx.Tree.Distribution.Types
import GHC.Generics (Generic)
import qualified Statistics.Distribution as D
data BirthDeathDistribution = BDD
{
bddTOr :: Time,
bddLa :: Rate,
bddMu :: Rate
}
deriving (Eq, Typeable, Data, Generic)
instance D.Distribution BirthDeathDistribution where
cumulative = cumulative
cumulative :: BirthDeathDistribution -> Time -> Double
cumulative (BDD t l m) x
| x <= 0 = 0
| x > t = 1
| otherwise = t1 * t2
where
d = l - m
t1 = (1.0 - exp (- d * x)) / (l - m * exp (- d * x))
t2 = (l - m * exp (- d * t)) / (1.0 - exp (- d * t))
instance D.ContDistr BirthDeathDistribution where
density = density
quantile = quantile
density :: BirthDeathDistribution -> Time -> Double
density (BDD t l m) x
| x < 0 = 0
| x > t = 0
| otherwise = d ** 2 * t1 * t2
where
d = l - m
t1 = exp (- d * x) / ((l - m * exp (- d * x)) ** 2)
t2 = (l - m * exp (- d * t)) / (1.0 - exp (- d * t))
quantile :: BirthDeathDistribution -> Double -> Time
quantile (BDD t l m) p
| p >= 0 && p <= 1 =
res
| otherwise =
error $
"PointProcess.quantile: p must be in range [0,1] but got "
++ show p
++ "."
where
d = l - m
t2 = (l - m * exp (- d * t)) / (1.0 - exp (- d * t))
res = (-1.0 / d) * log ((1.0 - p * l / t2) / (1.0 - p * m / t2))
instance D.ContGen BirthDeathDistribution where
genContVar = D.genContinuous