Bayesian Posterior convergence

Diagnostics
1

Hi,
I applied NUTS on a 2 cmt model and obtain the posterior samples with trace plot pasted below. I am wondering why the samples “stuck” at a estimated value after the warm-up period. And from the Bayesian example I have, they showed the same behavior.
NUTS returned saying that the average acceptance rate is 0.88, but according to the trace plot, it seems the acceptance rate is 0 after warm-up.

Could I check whether this trace plot shows convergence or something is wrong with that?

Thanks!
Anqi
Picture1

2

It does look odd. Would you be able to share the model?

3

Hi Andreas,

Here attached the codes. And I didn’t have any error or warning message.

Thanks!
Anqi

@time using Distributions, Pumas, CSV, TableView, StatsPlots, DataFrames, LinearAlgebra

ped_mod = @model begin

    @param   begin
      θ ~ Constrained(MvNormal([16.6,13.1,4.05,4.89,0.623,0,0,0],
      Matrix{Float64}([0.0895	0.0167	0.00318	0.0033	0.000907	0	0	0;
                       0.0167	0.0655	-0.00765	0.00333	-0.000917	0	0	0;
                       0.00318	-0.00765	0.0238	0.0147	-0.000225	0	0	0;
                       0.0033	0.00333	0.0147	0.0161	-0.000243	0	0	0;
                       0.000907	-0.000917	-0.000225	-0.000243	0.0016	0	0	0;
                       0	0	0	0	0	1000000	0	0;
                       0	0	0	0	0	0	1000000	0;
                       0	0	0	0	0	0	0	1000000])),lower=[0,0,0,0,0,0,0,0])

      Ω ~ InverseWishart(3, diagm([0.282, 0.149, 0.0377]))
      σ_prop ∈ RealDomain(lower = 0)
    end

    @random begin
      η ~ MvNormal(Ω)
    end

    @covariates WTKG  CRCL

    @pre begin
    #--INSERT COVARIATE EFFECTS
    COV1=(CRCL/100)^θ[5]
    COV2=(WTKG/70)^θ[6]
    COV3=(WTKG/70)^θ[7]
    COV4=(WTKG/70)^θ[8]

    TVCLI = θ[1]*COV1*COV4
    TVCL  = TVCLI

    TVV1I = θ[2]*COV2
    TVV1  = TVV1I

    TVQI  = θ[3]
    TVQ   = TVQI

    TVV2I = θ[4]*COV3
    TVV2  = TVV2I

    CL  = TVCL*exp(η[1])
    V1  = TVV1*exp(η[2])
    Q   = TVQ
    V2  = TVV2*exp(η[3])

    Vc = V1
    Vp = V2

    S1 = V1

    #CALCULATION OF SECONDARY PARAMETERS
    KE  = CL/V1
    K12 = Q/V1
    K21 = Q/V2
    AA = KE+K12+K21
    ALPH = (AA+sqrt(AA*AA-4*KE*K21))/2
    BETA = (AA-sqrt(AA*AA-4*KE*K21))/2
    end

    @dynamics Central1Periph1 #a two compartment model

    @derived begin
        cp = @. 1000*(Central/Vc)
        DV ~ @. Normal(cp,sqrt((cp^2*σ_prop)))
      end
  end

  ## Read Dataset
  inputDataset = CSV.read("dat.csv")

  ## Pumas expects a dvs and id column. So map the names of the columns
  df = read_pumas(inputDataset, id = :ID, dvs =[:DV],  cvs=[:WTKG, :CRCL], evid=:EVID, amt=:AMT,cmt=:CMT, rate=:RATE, time=:TIME)

#Initla value for parameters
param = (
    θ=([18,14.2,2.13,4.29,1,1,1,1]),
    Ω  = diagm([0.25,0.25,0.25]),
    σ_prop = 0.04)

  pkres = fit(ped_mod, df, param, Pumas.BayesMCMC(),nsamples=10000)
4

Hi!

The covariance matrix you are using for the prior of θ looks pretty ill-conditioned. This means that the pdf curve is quite narrow and high along the dimensions with little variance but quite wide and low along the dimensions with high variance. This will make life difficult for NUTS when sampling. Consider re-parameterizing the model by separating the first 5 θ variables from the last 3. Also, it’s a good idea to use a truncated unit Gaussian prior for the last 3 θ variables, scaling them by 1000 after sampling. This makes the sampling space “nicer” as the random variables are then similarly scaled. If none of this works, please report again and I can take a closer look :slight_smile:

2 Likes
5

I see, thanks alot for the recommendations, I will try each of them to see how they work :slightly_smiling_face:

Many thanks again!

Best,
Anqi

1 Like
6

The other reason why this is occurring is because the target accept rate is hard coded to 0.8 but we should really make that flexible. I opened an issue.