 # Bayesian Posterior convergence

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 It does look odd. Would you be able to share the model?

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)^θ
COV2=(WTKG/70)^θ
COV3=(WTKG/70)^θ
COV4=(WTKG/70)^θ

TVCLI = θ*COV1*COV4
TVCL  = TVCLI

TVV1I = θ*COV2
TVV1  = TVV1I

TVQI  = θ
TVQ   = TVQI

TVV2I = θ*COV3
TVV2  = TVV2I

CL  = TVCL*exp(η)
V1  = TVV1*exp(η)
Q   = TVQ
V2  = TVV2*exp(η)

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

## 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)
``````

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 2 Likes

I see, thanks alot for the recommendations, I will try each of them to see how they work Many thanks again!

Best,
Anqi

1 Like

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.