Unable to specify a correct prior distribution

estimation bayesian
#1

Hi,
I am having difficulty specifying a prior distribution in Pumas with lower bound and initials, my codes is:

    @param   begin

      θ ~ 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,-Inf,-Inf,-Inf,-Inf],init=[18,14.2,2.13,4.29,1,1,1,1])

      Ω  ~ Wishart(3, Matrix{Float64}([0.282	0	0; 0	0.149	0; 0	0	0.0377]))

      σ_prop ∈ RealDomain(lower = 0)

    end

I mainly had 2 error message:
1.

LoadError: MethodError: no method matching MvNormal(::Array{Float64,1}, ::Array{Float64,2}; lower=[0.0, 0.0, 0.0, 0.0, -Inf, -Inf, -Inf, -Inf])

LoadError: ArgumentError: no Domain(::Wishart{Float64,PDMats.PDMat{Float64,Array{Float64,2}}}) constructor defined

I have google both of them and didn’t find anything matched…

Basically I am trying to specify a multivariate normal distribution for the Thetas with lower bound and initials, a Wishart distribution for the Omegas with initials.

May I have any recommendations on which part I missed in my codes?

Thanks!
Anqi

Unable to interpret the fitting error
#2

Hello,
Just follow up with the complete coding for this topic:

@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,-Inf,-Inf,-Inf,-Inf])


      Ω  ~ 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

inputDataset = CSV.read("dat.csv")
df = read_pumas(inputDataset, id = :ID, dvs =[:DV],  cvs=[:WTKG, :CRCL], evid=:EVID, amt=:AMT,cmt=:CMT, rate=:RATE, time=:TIME)


param = (
    θ=([18,14.2,2.13,4.29,1,1,1,1]),
    Ω  = Diagonal([0.25,0.25,0.25]),
    σ_prop = 0.04)

 pkres = fit(ped_mod, df, param, Pumas.FOCEI(),
              optimize_fn=Pumas.DefaultOptimizeFN(show_trace=true, extended_trace=false))

The error is saying no method matching inverse:

LoadError: MethodError: no method matching inverse!(::SubArray{Float64,1,Array{Float64,1},Tuple{UnitRange{Int64}},true}, ::Pumas.PSDTransform, ::Diagonal{Float64,Array{Float64,1}})

Closest candidates are:

inverse!(::AbstractArray{T,1} where T<:Real, !Matched::TransformVariables.ArrayTransform, ::AbstractArray) at C:\Users\.juliapro\JuliaPro_v1.4.2-2\packages\TransformVariables\a4AMY\src\aggregation.jl:79

inverse!(::AbstractArray, !Matched::Pumas.ElementArrayTransform, ::AbstractArray) at C:\Users\.juliapro\JuliaPro_v1.4.2-2\packages\Pumas\NCmSe\src\estimation\transforms.jl:25

inverse!(::AbstractArray, !Matched::Pumas.ConstantTransform, ::Any) at C:\Users\juliapro\JuliaPro_v1.4.2-2\packages\Pumas\NCmSe\src\estimation\transforms.jl:54

...

in expression starting at C:\Users\118

_inverse!_tuple(::Array{Float64,1}, ::Tuple{Pumas.ElementArrayTransform{TransformVariables.ShiftedExp{true,Float64},1},Pumas.PSDTransform,TransformVariables.ShiftedExp{true,Int64}}, ::Tuple{Array{Float64,1},Diagonal{Float64,Array{Float64,1}},Float64}) at aggregation.jl:196

inverse!(::Array{Float64,1}, ::TransformVariables.TransformTuple{NamedTuple{(:θ, :Ω, :σ_prop),Tuple{Pumas.ElementArrayTransform{TransformVariables.ShiftedExp{true,Float64},1},Pumas.PSDTransform,TransformVariables.ShiftedExp{true,Int64}}}}, ::NamedTuple{(:θ, :Ω, :σ_prop),Tuple{Array{Float64,1},Diagonal{Float64,Array{Float64,1}},Float64}}) at aggregation.jl:237

inverse(::TransformVariables.TransformTuple{NamedTuple{(:θ, :Ω, :σ_prop),Tuple{Pumas.ElementArrayTransform{TransformVariables.ShiftedExp{true,Float64},1},Pumas.PSDTransform,TransformVariables.ShiftedExp{true,Int64}}}}, ::NamedTuple{(:θ, :Ω, :σ_prop),Tuple{Array{Float64,1},Diagonal{Float64,Array{Float64,1}},Float64}}) at generic.jl:206

