Edit on GitHub

Sample from the prior using rjMCMC

Running an MCMC algorithm without data should result in a sample from the prior. This is commonly used in complicated Bayesian inference settings such as those arising in phylogenetics to verify correct implemetation of an MCMC algorithm. Of course, this cannot indicate problems in the likelihood implementation! (which can also be quite complicated in phylogenetic applications, at least in my humble opinion).

using Distributions, Beluga, Plots, StatsPlots, LaTeXStrings, Random
Random.seed!(230364)
Random.MersenneTwister(UInt32[0x000383dc], Random.DSFMT.DSFMT_state(Int32[1737833956, 1073610323, -1366922390, 1073425006, 2022353993, 1073475382, 140623403, 1073014326, -806344190, 1073230602  …  -1107069248, 1072895613, -717247506, 1073322432, -293592639, 600479730, -805356924, -698342859, 382, 0]), [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], UInt128[0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000  …  0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000], 1002, 0)

Obtain the species tree

nw = readline(joinpath(@__DIR__, "../../example/9dicots/9dicots.nw"))
d, p = DLWGD(nw, 1., 1., 0.9)
(model = DLWGD{Float64,NewickTree.TreeNode{Beluga.Branch{Float64}}}(17), data = PArray{Float64}(1))

Calling DLWGD without a data frame with gene counts will result in a 'mock' phylogenetic profile for which the likelihood defaults to 1 (or log-likelihood to 0).

logpdf!(d, p)
0.0

The prior is specified here

prior = IRRevJumpPrior(
    Ψ=[1 0.0 ; 0.0 1],
    X₀=MvNormal(log.([3,3]), [0.2 0.; 0. 0.2]),
    πK=DiscreteUniform(0, 10),
    πq=Beta(1,3),
    πη=Beta(3,1),
    Tl=treelength(d))
IRRevJumpPrior{Distributions.MvNormal{Float64,PDMats.PDMat{Float64,Array{Float64,2}},Array{Float64,1}},Distributions.Beta{Float64},Distributions.Beta{Float64},Distributions.DiscreteUniform,Nothing}
  Ψ: Array{Float64}((2, 2)) [1.0 0.0; 0.0 1.0]
  X₀: Distributions.MvNormal{Float64,PDMats.PDMat{Float64,Array{Float64,2}},Array{Float64,1}}
  πη: Distributions.Beta{Float64}
  πq: Distributions.Beta{Float64}
  πK: Distributions.DiscreteUniform
  πE: Nothing nothing
  equal: Bool false
  Tl: Float64 0.8565633076700001

Ψ is the prior covariance matrix for the Inverse-Wishart distribution. X₀ is the multivariate Normal prior on the mean duplication and loss rates. πK is the prior on the number of WGDs (i.e. the model indicator). πq is the Beta prior on the retention rates (iid). πη is the hyperprior on the parameter of the geometric distribution on the number of ancestral lineages at the root. Tl is the tree length (and is used for the prior on the WGD ages).

We have to set the rjMCMC kernel,

kernel = SimpleKernel(qkernel=Beta(1,3))
SimpleKernel
  qkernel: Distributions.Beta{Float64}
  accepted: Int64 0

construct the chain

chain = RevJumpChain(data=p, model=d, prior=prior, kernel=kernel)
init!(chain)

and sample

@time rjmcmc!(chain, 10000, trace=5, show=1000)
|    gen, pjmp,  k,k2,k3,    λ1,    λ2,    λ3,    μ1,    μ2,    μ3 …       logp,    logπ,     η
|   1000, 0.89, 10, 0, 3, 3.679, 2.474, 1.971, 2.619, 4.628, 3.272 ⋯      0.000, -42.678, 0.750
|   2000, 0.90,  7, 0, 1, 3.366,16.522,17.422, 2.809, 2.356, 0.795 ⋯      0.000, -66.850, 0.857
|   3000, 0.90,  6, 0, 0, 3.269, 0.719, 3.162, 4.603, 4.007, 4.048 ⋯      0.000, -43.807, 0.233
|   4000, 0.90,  4, 0, 1, 5.436, 4.459, 4.800, 4.650, 4.065, 4.163 ⋯      0.000, -31.587, 0.502
|   5000, 0.90,  8, 0, 1, 2.046, 3.237, 2.665, 4.460, 5.363, 2.958 ⋯      0.000, -21.571, 0.878
|   6000, 0.91,  8, 0, 0, 3.493, 4.696, 8.115, 2.451,14.995, 1.378 ⋯      0.000, -52.589, 0.488
|   7000, 0.91,  5, 0, 0, 2.147, 1.040, 3.344, 2.449, 5.052, 4.167 ⋯      0.000, -51.273, 0.894
|   8000, 0.90,  6, 0, 3, 1.998,211.492, 5.218, 4.013, 0.072, 3.319 ⋯      0.000, -68.781, 0.572
|   9000, 0.90,  1, 0, 0, 3.396, 2.900, 4.939, 3.159, 3.066, 3.040 ⋯      0.000, -29.138, 0.645
|  10000, 0.90,  3, 0, 0, 4.702, 2.710, 7.086, 2.133, 5.059, 1.350 ⋯      0.000, -41.516, 0.737
 76.195790 seconds (93.53 M allocations: 7.661 GiB, 2.77% gc time)

Now plot some stuff

p1 = bar(Beluga.freqmap(chain.trace[!,:k]), color=:white, title=L"k")
plot!(p1, prior.πK, color=:black, marker=nothing)
p2 = stephist(log.(chain.trace[!,:λ1]), fill=true, color=:black, alpha=0.2, normalize=true)
plot!(p2, Normal(log(3),√0.2), color=:black, linewidth=2, title=L"\lambda")
stephist!(log.(chain.trace[!,:λ9]), fill=true, color=:firebrick, alpha=0.2, normalize=true)
p3 = stephist(log.(chain.trace[!,:μ1]), fill=true, color=:black, alpha=0.2, normalize=true)
plot!(p3, Normal(log(3),√0.2), color=:black, linewidth=2, title=L"\mu")
stephist!(log.(chain.trace[!,:μ4]), fill=true, color=:firebrick, alpha=0.2, normalize=true)
p4 = stephist(chain.trace[!,:η1], color=:black, fill=true, alpha=0.2, normalize=true)
plot!(prior.πη, linewidth=2, color=:black, title=L"\eta")
plot(p1, p2, p3, p4, layout=(2,2), legend=false, grid=false, title_loc=:left)

Note that the prior for duplication and loss rates can have very heavy tails under the molecular clock priors (see the distributions in red as examples).

This page was generated using Literate.jl.