module ELynx.Data.MarkovProcess.GammaRateHeterogeneity
( summarizeGammaRateHeterogeneity
, expand
) where
import Control.Lens
import qualified Data.ByteString.Lazy.Char8 as L
import Numeric.Integration.TanhSinh
import Statistics.Distribution
import Statistics.Distribution.Gamma
import qualified ELynx.Data.MarkovProcess.MixtureModel as M
import qualified ELynx.Data.MarkovProcess.PhyloModel as P
import qualified ELynx.Data.MarkovProcess.SubstitutionModel as S
summarizeGammaRateHeterogeneity :: Int -> Double -> [L.ByteString]
summarizeGammaRateHeterogeneity n alpha = map L.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 -> [S.SubstitutionModel]
splitSubstitutionModel n alpha sm = renamedSMs
where
means = getMeans n alpha
scaledSMs = map (`S.scale` sm) means
names = map (("; gamma rate category " ++) . show) [1 :: Int ..]
renamedSMs = zipWith S.appendName names scaledSMs
expandSubstitutionModel :: Int -> Double -> S.SubstitutionModel -> M.MixtureModel
expandSubstitutionModel n alpha sm = M.fromSubstitutionModels name (repeat 1.0) sms
where
name = sm ^. S.name <> getName n alpha
sms = splitSubstitutionModel n alpha sm
expandMixtureModel :: Int -> Double -> M.MixtureModel -> M.MixtureModel
expandMixtureModel n alpha mm = M.concatenate name renamedMMs
where
name = mm ^. M.name <> getName n alpha
means = getMeans n alpha
scaledMMs = map (`M.scale` mm) means
names = map (("; gamma rate category " ++) . show) [1 :: Int ..]
renamedMMs = zipWith M.appendName names scaledMMs
getMeans :: Int -> Double -> [Double]
getMeans n alpha = means ++ 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-6
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