fit(::PumasModel{ParamSet{NamedTuple{(:θ, :Ω, :σ_prop),Tuple{Constrained{MvNormal{Float64,PDMats.PDMat{Float64,Array{Float64,2}},Array{Float64,1}},VectorDomain{Array{Float64,1},Array{TransformVariables.Infinity{true},1},Array{Float64,1}}},InverseWishart{Float64,PDMats.PDMat{Float64,Array{Float64,2}}},RealDomain{Int64,TransformVariables.Infinity{true},Float64}}}},var"#548#603",var"#549#604",var"#550#606",Central1Periph1,var"#551#607",var"#577#633"}, ::Array{Subject{NamedTuple{(:DV,),Tuple{Array{Union{Missing, Float64},1}}},Pumas.ConstantCovar{NamedTuple{(:WTKG, :CRCL),Tuple{Float64,Float64}}},Array{Pumas.Event{Float64,Float64,Float64,Float64,Float64,Float64,Int64},1},Array{Float64,1},Float64},1}, ::NamedTuple{(:θ, :Ω, :σ_prop),Tuple{Array{Float64,1},Diagonal{Float64,Array{Float64,1}},Float64}}, ::Pumas.FOCEI; optimize_fn::Pumas.DefaultOptimizeFN{Nothing,NamedTuple{(:show_trace, :store_trace, :extended_trace, :g_tol, :allow_f_increases),Tuple{Bool,Bool,Bool,Float64,Bool}}}, constantcoef::NamedTuple{,Tuple{}}, omegas::Tuple{}, ensemblealg::EnsembleSerial, kwargs::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{,Tuple{}}}) at likelihoods.jl:1330

(::StatsBase.var"#fit##kw")(::NamedTuple{(:optimize_fn,),Tuple{Pumas.DefaultOptimizeFN{Nothing,NamedTuple{(:show_trace, :store_trace, :extended_trace, :g_tol, :allow_f_increases),Tuple{Bool,Bool,Bool,Float64,Bool}}}}}, ::typeof(fit), ::PumasModel{ParamSet{NamedTuple{(:θ, :Ω, :σ_prop),Tuple{Constrained{MvNormal{Float64,PDMats.PDMat{Float64,Array{Float64,2}},Array{Float64,1}},VectorDomain{Array{Float64,1},Array{TransformVariables.Infinity{true},1},Array{Float64,1}}},InverseWishart{Float64,PDMats.PDMat{Float64,Array{Float64,2}}},RealDomain{Int64,TransformVariables.Infinity{true},Float64}}}},var"#548#603",var"#549#604",var"#550#606",Central1Periph1,var"#551#607",var"#577#633"}, ::Array{Subject{NamedTuple{(:DV,),Tuple{Array{Union{Missing, Float64},1}}},Pumas.ConstantCovar{NamedTuple{(:WTKG, :CRCL),Tuple{Float64,Float64}}},Array{Pumas.Event{Float64,Float64,Float64,Float64,Float64,Float64,Int64},1},Array{Float64,1},Float64},1}, ::NamedTuple{(:θ, :Ω, :σ_prop),Tuple{Array{Float64,1},Diagonal{Float64,Array{Float64,1}},Float64}}, ::Pumas.FOCEI) at likelihoods.jl:1328

top-level scope at tp4_5.jl:118

Thanks!
Anqi

#3

Error 2. in your first post is a missing feature in Pumas but it’s tiny so we might be able to consider it a bugfix and include it Pumas 1.0. Right now, we only support the InverseWishart but extending the support to cover Wishart should be very easy.

The problem in your followup post is Ω = Diagonal([0.25,0.25,0.25]) because it restricts Ω do a diagonal matrix but the InverseWishart isn’t restricted to diagonal matrices so you should be able to make this work by setting Ω = diagm([0.25,0.25,0.25]) instead. We might be able to make your version to work as well. I’ve filed an issue.

#4

Thank you! The error is fixed after I changed the initial to a diagm.

And may I just confirm that, earlier this morning you mentioned the prior is just a domain(at least for the current version), which sounds it won’t enter the posterior distribution like a real prior.
Is this the case for both FOCE and HMC-NUTS?
Or the distributional domain(prior) does enter the posterior as a real prior in HMC-NUTS?

Thanks!
Anqi

#5

And may I just confirm that, earlier this morning you mentioned the prior is just a domain(at least for the current version), which sounds it won’t enter the posterior distribution like a real prior.
Is this the case for both FOCE and HMC-NUTS?

Yes. That is indeed the case for FOCE(I) (and LaplaceI) but for the NUTS sampler aka BayesMCMC the likelihood contribution of the prior is applied as you would expect.

#6

I see, thank you for the confirmation!

#7

I do have a further question regard the same model. Is there anyway that I could fix the off-diagonal elements of OMEGA matrix to 0?
I tried to use Diagonal() to force it as diagonal matrix both in @param and the initials, but it doesn’t run through.

Thank you!
Anqi

#8

I think the only way to do that would be to specify each of the diagonal elements separately, i.e. something like

@random begin
  ηCL  ~ InverseGamma(...)
  ηV1  ~ InverseGamma(...)
  ηV2  ~ InverseGamma(...)
end

I chose InverseGamma here since it’s equivalent to a scalar InverseWishart but it might be better to use the Gamma since it has support at zero.

1 Like