{-# LANGUAGE DeriveDataTypeable #-}
{-# LANGUAGE DeriveGeneric #-}
module ELynx.Distribution.BirthDeathNearlyCritical
( BirthDeathNearlyCriticalDistribution(..)
, cumulative
, density
, quantile
) where
import Data.Data (Data, Typeable)
import GHC.Generics (Generic)
import qualified Statistics.Distribution as D
import ELynx.Distribution.Types
data BirthDeathNearlyCriticalDistribution = BDNCD
{ bdncdTOr :: Time
, bdncdLa :: Rate
, bdncdMu :: Rate
} deriving (Eq, Typeable, Data, Generic)
instance D.Distribution BirthDeathNearlyCriticalDistribution where
cumulative = cumulative
cumulative :: BirthDeathNearlyCriticalDistribution -> Time -> Double
cumulative (BDNCD t l m) s
| s <= 0 = 0
| s > t = 1
| otherwise = o0 + o1
where o0 = s * (1.0 + t*l) / t / (1.0 + s*l)
o1 = (-s*s + s*t) * (m - l) / (2.0*t * (1.0 + s*l)**2)
instance D.ContDistr BirthDeathNearlyCriticalDistribution where
density = density
quantile = quantile
density :: BirthDeathNearlyCriticalDistribution -> Time -> Double
density (BDNCD t l m) s
| s < 0 = 0
| s > t = 0
| otherwise = o0 + o1
where
o0 = (1.0 + t*l) / (t * (1.0 + s*l)**2)
o1 = (-2.0*s + t - s*t*l) * (m - l) / (2.0*t * (1.0 + s*l)**3)
quantile :: BirthDeathNearlyCriticalDistribution -> Double -> Time
quantile (BDNCD t l m) p
| p >= 0 && p <= 1 = res
| otherwise =
error $ "PointProcess.quantile: p must be in [0,1] range. Got: " ++ show p ++ "."
where
den = l*(-3.0 + 2.0*t*(-1.0+p)*l)+m
t1 = (2.0 + t*(l - 4.0*p*l + m)) / den
t2Nom = 4.0 + t*(l*(4.0 + t*l + 8.0*p*(1.0 + t*l)) + 2.0*(2.0 + t*l - 4.0*p*(1.0 + t*l))*m + t*m*m)
t2 = t2Nom / (den**2)
res = 0.5 * (t1 + sqrt t2)
instance D.ContGen BirthDeathNearlyCriticalDistribution where
genContVar = D.genContinuous