{-# LANGUAGE DeriveDataTypeable #-}
{-# LANGUAGE DeriveGeneric #-}
module ELynx.Distribution.BirthDeath
( BirthDeathDistribution(..)
, cumulative
, density
, quantile
) where
import Data.Data (Data, Typeable)
import GHC.Generics (Generic)
import qualified Statistics.Distribution as D
import ELynx.Distribution.Types
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 [0,1] range. 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