# Bayesian Posterior convergence

**URL:** <https://discourse.pumas.ai/t/bayesian-posterior-convergence/381>\
**Category:** Diagnostics\
**Created:** [August 4, 2020, 2:46pm UTC](https://discourse.pumas.ai/t/bayesian-posterior-convergence/381 "2020-08-04T14:46:21Z")\
**Posts on this page:** 6\
**Page:** 1

<div class="post-metadata">

**Author:** ![anqipan](https://avatars.discourse-cdn.com/v4/letter/a/ecae2f/32.png) [@anqipan](https://discourse.pumas.ai/u/anqipan)\
**Post date:** [August 4, 2020, 2:46pm UTC](https://discourse.pumas.ai/t/bayesian-posterior-convergence/381/1 "2020-08-04T14:46:21Z")

</div>

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](https://canada1.discourse-cdn.com/flex030/uploads/pumas/original/1X/858346db57d89ff86b0aa96ed97d244619fc84c3.png)

---

<div class="post-metadata">

**Author:** ![andreasnoack](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.pumas.ai/andreasnoack/32/22_2.png) [@andreasnoack](https://discourse.pumas.ai/u/andreasnoack)\
**Post date:** [August 5, 2020, 9:21am UTC](https://discourse.pumas.ai/t/bayesian-posterior-convergence/381/2 "2020-08-05T09:21:57Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![anqipan](https://avatars.discourse-cdn.com/v4/letter/a/ecae2f/32.png) [@anqipan](https://discourse.pumas.ai/u/anqipan)\
**Post date:** [August 5, 2020, 2:32pm UTC](https://discourse.pumas.ai/t/bayesian-posterior-convergence/381/3 "2020-08-05T14:32:12Z")

</div>

Hi Andreas,

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

Thanks!  
Anqi

```auto
@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)

```

---

<div class="post-metadata">

**Author:** ![mohamed82008](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.pumas.ai/mohamed82008/32/127_2.png) [@mohamed82008](https://discourse.pumas.ai/u/mohamed82008)\
**Post date:** [August 6, 2020, 2:05pm UTC](https://discourse.pumas.ai/t/bayesian-posterior-convergence/381/4 "2020-08-06T14:05:22Z")

</div>

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 🙂

---

<div class="post-metadata">

**Author:** ![anqipan](https://avatars.discourse-cdn.com/v4/letter/a/ecae2f/32.png) [@anqipan](https://discourse.pumas.ai/u/anqipan)\
**Post date:** [August 7, 2020, 6:00pm UTC](https://discourse.pumas.ai/t/bayesian-posterior-convergence/381/5 "2020-08-07T18:00:36Z")

</div>

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

Many thanks again!

Best,  
Anqi

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.pumas.ai/chrisrackauckas/32/14_2.png) [@ChrisRackauckas](https://discourse.pumas.ai/u/ChrisRackauckas)\
**Post date:** [August 12, 2020, 6:46am UTC](https://discourse.pumas.ai/t/bayesian-posterior-convergence/381/6 "2020-08-12T06:46:05Z")

</div>

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.