{-# LANGUAGE DeriveDataTypeable #-}
{-# LANGUAGE DeriveGeneric #-}

-- |
--   Module      :  ELynx.Tree.Distribution.TimeOfOriginNearCritical
--   Description :  Distribution of time of origin for birth and death trees
--   Copyright   :  (c) Dominik Schrempf 2021
--   License     :  GPL-3.0-or-later
--
--   Maintainer  :  dominik.schrempf@gmail.com
--   Stability   :  unstable
--   Portability :  portable
--
-- Creation date: Tue Feb 13 13:16:18 2018.
--
-- See Gernhard, T. (2008). The conditioned reconstructed process. Journal of
-- Theoretical Biology, 253(4), 769–778. http://doi.org/10.1016/j.jtbi.2008.04.005.
--
-- Distribution of the time of origin for birth and death trees. See corollary 3.3
-- in the paper cited above.
module ELynx.Tree.Distribution.TimeOfOriginNearCritical
  ( TimeOfOriginNearCriticalDistribution (..),
    cumulative,
    density,
    quantile,
  )
where

import Data.Data
  ( Data,
    Typeable,
  )
import ELynx.Tree.Distribution.Types
import GHC.Generics (Generic)
import qualified Statistics.Distribution as D

-- | Distribution of the time of origin for a phylogenetic tree evolving under
-- the birth and death process and conditioned on observing n leaves today.
data TimeOfOriginNearCriticalDistribution = TONCD
  { -- | Number of leaves of the tree.
    TimeOfOriginNearCriticalDistribution -> Int
todTN :: Int,
    -- | Birth rate.
    TimeOfOriginNearCriticalDistribution -> Rate
todLa :: Rate,
    -- | Death rate.
    TimeOfOriginNearCriticalDistribution -> Rate
todMu :: Rate
  }
  deriving (TimeOfOriginNearCriticalDistribution
-> TimeOfOriginNearCriticalDistribution -> Bool
(TimeOfOriginNearCriticalDistribution
 -> TimeOfOriginNearCriticalDistribution -> Bool)
-> (TimeOfOriginNearCriticalDistribution
    -> TimeOfOriginNearCriticalDistribution -> Bool)
-> Eq TimeOfOriginNearCriticalDistribution
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
/= :: TimeOfOriginNearCriticalDistribution
-> TimeOfOriginNearCriticalDistribution -> Bool
$c/= :: TimeOfOriginNearCriticalDistribution
-> TimeOfOriginNearCriticalDistribution -> Bool
== :: TimeOfOriginNearCriticalDistribution
-> TimeOfOriginNearCriticalDistribution -> Bool
$c== :: TimeOfOriginNearCriticalDistribution
-> TimeOfOriginNearCriticalDistribution -> Bool
Eq, Typeable, Typeable TimeOfOriginNearCriticalDistribution
DataType
Constr
Typeable TimeOfOriginNearCriticalDistribution
-> (forall (c :: * -> *).
    (forall d b. Data d => c (d -> b) -> d -> c b)
    -> (forall g. g -> c g)
    -> TimeOfOriginNearCriticalDistribution
    -> c TimeOfOriginNearCriticalDistribution)
-> (forall (c :: * -> *).
    (forall b r. Data b => c (b -> r) -> c r)
    -> (forall r. r -> c r)
    -> Constr
    -> c TimeOfOriginNearCriticalDistribution)
-> (TimeOfOriginNearCriticalDistribution -> Constr)
-> (TimeOfOriginNearCriticalDistribution -> DataType)
-> (forall (t :: * -> *) (c :: * -> *).
    Typeable t =>
    (forall d. Data d => c (t d))
    -> Maybe (c TimeOfOriginNearCriticalDistribution))
-> (forall (t :: * -> * -> *) (c :: * -> *).
    Typeable t =>
    (forall d e. (Data d, Data e) => c (t d e))
    -> Maybe (c TimeOfOriginNearCriticalDistribution))
-> ((forall b. Data b => b -> b)
    -> TimeOfOriginNearCriticalDistribution
    -> TimeOfOriginNearCriticalDistribution)
-> (forall r r'.
    (r -> r' -> r)
    -> r
    -> (forall d. Data d => d -> r')
    -> TimeOfOriginNearCriticalDistribution
    -> r)
-> (forall r r'.
    (r' -> r -> r)
    -> r
    -> (forall d. Data d => d -> r')
    -> TimeOfOriginNearCriticalDistribution
    -> r)
-> (forall u.
    (forall d. Data d => d -> u)
    -> TimeOfOriginNearCriticalDistribution -> [u])
-> (forall u.
    Int
    -> (forall d. Data d => d -> u)
    -> TimeOfOriginNearCriticalDistribution
    -> u)
-> (forall (m :: * -> *).
    Monad m =>
    (forall d. Data d => d -> m d)
    -> TimeOfOriginNearCriticalDistribution
    -> m TimeOfOriginNearCriticalDistribution)
-> (forall (m :: * -> *).
    MonadPlus m =>
    (forall d. Data d => d -> m d)
    -> TimeOfOriginNearCriticalDistribution
    -> m TimeOfOriginNearCriticalDistribution)
-> (forall (m :: * -> *).
    MonadPlus m =>
    (forall d. Data d => d -> m d)
    -> TimeOfOriginNearCriticalDistribution
    -> m TimeOfOriginNearCriticalDistribution)
-> Data TimeOfOriginNearCriticalDistribution
TimeOfOriginNearCriticalDistribution -> DataType
TimeOfOriginNearCriticalDistribution -> Constr
(forall b. Data b => b -> b)
-> TimeOfOriginNearCriticalDistribution
-> TimeOfOriginNearCriticalDistribution
(forall d b. Data d => c (d -> b) -> d -> c b)
-> (forall g. g -> c g)
-> TimeOfOriginNearCriticalDistribution
-> c TimeOfOriginNearCriticalDistribution
(forall b r. Data b => c (b -> r) -> c r)
-> (forall r. r -> c r)
-> Constr
-> c TimeOfOriginNearCriticalDistribution
forall a.
Typeable a
-> (forall (c :: * -> *).
    (forall d b. Data d => c (d -> b) -> d -> c b)
    -> (forall g. g -> c g) -> a -> c a)
-> (forall (c :: * -> *).
    (forall b r. Data b => c (b -> r) -> c r)
    -> (forall r. r -> c r) -> Constr -> c a)
-> (a -> Constr)
-> (a -> DataType)
-> (forall (t :: * -> *) (c :: * -> *).
    Typeable t =>
    (forall d. Data d => c (t d)) -> Maybe (c a))
-> (forall (t :: * -> * -> *) (c :: * -> *).
    Typeable t =>
    (forall d e. (Data d, Data e) => c (t d e)) -> Maybe (c a))
-> ((forall b. Data b => b -> b) -> a -> a)
-> (forall r r'.
    (r -> r' -> r) -> r -> (forall d. Data d => d -> r') -> a -> r)
-> (forall r r'.
    (r' -> r -> r) -> r -> (forall d. Data d => d -> r') -> a -> r)
-> (forall u. (forall d. Data d => d -> u) -> a -> [u])
-> (forall u. Int -> (forall d. Data d => d -> u) -> a -> u)
-> (forall (m :: * -> *).
    Monad m =>
    (forall d. Data d => d -> m d) -> a -> m a)
-> (forall (m :: * -> *).
    MonadPlus m =>
    (forall d. Data d => d -> m d) -> a -> m a)
-> (forall (m :: * -> *).
    MonadPlus m =>
    (forall d. Data d => d -> m d) -> a -> m a)
-> Data a
forall u.
Int
-> (forall d. Data d => d -> u)
-> TimeOfOriginNearCriticalDistribution
-> u
forall u.
(forall d. Data d => d -> u)
-> TimeOfOriginNearCriticalDistribution -> [u]
forall r r'.
(r -> r' -> r)
-> r
-> (forall d. Data d => d -> r')
-> TimeOfOriginNearCriticalDistribution
-> r
forall r r'.
(r' -> r -> r)
-> r
-> (forall d. Data d => d -> r')
-> TimeOfOriginNearCriticalDistribution
-> r
forall (m :: * -> *).
Monad m =>
(forall d. Data d => d -> m d)
-> TimeOfOriginNearCriticalDistribution
-> m TimeOfOriginNearCriticalDistribution
forall (m :: * -> *).
MonadPlus m =>
(forall d. Data d => d -> m d)
-> TimeOfOriginNearCriticalDistribution
-> m TimeOfOriginNearCriticalDistribution
forall (c :: * -> *).
(forall b r. Data b => c (b -> r) -> c r)
-> (forall r. r -> c r)
-> Constr
-> c TimeOfOriginNearCriticalDistribution
forall (c :: * -> *).
(forall d b. Data d => c (d -> b) -> d -> c b)
-> (forall g. g -> c g)
-> TimeOfOriginNearCriticalDistribution
-> c TimeOfOriginNearCriticalDistribution
forall (t :: * -> *) (c :: * -> *).
Typeable t =>
(forall d. Data d => c (t d))
-> Maybe (c TimeOfOriginNearCriticalDistribution)
forall (t :: * -> * -> *) (c :: * -> *).
Typeable t =>
(forall d e. (Data d, Data e) => c (t d e))
-> Maybe (c TimeOfOriginNearCriticalDistribution)
$cTONCD :: Constr
$tTimeOfOriginNearCriticalDistribution :: DataType
gmapMo :: (forall d. Data d => d -> m d)
-> TimeOfOriginNearCriticalDistribution
-> m TimeOfOriginNearCriticalDistribution
$cgmapMo :: forall (m :: * -> *).
MonadPlus m =>
(forall d. Data d => d -> m d)
-> TimeOfOriginNearCriticalDistribution
-> m TimeOfOriginNearCriticalDistribution
gmapMp :: (forall d. Data d => d -> m d)
-> TimeOfOriginNearCriticalDistribution
-> m TimeOfOriginNearCriticalDistribution
$cgmapMp :: forall (m :: * -> *).
MonadPlus m =>
(forall d. Data d => d -> m d)
-> TimeOfOriginNearCriticalDistribution
-> m TimeOfOriginNearCriticalDistribution
gmapM :: (forall d. Data d => d -> m d)
-> TimeOfOriginNearCriticalDistribution
-> m TimeOfOriginNearCriticalDistribution
$cgmapM :: forall (m :: * -> *).
Monad m =>
(forall d. Data d => d -> m d)
-> TimeOfOriginNearCriticalDistribution
-> m TimeOfOriginNearCriticalDistribution
gmapQi :: Int
-> (forall d. Data d => d -> u)
-> TimeOfOriginNearCriticalDistribution
-> u
$cgmapQi :: forall u.
Int
-> (forall d. Data d => d -> u)
-> TimeOfOriginNearCriticalDistribution
-> u
gmapQ :: (forall d. Data d => d -> u)
-> TimeOfOriginNearCriticalDistribution -> [u]
$cgmapQ :: forall u.
(forall d. Data d => d -> u)
-> TimeOfOriginNearCriticalDistribution -> [u]
gmapQr :: (r' -> r -> r)
-> r
-> (forall d. Data d => d -> r')
-> TimeOfOriginNearCriticalDistribution
-> r
$cgmapQr :: forall r r'.
(r' -> r -> r)
-> r
-> (forall d. Data d => d -> r')
-> TimeOfOriginNearCriticalDistribution
-> r
gmapQl :: (r -> r' -> r)
-> r
-> (forall d. Data d => d -> r')
-> TimeOfOriginNearCriticalDistribution
-> r
$cgmapQl :: forall r r'.
(r -> r' -> r)
-> r
-> (forall d. Data d => d -> r')
-> TimeOfOriginNearCriticalDistribution
-> r
gmapT :: (forall b. Data b => b -> b)
-> TimeOfOriginNearCriticalDistribution
-> TimeOfOriginNearCriticalDistribution
$cgmapT :: (forall b. Data b => b -> b)
-> TimeOfOriginNearCriticalDistribution
-> TimeOfOriginNearCriticalDistribution
dataCast2 :: (forall d e. (Data d, Data e) => c (t d e))
-> Maybe (c TimeOfOriginNearCriticalDistribution)
$cdataCast2 :: forall (t :: * -> * -> *) (c :: * -> *).
Typeable t =>
(forall d e. (Data d, Data e) => c (t d e))
-> Maybe (c TimeOfOriginNearCriticalDistribution)
dataCast1 :: (forall d. Data d => c (t d))
-> Maybe (c TimeOfOriginNearCriticalDistribution)
$cdataCast1 :: forall (t :: * -> *) (c :: * -> *).
Typeable t =>
(forall d. Data d => c (t d))
-> Maybe (c TimeOfOriginNearCriticalDistribution)
dataTypeOf :: TimeOfOriginNearCriticalDistribution -> DataType
$cdataTypeOf :: TimeOfOriginNearCriticalDistribution -> DataType
toConstr :: TimeOfOriginNearCriticalDistribution -> Constr
$ctoConstr :: TimeOfOriginNearCriticalDistribution -> Constr
gunfold :: (forall b r. Data b => c (b -> r) -> c r)
-> (forall r. r -> c r)
-> Constr
-> c TimeOfOriginNearCriticalDistribution
$cgunfold :: forall (c :: * -> *).
(forall b r. Data b => c (b -> r) -> c r)
-> (forall r. r -> c r)
-> Constr
-> c TimeOfOriginNearCriticalDistribution
gfoldl :: (forall d b. Data d => c (d -> b) -> d -> c b)
-> (forall g. g -> c g)
-> TimeOfOriginNearCriticalDistribution
-> c TimeOfOriginNearCriticalDistribution
$cgfoldl :: forall (c :: * -> *).
(forall d b. Data d => c (d -> b) -> d -> c b)
-> (forall g. g -> c g)
-> TimeOfOriginNearCriticalDistribution
-> c TimeOfOriginNearCriticalDistribution
$cp1Data :: Typeable TimeOfOriginNearCriticalDistribution
Data, (forall x.
 TimeOfOriginNearCriticalDistribution
 -> Rep TimeOfOriginNearCriticalDistribution x)
-> (forall x.
    Rep TimeOfOriginNearCriticalDistribution x
    -> TimeOfOriginNearCriticalDistribution)
-> Generic TimeOfOriginNearCriticalDistribution
forall x.
Rep TimeOfOriginNearCriticalDistribution x
-> TimeOfOriginNearCriticalDistribution
forall x.
TimeOfOriginNearCriticalDistribution
-> Rep TimeOfOriginNearCriticalDistribution x
forall a.
(forall x. a -> Rep a x) -> (forall x. Rep a x -> a) -> Generic a
$cto :: forall x.
Rep TimeOfOriginNearCriticalDistribution x
-> TimeOfOriginNearCriticalDistribution
$cfrom :: forall x.
TimeOfOriginNearCriticalDistribution
-> Rep TimeOfOriginNearCriticalDistribution x
Generic)

instance D.Distribution TimeOfOriginNearCriticalDistribution where
  cumulative :: TimeOfOriginNearCriticalDistribution -> Rate -> Rate
cumulative = TimeOfOriginNearCriticalDistribution -> Rate -> Rate
cumulative

-- | Cumulative distribution function; see Mathematica notebook.
cumulative :: TimeOfOriginNearCriticalDistribution -> Time -> Double
cumulative :: TimeOfOriginNearCriticalDistribution -> Rate -> Rate
cumulative (TONCD Int
n' Rate
l Rate
m) Rate
t
  | Rate
t Rate -> Rate -> Bool
forall a. Ord a => a -> a -> Bool
<= Rate
0 = Rate
0
  | Bool
otherwise = Rate
t1 Rate -> Rate -> Rate
forall a. Num a => a -> a -> a
+ Rate
t2
  where
    d :: Rate
d = Rate
l Rate -> Rate -> Rate
forall a. Num a => a -> a -> a
- Rate
m
    n :: Rate
n = Int -> Rate
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n'
    t1 :: Rate
t1 = (Rate
t Rate -> Rate -> Rate
forall a. Num a => a -> a -> a
* Rate
l Rate -> Rate -> Rate
forall a. Fractional a => a -> a -> a
/ (Rate
1.0 Rate -> Rate -> Rate
forall a. Num a => a -> a -> a
+ Rate
t Rate -> Rate -> Rate
forall a. Num a => a -> a -> a
* Rate
l)) Rate -> Rate -> Rate
forall a. Floating a => a -> a -> a
** Rate
n
    t2 :: Rate
t2 = (Rate
n Rate -> Rate -> Rate
forall a. Num a => a -> a -> a
* Rate
t Rate -> Rate -> Rate
forall a. Num a => a -> a -> a
* Rate
t1) Rate -> Rate -> Rate
forall a. Num a => a -> a -> a
* Rate
d Rate -> Rate -> Rate
forall a. Fractional a => a -> a -> a
/ (Rate
2.0 Rate -> Rate -> Rate
forall a. Num a => a -> a -> a
* (Rate
1.0 Rate -> Rate -> Rate
forall a. Num a => a -> a -> a
+ Rate
t Rate -> Rate -> Rate
forall a. Num a => a -> a -> a
* Rate
l))

instance D.ContDistr TimeOfOriginNearCriticalDistribution where
  density :: TimeOfOriginNearCriticalDistribution -> Rate -> Rate
density = TimeOfOriginNearCriticalDistribution -> Rate -> Rate
density
  quantile :: TimeOfOriginNearCriticalDistribution -> Rate -> Rate
quantile = TimeOfOriginNearCriticalDistribution -> Rate -> Rate
quantile

-- | The density function Eq. (5).
density :: TimeOfOriginNearCriticalDistribution -> Time -> Double
density :: TimeOfOriginNearCriticalDistribution -> Rate -> Rate
density (TONCD Int
n' Rate
l Rate
m) Rate
t
  | Rate
t Rate -> Rate -> Bool
forall a. Ord a => a -> a -> Bool
< Rate
0 = Rate
0
  | Bool
otherwise = Rate
nom Rate -> Rate -> Rate
forall a. Fractional a => a -> a -> a
/ Rate
den
  where
    n :: Rate
n = Int -> Rate
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n'
    nom :: Rate
nom =
      Rate
n Rate -> Rate -> Rate
forall a. Num a => a -> a -> a
* (Rate
t Rate -> Rate -> Rate
forall a. Num a => a -> a -> a
* Rate
l Rate -> Rate -> Rate
forall a. Fractional a => a -> a -> a
/ (Rate
1 Rate -> Rate -> Rate
forall a. Num a => a -> a -> a
+ Rate
t Rate -> Rate -> Rate
forall a. Num a => a -> a -> a
* Rate
l)) Rate -> Rate -> Rate
forall a. Floating a => a -> a -> a
** Rate
n Rate -> Rate -> Rate
forall a. Num a => a -> a -> a
* (Rate
2 Rate -> Rate -> Rate
forall a. Num a => a -> a -> a
+ (Rate
3 Rate -> Rate -> Rate
forall a. Num a => a -> a -> a
+ Rate
n) Rate -> Rate -> Rate
forall a. Num a => a -> a -> a
* Rate
t Rate -> Rate -> Rate
forall a. Num a => a -> a -> a
* Rate
l Rate -> Rate -> Rate
forall a. Num a => a -> a -> a
- (Rate
1 Rate -> Rate -> Rate
forall a. Num a => a -> a -> a
+ Rate
n) Rate -> Rate -> Rate
forall a. Num a => a -> a -> a
* Rate
t Rate -> Rate -> Rate
forall a. Num a => a -> a -> a
* Rate
m)
    den :: Rate
den = Rate
2 Rate -> Rate -> Rate
forall a. Num a => a -> a -> a
* Rate
t Rate -> Rate -> Rate
forall a. Num a => a -> a -> a
* (Rate
1 Rate -> Rate -> Rate
forall a. Num a => a -> a -> a
+ Rate
t Rate -> Rate -> Rate
forall a. Num a => a -> a -> a
* Rate
l) Rate -> Rate -> Rate
forall a. Floating a => a -> a -> a
** Rate
2

-- | The inverted cumulative probability distribution 'cumulative'. See also
-- 'D.ContDistr'.
quantile :: TimeOfOriginNearCriticalDistribution -> Double -> Time
quantile :: TimeOfOriginNearCriticalDistribution -> Rate -> Rate
quantile (TONCD Int
n' Rate
l Rate
m) Rate
p
  | Rate
p Rate -> Rate -> Bool
forall a. Ord a => a -> a -> Bool
>= Rate
0 Bool -> Bool -> Bool
&& Rate
p Rate -> Rate -> Bool
forall a. Ord a => a -> a -> Bool
<= Rate
1 =
    Rate
t1 Rate -> Rate -> Rate
forall a. Num a => a -> a -> a
+ Rate
t2nom Rate -> Rate -> Rate
forall a. Fractional a => a -> a -> a
/ Rate
t2den
  | Bool
otherwise =
    [Char] -> Rate
forall a. HasCallStack => [Char] -> a
error ([Char] -> Rate) -> [Char] -> Rate
forall a b. (a -> b) -> a -> b
$
      [Char]
"PointProcess.quantile: p must be in [0,1] range. Got: "
        [Char] -> [Char] -> [Char]
forall a. [a] -> [a] -> [a]
++ Rate -> [Char]
forall a. Show a => a -> [Char]
show Rate
p
        [Char] -> [Char] -> [Char]
forall a. [a] -> [a] -> [a]
++ [Char]
"."
  where
    n :: Rate
n = Int -> Rate
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n'
    t1 :: Rate
t1 = - Rate
p Rate -> Rate -> Rate
forall a. Floating a => a -> a -> a
** (Rate
1 Rate -> Rate -> Rate
forall a. Fractional a => a -> a -> a
/ Rate
n) Rate -> Rate -> Rate
forall a. Fractional a => a -> a -> a
/ ((-Rate
1 Rate -> Rate -> Rate
forall a. Num a => a -> a -> a
+ Rate
p Rate -> Rate -> Rate
forall a. Floating a => a -> a -> a
** (Rate
1 Rate -> Rate -> Rate
forall a. Fractional a => a -> a -> a
/ Rate
n)) Rate -> Rate -> Rate
forall a. Num a => a -> a -> a
* Rate
l)
    t2nom :: Rate
t2nom = Rate
p Rate -> Rate -> Rate
forall a. Floating a => a -> a -> a
** (Rate
2 Rate -> Rate -> Rate
forall a. Fractional a => a -> a -> a
/ Rate
n) Rate -> Rate -> Rate
forall a. Num a => a -> a -> a
* (Rate
m Rate -> Rate -> Rate
forall a. Num a => a -> a -> a
- Rate
l)
    t2den :: Rate
t2den = Rate
2 Rate -> Rate -> Rate
forall a. Num a => a -> a -> a
* (-Rate
1 Rate -> Rate -> Rate
forall a. Num a => a -> a -> a
+ Rate
p Rate -> Rate -> Rate
forall a. Floating a => a -> a -> a
** (Rate
1 Rate -> Rate -> Rate
forall a. Fractional a => a -> a -> a
/ Rate
n)) Rate -> Rate -> Rate
forall a. Floating a => a -> a -> a
** Rate
2 Rate -> Rate -> Rate
forall a. Num a => a -> a -> a
* Rate
l Rate -> Rate -> Rate
forall a. Floating a => a -> a -> a
** Rate
2

instance D.ContGen TimeOfOriginNearCriticalDistribution where
  genContVar :: TimeOfOriginNearCriticalDistribution -> Gen (PrimState m) -> m Rate
genContVar = TimeOfOriginNearCriticalDistribution -> Gen (PrimState m) -> m Rate
forall d (m :: * -> *).
(ContDistr d, PrimMonad m) =>
d -> Gen (PrimState m) -> m Rate
D.genContinuous