module ELynx.Data.MarkovProcess.GammaRateHeterogeneity
( summarizeGammaRateHeterogeneity,
expand,
)
where
import qualified Data.ByteString.Lazy.Char8 as BL
import Data.List.NonEmpty (NonEmpty)
import qualified Data.List.NonEmpty as N
import qualified ELynx.Data.MarkovProcess.MixtureModel as M
import qualified ELynx.Data.MarkovProcess.PhyloModel as P
import qualified ELynx.Data.MarkovProcess.SubstitutionModel as S
import Numeric.Integration.TanhSinh
import Statistics.Distribution
import Statistics.Distribution.Gamma
import Prelude hiding (repeat)
summarizeGammaRateHeterogeneity :: Int -> Double -> [BL.ByteString]
summarizeGammaRateHeterogeneity n alpha =
map
BL.pack
[ "Discrete gamma rate heterogeneity.",
"Number of categories: " ++ show n,
"Shape parameter of gamma distribution: " ++ show alpha,
"Rates: " ++ show (getMeans n alpha)
]
expand :: Int -> Double -> P.PhyloModel -> P.PhyloModel
expand n alpha (P.SubstitutionModel sm) =
P.MixtureModel $ expandSubstitutionModel n alpha sm
expand n alpha (P.MixtureModel mm) =
P.MixtureModel $ expandMixtureModel n alpha mm
getName :: Int -> Double -> String
getName n alpha =
" with discrete gamma rate heterogeneity; "
++ show n
++ " categories; "
++ "shape parameter "
++ show alpha
splitSubstitutionModel ::
Int -> Double -> S.SubstitutionModel -> NonEmpty S.SubstitutionModel
splitSubstitutionModel n alpha sm = renamedSMs
where
means = getMeans n alpha
scaledSMs = N.map (`S.scale` sm) means
names = N.fromList $ map (("; gamma rate category " ++) . show) [1 :: Int ..]
renamedSMs = N.zipWith S.appendName names scaledSMs
expandSubstitutionModel ::
Int -> Double -> S.SubstitutionModel -> M.MixtureModel
expandSubstitutionModel n alpha sm = M.fromSubstitutionModels name ws sms
where
name = S.name sm <> getName n alpha
ws = N.repeat 1.0
sms = splitSubstitutionModel n alpha sm
expandMixtureModel :: Int -> Double -> M.MixtureModel -> M.MixtureModel
expandMixtureModel n alpha mm = M.concatenate name renamedMMs
where
name = M.name mm <> getName n alpha
means = getMeans n alpha
scaledMMs = N.map (`M.scale` mm) means
names = N.fromList $ map (("; gamma rate category " ++) . show) [1 :: Int ..]
renamedMMs = N.zipWith M.appendNameComponents names scaledMMs
getMeans :: Int -> Double -> NonEmpty Double
getMeans n alpha = N.fromList $ means <> pure lastMean
where
gamma = gammaDistr alpha (1.0 / alpha)
quantiles =
[quantile gamma (fromIntegral i / fromIntegral n) | i <- [0 .. n]]
meanFunc x = fromIntegral n * x * density gamma x
means =
[ integralAToB meanFunc (quantiles !! i) (quantiles !! (i + 1))
| i <- [0 .. n - 2]
]
lastMean = integralAToInf meanFunc (quantiles !! (n - 1))
eps :: Double
eps = 1e-8
method :: (Double -> Double) -> Double -> Double -> [Result]
method = parSimpson
integralAToB :: (Double -> Double) -> Double -> Double -> Double
integralAToB f a b = result . absolute eps $ method f a b
integralAToInf :: (Double -> Double) -> Double -> Double
integralAToInf f a =
(result . absolute eps $ nonNegative method f) - integralAToB f 0 a