# Unable to specify a correct prior distribution

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

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

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

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.

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

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.

I see, thank you for the confirmation!

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

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

This should be fixed on the new release @anqipan