Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

How to transform the parameters? #28

Closed
y1my1 opened this issue Apr 17, 2020 · 8 comments
Closed

How to transform the parameters? #28

y1my1 opened this issue Apr 17, 2020 · 8 comments

Comments

@y1my1
Copy link

y1my1 commented Apr 17, 2020

Hi Tamas,

Sorry to trouble you again. I am working on a problem that is defined this way

struct REProblem{Td, Tλ}
    "data of the model"
    data::Td
    "prior of the model"
    λ::Tλ
end

Right now I define the callable as the following

function (problem::REProblem)(θ)
    @unpack data, λ = problem
    @unpack Θ = θ
    return LogLikelihood(data, Θ) + LogPriorProbability(Θ, λ)
end

where the function LogLikelihood takes the data and Θ to calculate the log likelihood in a somewhat complicate process. And I also wrap up the calculation of prior probability into LogPriorProbability.

Since I have a lot of parameters to be estimated, I wrap up all the parameters as a struct

@with_kw struct DynamicParameters{F<:AbstractFloat} @deftype F
    ρ_z
    σ_μ

    γᴵ
    γᴰ
    γᶠ
    γˢ

    a₁
    a₂
    a₃
    a₄

    b₁
    b₂
    b₃
    b₄

    Φˣ
end

and Θ is an instance of DynamicParameters.
I wonder how can I transform the parameters in this model?

Or I could define the callable as this more expanded way

function (problem::REProblem)(θ)
    @unpack data, λ = problem
    @unpack ρ_z, σ_μ, γᴵ, γᴰ, γᶠ, γˢ = θ
    @unpack a₁, a₂, a₃, a₄, b₁, b₂, b₃, b₄, Φˣ = θ
    return LogLikelihood(data, ρ_z, σ_μ, γᴵ, γᴰ, γᶠ, γˢ,a₁, a₂, a₃, a₄, b₁, b₂, b₃, b₄, Φˣ) + 
LogPriorProbability(ρ_z, σ_μ, γᴵ, γᴰ, γᶠ, γˢ,a₁, a₂, a₃, a₄, b₁, b₂, b₃, b₄, Φˣ, λ)
end

Then, I could transform the parameters as

trans = as((ρ_z=as(Real, -1.0, 1.0), σ_μ = asℝ₊, ..., Φˣ=asℝ))

But I really want to define the callable in a more compact way, can you please give me some suggestions? Thanks a lot for your help.

@y1my1
Copy link
Author

y1my1 commented Apr 17, 2020

As I look at the problem again, or I may define the callable as this

function (problem::REProblem)(θ)
    @unpack data, λ = problem
    @unpack ρ_z, σ_μ, γᴵ, γᴰ, γᶠ, γˢ = θ
    @unpack a₁, a₂, a₃, a₄, b₁, b₂, b₃, b₄, Φˣ = θ
    return LogLikelihood(data, θ) + LogPriorProbability(θ, λ)
end

Then use the transformation above. Or you may have better suggestions? Thanks again.

@y1my1
Copy link
Author

y1my1 commented Apr 17, 2020

After learning from your model here, I redefine the callable this way

function (problem:: REProblem)(θ)
    @unpack data, ie, λ = problem                  # extract the data and priors
    @unpack ρ_z, σ_μ, γᴵ, γᴰ, γᶠ, γˢ = θ           # extract the dynamic parameters (part 1)
    @unpack a₀, a₁, a₂, a₃, b₀, b₁, b₂, b₃, Φˣ = θ  # extract the dynamic parameters (part 2)
    Θ = DynamicParameters(ρ_z,σ_μ,γᴵ,γᴰ,γᶠ,γˢ,a₀,a₁,a₂,a₃,b₀,b₁,b₂,b₃,Φˣ)
    return LogLikelihood(data, Θ, ie) + LogPriorProbability(Θ, λ)
end

And my problem definition has changed a little (adding one component)

struct REProblem{Td, Ti, Tλ}
    "data of the model"
    data::Td
    "intergration elements need for calculation the likelihood"
    ie::Ti
    "prior of the model"
    λ::Tλ
end

Then I forward like this way

   # define a model using data and priors
    p = RdExpProductivity(rdexp_data, ie, prior)
    # transform variables
    trans = as((ρ_z=as𝕀, σ_μ=asℝ, γᴵ=asℝ, γᴰ=asℝ, γᶠ=asℝ, γˢ=asℝ,
                a₀=asℝ, a₁=asℝ, a₂=asℝ, a₃=asℝ,
                b₀=asℝ, b₁=asℝ, b₂=asℝ, b₃=asℝ, Φˣ=asℝ))
    # transformed probability
    tP = TransformedLogDensity(trans, p)
    # use ForwardDiff for automatic differentiation (AD)
    ∇tP = ADgradient(:ForwardDiff, tP)

I tried to evaluate the callable p, it works fine

p((ρ_z=.76, σ_μ=-.21, γᴵ=70.0, γᴰ=380.0, γᶠ=10.0, γˢ=50.0, a₀=-4.0059, a₁=1.9955, a₂=0.1724, 
a₃=0.2710, b₀=-5.4041, b₁=2.2519, b₂=0.1270, b₃=0.3305, Φˣ=3.80))
-5974.952402085821

Then I ran my HMC

    # sampling using Dynamic Hamiltonian Monte Carlo, random stating point.
    results = mcmc_with_warmup(Random.GLOBAL_RNG, ∇tP, 1_000)

but got the following error

ERROR: MethodError: no method matching DynamicParameters(::ForwardDiff.Dual{ForwardDiff.Tag{LogDensityProblems.var"#34#35"{TransformedLogDensity{TransformVariables.TransformTuple{NamedTuple{(:ρ_z, :σ_μ, :γᴵ, :γᴰ, :γᶠ, :γˢ, :a₀, :a₁, :a₂, :a₃, :b₀, :b₁, :b₂, :b₃, :Φˣ),NTuple{15,TransformVariables.Identity}}},RdExpProductivity{RdExpData{Float64,Int64},IntegrationElements{Float64},Priors{Float64}}}},Float64},Float64,8}, ::ForwardDiff.Dual{ForwardDiff.Tag{LogDensityProblems.var"#34#35"{TransformedLogDensity{TransformVariables.TransformTuple{NamedTuple{(:ρ_z, :σ_μ, :γᴵ, :γᴰ, :γᶠ, :γˢ, :a₀, :a₁, :a₂, :a₃, :b₀, :b₁, :b₂, :b₃, :Φˣ),NTuple{15,TransformVariables.Identity}}},RdExpProductivity{RdExpData{Float64,Int64},IntegrationElements{Float64},Priors{Float64}}}},Float64},Float64,8}, ::ForwardDiff.Dual{ForwardDiff.Tag{LogDensityProblems.var"#34#35"{TransformedLogDensity{TransformVariables.TransformTuple{NamedTuple{(:ρ_z, :σ_μ, :γᴵ, :γᴰ, :γᶠ, :γˢ, :a₀, :a₁, :a₂, :a₃, :b₀, :b₁, :b₂, :b₃, :Φˣ),NTuple{15,TransformVariables.Identity}}},RdExpProductivity{RdExpData{Float64,Int64},IntegrationElements{Float64},Priors{Float64}}}},Float64},Float64,8}, ::ForwardDiff.Dual{ForwardDiff.Tag{LogDensityProblems.var"#34#35"{TransformedLogDensity{TransformVariables.TransformTuple{NamedTuple{(:ρ_z, :σ_μ, :γᴵ, :γᴰ, :γᶠ, :γˢ, :a₀, :a₁, :a₂, :a₃, :b₀, :b₁, :b₂, :b₃, :Φˣ),NTuple{15,TransformVariables.Identity}}},RdExpProductivity{RdExpData{Float64,Int64},IntegrationElements{Float64},Priors{Float64}}}},Float64},Float64,8}, ::ForwardDiff.Dual{ForwardDiff.Tag{LogDensityProblems.var"#34#35"{TransformedLogDensity{TransformVariables.TransformTuple{NamedTuple{(:ρ_z, :σ_μ, :γᴵ, :γᴰ, :γᶠ, :γˢ, :a₀, :a₁, :a₂, :a₃, :b₀, :b₁, :b₂, :b₃, :Φˣ),NTuple{15,TransformVariables.Identity}}},RdExpProductivity{RdExpData{Float64,Int64},IntegrationElements{Float64},Priors{Float64}}}},Float64},Float64,8}, ::ForwardDiff.Dual{ForwardDiff.Tag{LogDensityProblems.var"#34#35"{TransformedLogDensity{TransformVariables.TransformTuple{NamedTuple{(:ρ_z, :σ_μ, :γᴵ, :γᴰ, :γᶠ, :γˢ, :a₀, :a₁, :a₂, :a₃, :b₀, :b₁, :b₂, :b₃, :Φˣ),NTuple{15,TransformVariables.Identity}}},RdExpProductivity{RdExpData{Float64,Int64},IntegrationElements{Float64},Priors{Float64}}}},Float64},Float64,8}, ::ForwardDiff.Dual{ForwardDiff.Tag{LogDensityProblems.var"#34#35"{TransformedLogDensity{TransformVariables.TransformTuple{NamedTuple{(:ρ_z, :σ_μ, :γᴵ, :γᴰ, :γᶠ, :γˢ, :a₀, :a₁, :a₂, :a₃, :b₀, :b₁, :b₂, :b₃, :Φˣ),NTuple{15,TransformVariables.Identity}}},RdExpProductivity{RdExpData{Float64,Int64},IntegrationElements{Float64},Priors{Float64}}}},Float64},Float64,8}, ::ForwardDiff.Dual{ForwardDiff.Tag{LogDensityProblems.var"#34#35"{TransformedLogDensity{TransformVariables.TransformTuple{NamedTuple{(:ρ_z, :σ_μ, :γᴵ, :γᴰ, :γᶠ, :γˢ, :a₀, :a₁, :a₂, :a₃, :b₀, :b₁, :b₂, :b₃, :Φˣ),NTuple{15,TransformVariables.Identity}}},RdExpProductivity{RdExpData{Float64,Int64},IntegrationElements{Float64},Priors{Float64}}}},Float64},Float64,8}, ::ForwardDiff.Dual{ForwardDiff.Tag{LogDensityProblems.var"#34#35"{TransformedLogDensity{TransformVariables.TransformTuple{NamedTuple{(:ρ_z, :σ_μ, :γᴵ, :γᴰ, :γᶠ, :γˢ, :a₀, :a₁, :a₂, :a₃, :b₀, :b₁, :b₂, :b₃, :Φˣ),NTuple{15,TransformVariables.Identity}}},RdExpProductivity{RdExpData{Float64,Int64},IntegrationElements{Float64},Priors{Float64}}}},Float64},Float64,8}, ::ForwardDiff.Dual{ForwardDiff.Tag{LogDensityProblems.var"#34#35"{TransformedLogDensity{TransformVariables.TransformTuple{NamedTuple{(:ρ_z, :σ_μ, :γᴵ, :γᴰ, :γᶠ, :γˢ, :a₀, :a₁, :a₂, :a₃, :b₀, :b₁, :b₂, :b₃, :Φˣ),NTuple{15,TransformVariables.Identity}}},RdExpProductivity{RdExpData{Float64,Int64},IntegrationElements{Float64},Priors{Float64}}}},Float64},Float64,8}, ::ForwardDiff.Dual{ForwardDiff.Tag{LogDensityProblems.var"#34#35"{TransformedLogDensity{TransformVariables.TransformTuple{NamedTuple{(:ρ_z, :σ_μ, :γᴵ, :γᴰ, :γᶠ, :γˢ, :a₀, :a₁, :a₂, :a₃, :b₀, :b₁, :b₂, :b₃, :Φˣ),NTuple{15,TransformVariables.Identity}}},RdExpProductivity{RdExpData{Float64,Int64},IntegrationElements{Float64},Priors{Float64}}}},Float64},Float64,8}, ::ForwardDiff.Dual{ForwardDiff.Tag{LogDensityProblems.var"#34#35"{TransformedLogDensity{TransformVariables.TransformTuple{NamedTuple{(:ρ_z, :σ_μ, :γᴵ, :γᴰ, :γᶠ, :γˢ, :a₀, :a₁, :a₂, :a₃, :b₀, :b₁, :b₂, :b₃, :Φˣ),NTuple{15,TransformVariables.Identity}}},RdExpProductivity{RdExpData{Float64,Int64},IntegrationElements{Float64},Priors{Float64}}}},Float64},Float64,8}, ::ForwardDiff.Dual{ForwardDiff.Tag{LogDensityProblems.var"#34#35"{TransformedLogDensity{TransformVariables.TransformTuple{NamedTuple{(:ρ_z, :σ_μ, :γᴵ, :γᴰ, :γᶠ, :γˢ, :a₀, :a₁, :a₂, :a₃, :b₀, :b₁, :b₂, :b₃, :Φˣ),NTuple{15,TransformVariables.Identity}}},RdExpProductivity{RdExpData{Float64,Int64},IntegrationElements{Float64},Priors{Float64}}}},Float64},Float64,8}, ::ForwardDiff.Dual{ForwardDiff.Tag{LogDensityProblems.var"#34#35"{TransformedLogDensity{TransformVariables.TransformTuple{NamedTuple{(:ρ_z, :σ_μ, :γᴵ, :γᴰ, :γᶠ, :γˢ, :a₀, :a₁, :a₂, :a₃, :b₀, :b₁, :b₂, :b₃, :Φˣ),NTuple{15,TransformVariables.Identity}}},RdExpProductivity{RdExpData{Float64,Int64},IntegrationElements{Float64},Priors{Float64}}}},Float64},Float64,8}, ::ForwardDiff.Dual{ForwardDiff.Tag{LogDensityProblems.var"#34#35"{TransformedLogDensity{TransformVariables.TransformTuple{NamedTuple{(:ρ_z, :σ_μ, :γᴵ, :γᴰ, :γᶠ, :γˢ, :a₀, :a₁, :a₂, :a₃, :b₀, :b₁, :b₂, :b₃, :Φˣ),NTuple{15,TransformVariables.Identity}}},RdExpProductivity{RdExpData{Float64,Int64},IntegrationElements{Float64},Priors{Float64}}}},Float64},Float64,8})
Stacktrace:
 [1] (::RdExpProductivity{RdExpData{Float64,Int64},IntegrationElements{Float64},Priors{Float64}})(::NamedTuple{(:ρ_z, :σ_μ, :γᴵ, :γᴰ, :γᶠ, :γˢ, :a₀, :a₁, :a₂, :a₃, :b₀, :b₁, :b₂, :b₃, :Φˣ),NTuple{15,ForwardDiff.Dual{ForwardDiff.Tag{LogDensityProblems.var"#34#35"{TransformedLogDensity{TransformVariables.TransformTuple{NamedTuple{(:ρ_z, :σ_μ, :γᴵ, :γᴰ, :γᶠ, :γˢ, :a₀, :a₁, :a₂, :a₃, :b₀, :b₁, :b₂, :b₃, :Φˣ),NTuple{15,TransformVariables.Identity}}},RdExpProductivity{RdExpData{Float64,Int64},IntegrationElements{Float64},Priors{Float64}}}},Float64},Float64,8}}}) at /Users/jack/Programming/rdexp_hmc/code/mcmc_hamiltonian.jl:30
 [2] transform_logdensity(::TransformVariables.TransformTuple{NamedTuple{(:ρ_z, :σ_μ, :γᴵ, :γᴰ, :γᶠ, :γˢ, :a₀, :a₁, :a₂, :a₃, :b₀, :b₁, :b₂, :b₃, :Φˣ),NTuple{15,TransformVariables.Identity}}}, ::RdExpProductivity{RdExpData{Float64,Int64},IntegrationElements{Float64},Priors{Float64}}, ::Array{ForwardDiff.Dual{ForwardDiff.Tag{LogDensityProblems.var"#34#35"{TransformedLogDensity{TransformVariables.TransformTuple{NamedTuple{(:ρ_z, :σ_μ, :γᴵ, :γᴰ, :γᶠ, :γˢ, :a₀, :a₁, :a₂, :a₃, :b₀, :b₁, :b₂, :b₃, :Φˣ),NTuple{15,TransformVariables.Identity}}},RdExpProductivity{RdExpData{Float64,Int64},IntegrationElements{Float64},Priors{Float64}}}},Float64},Float64,8},1}) at /Users/jack/.julia/packages/TransformVariables/eSbzF/src/generic.jl:166
 [3] logdensity(::TransformedLogDensity{TransformVariables.TransformTuple{NamedTuple{(:ρ_z, :σ_μ, :γᴵ, :γᴰ, :γᶠ, :γˢ, :a₀, :a₁, :a₂, :a₃, :b₀, :b₁, :b₂, :b₃, :Φˣ),NTuple{15,TransformVariables.Identity}}},RdExpProductivity{RdExpData{Float64,Int64},IntegrationElements{Float64},Priors{Float64}}}, ::Array{ForwardDiff.Dual{ForwardDiff.Tag{LogDensityProblems.var"#34#35"{TransformedLogDensity{TransformVariables.TransformTuple{NamedTuple{(:ρ_z, :σ_μ, :γᴵ, :γᴰ, :γᶠ, :γˢ, :a₀, :a₁, :a₂, :a₃, :b₀, :b₁, :b₂, :b₃, :Φˣ),NTuple{15,TransformVariables.Identity}}},RdExpProductivity{RdExpData{Float64,Int64},IntegrationElements{Float64},Priors{Float64}}}},Float64},Float64,8},1}) at /Users/jack/.julia/packages/LogDensityProblems/7AvD8/src/LogDensityProblems.jl:168
 [4] (::LogDensityProblems.var"#34#35"{TransformedLogDensity{TransformVariables.TransformTuple{NamedTuple{(:ρ_z, :σ_μ, :γᴵ, :γᴰ, :γᶠ, :γˢ, :a₀, :a₁, :a₂, :a₃, :b₀, :b₁, :b₂, :b₃, :Φˣ),NTuple{15,TransformVariables.Identity}}},RdExpProductivity{RdExpData{Float64,Int64},IntegrationElements{Float64},Priors{Float64}}}})(::Array{ForwardDiff.Dual{ForwardDiff.Tag{LogDensityProblems.var"#34#35"{TransformedLogDensity{TransformVariables.TransformTuple{NamedTuple{(:ρ_z, :σ_μ, :γᴵ, :γᴰ, :γᶠ, :γˢ, :a₀, :a₁, :a₂, :a₃, :b₀, :b₁, :b₂, :b₃, :Φˣ),NTuple{15,TransformVariables.Identity}}},RdExpProductivity{RdExpData{Float64,Int64},IntegrationElements{Float64},Priors{Float64}}}},Float64},Float64,8},1}) at /Users/jack/.julia/packages/LogDensityProblems/7AvD8/src/AD_ForwardDiff.jl:21
 [5] chunk_mode_gradient!(::DiffResults.MutableDiffResult{1,Float64,Tuple{Array{Float64,1}}}, ::LogDensityProblems.var"#34#35"{TransformedLogDensity{TransformVariables.TransformTuple{NamedTuple{(:ρ_z, :σ_μ, :γᴵ, :γᴰ, :γᶠ, :γˢ, :a₀, :a₁, :a₂, :a₃, :b₀, :b₁, :b₂, :b₃, :Φˣ),NTuple{15,TransformVariables.Identity}}},RdExpProductivity{RdExpData{Float64,Int64},IntegrationElements{Float64},Priors{Float64}}}}, ::Array{Float64,1}, ::ForwardDiff.GradientConfig{ForwardDiff.Tag{LogDensityProblems.var"#34#35"{TransformedLogDensity{TransformVariables.TransformTuple{NamedTuple{(:ρ_z, :σ_μ, :γᴵ, :γᴰ, :γᶠ, :γˢ, :a₀, :a₁, :a₂, :a₃, :b₀, :b₁, :b₂, :b₃, :Φˣ),NTuple{15,TransformVariables.Identity}}},RdExpProductivity{RdExpData{Float64,Int64},IntegrationElements{Float64},Priors{Float64}}}},Float64},Float64,8,Array{ForwardDiff.Dual{ForwardDiff.Tag{LogDensityProblems.var"#34#35"{TransformedLogDensity{TransformVariables.TransformTuple{NamedTuple{(:ρ_z, :σ_μ, :γᴵ, :γᴰ, :γᶠ, :γˢ, :a₀, :a₁, :a₂, :a₃, :b₀, :b₁, :b₂, :b₃, :Φˣ),NTuple{15,TransformVariables.Identity}}},RdExpProductivity{RdExpData{Float64,Int64},IntegrationElements{Float64},Priors{Float64}}}},Float64},Float64,8},1}}) at /Users/jack/.julia/packages/ForwardDiff/cXTw0/src/gradient.jl:140
 [6] gradient!(::DiffResults.MutableDiffResult{1,Float64,Tuple{Array{Float64,1}}}, ::LogDensityProblems.var"#34#35"{TransformedLogDensity{TransformVariables.TransformTuple{NamedTuple{(:ρ_z, :σ_μ, :γᴵ, :γᴰ, :γᶠ, :γˢ, :a₀, :a₁, :a₂, :a₃, :b₀, :b₁, :b₂, :b₃, :Φˣ),NTuple{15,TransformVariables.Identity}}},RdExpProductivity{RdExpData{Float64,Int64},IntegrationElements{Float64},Priors{Float64}}}}, ::Array{Float64,1}, ::ForwardDiff.GradientConfig{ForwardDiff.Tag{LogDensityProblems.var"#34#35"{TransformedLogDensity{TransformVariables.TransformTuple{NamedTuple{(:ρ_z, :σ_μ, :γᴵ, :γᴰ, :γᶠ, :γˢ, :a₀, :a₁, :a₂, :a₃, :b₀, :b₁, :b₂, :b₃, :Φˣ),NTuple{15,TransformVariables.Identity}}},RdExpProductivity{RdExpData{Float64,Int64},IntegrationElements{Float64},Priors{Float64}}}},Float64},Float64,8,Array{ForwardDiff.Dual{ForwardDiff.Tag{LogDensityProblems.var"#34#35"{TransformedLogDensity{TransformVariables.TransformTuple{NamedTuple{(:ρ_z, :σ_μ, :γᴵ, :γᴰ, :γᶠ, :γˢ, :a₀, :a₁, :a₂, :a₃, :b₀, :b₁, :b₂, :b₃, :Φˣ),NTuple{15,TransformVariables.Identity}}},RdExpProductivity{RdExpData{Float64,Int64},IntegrationElements{Float64},Priors{Float64}}}},Float64},Float64,8},1}}, ::Val{true}) at /Users/jack/.julia/packages/ForwardDiff/cXTw0/src/gradient.jl:37
 [7] gradient!(::DiffResults.MutableDiffResult{1,Float64,Tuple{Array{Float64,1}}}, ::LogDensityProblems.var"#34#35"{TransformedLogDensity{TransformVariables.TransformTuple{NamedTuple{(:ρ_z, :σ_μ, :γᴵ, :γᴰ, :γᶠ, :γˢ, :a₀, :a₁, :a₂, :a₃, :b₀, :b₁, :b₂, :b₃, :Φˣ),NTuple{15,TransformVariables.Identity}}},RdExpProductivity{RdExpData{Float64,Int64},IntegrationElements{Float64},Priors{Float64}}}}, ::Array{Float64,1}, ::ForwardDiff.GradientConfig{ForwardDiff.Tag{LogDensityProblems.var"#34#35"{TransformedLogDensity{TransformVariables.TransformTuple{NamedTuple{(:ρ_z, :σ_μ, :γᴵ, :γᴰ, :γᶠ, :γˢ, :a₀, :a₁, :a₂, :a₃, :b₀, :b₁, :b₂, :b₃, :Φˣ),NTuple{15,TransformVariables.Identity}}},RdExpProductivity{RdExpData{Float64,Int64},IntegrationElements{Float64},Priors{Float64}}}},Float64},Float64,8,Array{ForwardDiff.Dual{ForwardDiff.Tag{LogDensityProblems.var"#34#35"{TransformedLogDensity{TransformVariables.TransformTuple{NamedTuple{(:ρ_z, :σ_μ, :γᴵ, :γᴰ, :γᶠ, :γˢ, :a₀, :a₁, :a₂, :a₃, :b₀, :b₁, :b₂, :b₃, :Φˣ),NTuple{15,TransformVariables.Identity}}},RdExpProductivity{RdExpData{Float64,Int64},IntegrationElements{Float64},Priors{Float64}}}},Float64},Float64,8},1}}) at /Users/jack/.julia/packages/ForwardDiff/cXTw0/src/gradient.jl:33
 [8] logdensity_and_gradient(::LogDensityProblems.ForwardDiffLogDensity{TransformedLogDensity{TransformVariables.TransformTuple{NamedTuple{(:ρ_z, :σ_μ, :γᴵ, :γᴰ, :γᶠ, :γˢ, :a₀, :a₁, :a₂, :a₃, :b₀, :b₁, :b₂, :b₃, :Φˣ),NTuple{15,TransformVariables.Identity}}},RdExpProductivity{RdExpData{Float64,Int64},IntegrationElements{Float64},Priors{Float64}}},ForwardDiff.GradientConfig{ForwardDiff.Tag{LogDensityProblems.var"#34#35"{TransformedLogDensity{TransformVariables.TransformTuple{NamedTuple{(:ρ_z, :σ_μ, :γᴵ, :γᴰ, :γᶠ, :γˢ, :a₀, :a₁, :a₂, :a₃, :b₀, :b₁, :b₂, :b₃, :Φˣ),NTuple{15,TransformVariables.Identity}}},RdExpProductivity{RdExpData{Float64,Int64},IntegrationElements{Float64},Priors{Float64}}}},Float64},Float64,8,Array{ForwardDiff.Dual{ForwardDiff.Tag{LogDensityProblems.var"#34#35"{TransformedLogDensity{TransformVariables.TransformTuple{NamedTuple{(:ρ_z, :σ_μ, :γᴵ, :γᴰ, :γᶠ, :γˢ, :a₀, :a₁, :a₂, :a₃, :b₀, :b₁, :b₂, :b₃, :Φˣ),NTuple{15,TransformVariables.Identity}}},RdExpProductivity{RdExpData{Float64,Int64},IntegrationElements{Float64},Priors{Float64}}}},Float64},Float64,8},1}}}, ::Array{Float64,1}) at /Users/jack/.julia/packages/LogDensityProblems/7AvD8/src/AD_ForwardDiff.jl:51
 [9] evaluate_ℓ(::LogDensityProblems.ForwardDiffLogDensity{TransformedLogDensity{TransformVariables.TransformTuple{NamedTuple{(:ρ_z, :σ_μ, :γᴵ, :γᴰ, :γᶠ, :γˢ, :a₀, :a₁, :a₂, :a₃, :b₀, :b₁, :b₂, :b₃, :Φˣ),NTuple{15,TransformVariables.Identity}}},RdExpProductivity{RdExpData{Float64,Int64},IntegrationElements{Float64},Priors{Float64}}},ForwardDiff.GradientConfig{ForwardDiff.Tag{LogDensityProblems.var"#34#35"{TransformedLogDensity{TransformVariables.TransformTuple{NamedTuple{(:ρ_z, :σ_μ, :γᴵ, :γᴰ, :γᶠ, :γˢ, :a₀, :a₁, :a₂, :a₃, :b₀, :b₁, :b₂, :b₃, :Φˣ),NTuple{15,TransformVariables.Identity}}},RdExpProductivity{RdExpData{Float64,Int64},IntegrationElements{Float64},Priors{Float64}}}},Float64},Float64,8,Array{ForwardDiff.Dual{ForwardDiff.Tag{LogDensityProblems.var"#34#35"{TransformedLogDensity{TransformVariables.TransformTuple{NamedTuple{(:ρ_z, :σ_μ, :γᴵ, :γᴰ, :γᶠ, :γˢ, :a₀, :a₁, :a₂, :a₃, :b₀, :b₁, :b₂, :b₃, :Φˣ),NTuple{15,TransformVariables.Identity}}},RdExpProductivity{RdExpData{Float64,Int64},IntegrationElements{Float64},Priors{Float64}}}},Float64},Float64,8},1}}}, ::Array{Float64,1}) at /Users/jack/.julia/packages/DynamicHMC/10IZC/src/hamiltonian.jl:193
 [10] initialize_warmup_state(::Random._GLOBAL_RNG, ::LogDensityProblems.ForwardDiffLogDensity{TransformedLogDensity{TransformVariables.TransformTuple{NamedTuple{(:ρ_z, :σ_μ, :γᴵ, :γᴰ, :γᶠ, :γˢ, :a₀, :a₁, :a₂, :a₃, :b₀, :b₁, :b₂, :b₃, :Φˣ),NTuple{15,TransformVariables.Identity}}},RdExpProductivity{RdExpData{Float64,Int64},IntegrationElements{Float64},Priors{Float64}}},ForwardDiff.GradientConfig{ForwardDiff.Tag{LogDensityProblems.var"#34#35"{TransformedLogDensity{TransformVariables.TransformTuple{NamedTuple{(:ρ_z, :σ_μ, :γᴵ, :γᴰ, :γᶠ, :γˢ, :a₀, :a₁, :a₂, :a₃, :b₀, :b₁, :b₂, :b₃, :Φˣ),NTuple{15,TransformVariables.Identity}}},RdExpProductivity{RdExpData{Float64,Int64},IntegrationElements{Float64},Priors{Float64}}}},Float64},Float64,8,Array{ForwardDiff.Dual{ForwardDiff.Tag{LogDensityProblems.var"#34#35"{TransformedLogDensity{TransformVariables.TransformTuple{NamedTuple{(:ρ_z, :σ_μ, :γᴵ, :γᴰ, :γᶠ, :γˢ, :a₀, :a₁, :a₂, :a₃, :b₀, :b₁, :b₂, :b₃, :Φˣ),NTuple{15,TransformVariables.Identity}}},RdExpProductivity{RdExpData{Float64,Int64},IntegrationElements{Float64},Priors{Float64}}}},Float64},Float64,8},1}}}; q::Array{Float64,1}, κ::GaussianKineticEnergy{Diagonal{Float64,Array{Float64,1}},Diagonal{Float64,Array{Float64,1}}}, ϵ::Nothing) at /Users/jack/.julia/packages/DynamicHMC/10IZC/src/mcmc.jl:109
 [11] initialize_warmup_state(::Random._GLOBAL_RNG, ::LogDensityProblems.ForwardDiffLogDensity{TransformedLogDensity{TransformVariables.TransformTuple{NamedTuple{(:ρ_z, :σ_μ, :γᴵ, :γᴰ, :γᶠ, :γˢ, :a₀, :a₁, :a₂, :a₃, :b₀, :b₁, :b₂, :b₃, :Φˣ),NTuple{15,TransformVariables.Identity}}},RdExpProductivity{RdExpData{Float64,Int64},IntegrationElements{Float64},Priors{Float64}}},ForwardDiff.GradientConfig{ForwardDiff.Tag{LogDensityProblems.var"#34#35"{TransformedLogDensity{TransformVariables.TransformTuple{NamedTuple{(:ρ_z, :σ_μ, :γᴵ, :γᴰ, :γᶠ, :γˢ, :a₀, :a₁, :a₂, :a₃, :b₀, :b₁, :b₂, :b₃, :Φˣ),NTuple{15,TransformVariables.Identity}}},RdExpProductivity{RdExpData{Float64,Int64},IntegrationElements{Float64},Priors{Float64}}}},Float64},Float64,8,Array{ForwardDiff.Dual{ForwardDiff.Tag{LogDensityProblems.var"#34#35"{TransformedLogDensity{TransformVariables.TransformTuple{NamedTuple{(:ρ_z, :σ_μ, :γᴵ, :γᴰ, :γᶠ, :γˢ, :a₀, :a₁, :a₂, :a₃, :b₀, :b₁, :b₂, :b₃, :Φˣ),NTuple{15,TransformVariables.Identity}}},RdExpProductivity{RdExpData{Float64,Int64},IntegrationElements{Float64},Priors{Float64}}}},Float64},Float64,8},1}}}) at /Users/jack/.julia/packages/DynamicHMC/10IZC/src/mcmc.jl:109
 [12] mcmc_keep_warmup(::Random._GLOBAL_RNG, ::LogDensityProblems.ForwardDiffLogDensity{TransformedLogDensity{TransformVariables.TransformTuple{NamedTuple{(:ρ_z, :σ_μ, :γᴵ, :γᴰ, :γᶠ, :γˢ, :a₀, :a₁, :a₂, :a₃, :b₀, :b₁, :b₂, :b₃, :Φˣ),NTuple{15,TransformVariables.Identity}}},RdExpProductivity{RdExpData{Float64,Int64},IntegrationElements{Float64},Priors{Float64}}},ForwardDiff.GradientConfig{ForwardDiff.Tag{LogDensityProblems.var"#34#35"{TransformedLogDensity{TransformVariables.TransformTuple{NamedTuple{(:ρ_z, :σ_μ, :γᴵ, :γᴰ, :γᶠ, :γˢ, :a₀, :a₁, :a₂, :a₃, :b₀, :b₁, :b₂, :b₃, :Φˣ),NTuple{15,TransformVariables.Identity}}},RdExpProductivity{RdExpData{Float64,Int64},IntegrationElements{Float64},Priors{Float64}}}},Float64},Float64,8,Array{ForwardDiff.Dual{ForwardDiff.Tag{LogDensityProblems.var"#34#35"{TransformedLogDensity{TransformVariables.TransformTuple{NamedTuple{(:ρ_z, :σ_μ, :γᴵ, :γᴰ, :γᶠ, :γˢ, :a₀, :a₁, :a₂, :a₃, :b₀, :b₁, :b₂, :b₃, :Φˣ),NTuple{15,TransformVariables.Identity}}},RdExpProductivity{RdExpData{Float64,Int64},IntegrationElements{Float64},Priors{Float64}}}},Float64},Float64,8},1}}}, ::Int64; initialization::Tuple{}, warmup_stages::Tuple{FindLocalOptimum{Float64},InitialStepsizeSearch,TuningNUTS{Nothing,DualAveraging{Float64}},TuningNUTS{Diagonal,DualAveraging{Float64}},TuningNUTS{Diagonal,DualAveraging{Float64}},TuningNUTS{Diagonal,DualAveraging{Float64}},TuningNUTS{Diagonal,DualAveraging{Float64}},TuningNUTS{Diagonal,DualAveraging{Float64}},TuningNUTS{Nothing,DualAveraging{Float64}}}, algorithm::DynamicHMC.NUTS{Val{:generalized}}, reporter::LogProgressReport{Nothing}) at /Users/jack/.julia/packages/DynamicHMC/10IZC/src/mcmc.jl:518
 [13] macro expansion at /Users/jack/.julia/packages/UnPack/7uwKF/src/UnPack.jl:100 [inlined]
 [14] mcmc_with_warmup(::Random._GLOBAL_RNG, ::LogDensityProblems.ForwardDiffLogDensity{TransformedLogDensity{TransformVariables.TransformTuple{NamedTuple{(:ρ_z, :σ_μ, :γᴵ, :γᴰ, :γᶠ, :γˢ, :a₀, :a₁, :a₂, :a₃, :b₀, :b₁, :b₂, :b₃, :Φˣ),NTuple{15,TransformVariables.Identity}}},RdExpProductivity{RdExpData{Float64,Int64},IntegrationElements{Float64},Priors{Float64}}},ForwardDiff.GradientConfig{ForwardDiff.Tag{LogDensityProblems.var"#34#35"{TransformedLogDensity{TransformVariables.TransformTuple{NamedTuple{(:ρ_z, :σ_μ, :γᴵ, :γᴰ, :γᶠ, :γˢ, :a₀, :a₁, :a₂, :a₃, :b₀, :b₁, :b₂, :b₃, :Φˣ),NTuple{15,TransformVariables.Identity}}},RdExpProductivity{RdExpData{Float64,Int64},IntegrationElements{Float64},Priors{Float64}}}},Float64},Float64,8,Array{ForwardDiff.Dual{ForwardDiff.Tag{LogDensityProblems.var"#34#35"{TransformedLogDensity{TransformVariables.TransformTuple{NamedTuple{(:ρ_z, :σ_μ, :γᴵ, :γᴰ, :γᶠ, :γˢ, :a₀, :a₁, :a₂, :a₃, :b₀, :b₁, :b₂, :b₃, :Φˣ),NTuple{15,TransformVariables.Identity}}},RdExpProductivity{RdExpData{Float64,Int64},IntegrationElements{Float64},Priors{Float64}}}},Float64},Float64,8},1}}}, ::Int64; initialization::Tuple{}, warmup_stages::Tuple{FindLocalOptimum{Float64},InitialStepsizeSearch,TuningNUTS{Nothing,DualAveraging{Float64}},TuningNUTS{Diagonal,DualAveraging{Float64}},TuningNUTS{Diagonal,DualAveraging{Float64}},TuningNUTS{Diagonal,DualAveraging{Float64}},TuningNUTS{Diagonal,DualAveraging{Float64}},TuningNUTS{Diagonal,DualAveraging{Float64}},TuningNUTS{Nothing,DualAveraging{Float64}}}, algorithm::DynamicHMC.NUTS{Val{:generalized}}, reporter::LogProgressReport{Nothing}) at /Users/jack/.julia/packages/DynamicHMC/10IZC/src/mcmc.jl:577
 [15] mcmc_with_warmup(::Random._GLOBAL_RNG, ::LogDensityProblems.ForwardDiffLogDensity{TransformedLogDensity{TransformVariables.TransformTuple{NamedTuple{(:ρ_z, :σ_μ, :γᴵ, :γᴰ, :γᶠ, :γˢ, :a₀, :a₁, :a₂, :a₃, :b₀, :b₁, :b₂, :b₃, :Φˣ),NTuple{15,TransformVariables.Identity}}},RdExpProductivity{RdExpData{Float64,Int64},IntegrationElements{Float64},Priors{Float64}}},ForwardDiff.GradientConfig{ForwardDiff.Tag{LogDensityProblems.var"#34#35"{TransformedLogDensity{TransformVariables.TransformTuple{NamedTuple{(:ρ_z, :σ_μ, :γᴵ, :γᴰ, :γᶠ, :γˢ, :a₀, :a₁, :a₂, :a₃, :b₀, :b₁, :b₂, :b₃, :Φˣ),NTuple{15,TransformVariables.Identity}}},RdExpProductivity{RdExpData{Float64,Int64},IntegrationElements{Float64},Priors{Float64}}}},Float64},Float64,8,Array{ForwardDiff.Dual{ForwardDiff.Tag{LogDensityProblems.var"#34#35"{TransformedLogDensity{TransformVariables.TransformTuple{NamedTuple{(:ρ_z, :σ_μ, :γᴵ, :γᴰ, :γᶠ, :γˢ, :a₀, :a₁, :a₂, :a₃, :b₀, :b₁, :b₂, :b₃, :Φˣ),NTuple{15,TransformVariables.Identity}}},RdExpProductivity{RdExpData{Float64,Int64},IntegrationElements{Float64},Priors{Float64}}}},Float64},Float64,8},1}}}, ::Int64) at /Users/jack/.julia/packages/DynamicHMC/10IZC/src/mcmc.jl:577
 [16] top-level scope at none:0

Can you please help me look at the error and give me some suggestions?

@y1my1
Copy link
Author

y1my1 commented Apr 18, 2020

Hi @tpapp , I just wonder if you would mind taking a look at this when you have time.
To better use your time, I have created an MWE which basically follows the linear regression example here in this repo

This following one runs smothly

using TransformVariables, LogDensityProblems, DynamicHMC, DynamicHMC.Diagnostics
using MCMCDiagnostics
using Parameters, Statistics, Random, Distributions, StatsFuns
import ForwardDiff

struct LinearRegProblem{Ty<:AbstractVector,TX<:AbstractMatrix,Tv<:Real}
    y::Ty
    X::TX
    v::Tv
end

function (problem::LinearRegProblem)(θ)
    @unpack y, X, v = problem
    @unpack β₀, β₁, σ = θ
    @info β₀, β₁, σ

    return LogLikelihood(y, X, β₀, β₁, σ) + LogPriorProbability(v, σ)
end

function LogLikelihood(y, X, β₀, β₁, σ)
    return loglikelihood(Normal(0,σ), y.-X*[β₀, β₁])
end

function LogPriorProbability(v, σ)
    return logpdf(TDist(v), σ)
end

N = 2000
X = hcat(ones(N), randn(N,1))
β₀ = 1.8
β₁ = 2.5
σ = .5
y = X*[β₀, β₁] .+ randn(N).*σ

p = LinearRegProblem(y, X, 1.0)
p((β₀=β₀, β₁=β₁, σ=σ))

trans = as((β₀=asℝ, β₁=asℝ, σ=asℝ₊))
tP = TransformedLogDensity(trans, p)
∇P = ADgradient(:ForwardDiff, tP);

results = mcmc_with_warmup(Random.GLOBAL_RNG, ∇P, 1000)

posterior = transform.(trans, results.chain)

post_β₁ = mean(first, posterior)
post_σ = mean(last, posterior)

This one will give the error

using TransformVariables, LogDensityProblems, DynamicHMC, DynamicHMC.Diagnostics
using MCMCDiagnostics
using Parameters, Statistics, Random, Distributions, StatsFuns
import ForwardDiff

@with_kw struct DynamicParameters{F<:AbstractFloat} @deftype F
    β₀
    β₁
    σ
end

struct LinearRegProblem{Ty<:AbstractVector,TX<:AbstractMatrix,Tv<:Real}
    y::Ty
    X::TX
    v::Tv
end

function (problem::LinearRegProblem)(θ)
    @unpack y, X, v = problem
    @unpack β₀, β₁, σ = θ
    @info β₀, β₁, σ

    Θ = DynamicParameters(β₀=β₀, β₁=β₁, σ=σ)
    @info Θ
    return LogLikelihood(y, X, Θ) + LogPriorProbability(v, Θ)
end

function LogLikelihood(y, X, Θ::DynamicParameters)
    @unpack β₀, β₁, σ = Θ
    return loglikelihood(Normal(0,σ), y.-X*[β₀, β₁])
end

function LogPriorProbability(v, Θ::DynamicParameters)
    @unpack σ = Θ
    return logpdf(TDist(v), σ)
end

N = 2000
X = hcat(ones(N), randn(N,1))
β₀ = 1.8
β₁ = 2.5
σ = .5
y = X*[β₀, β₁] .+ randn(N).*σ

p = LinearRegProblem(y, X, 1.0)
p((β₀=β₀, β₁=β₁, σ=σ))

trans = as((β₀=asℝ, β₁=asℝ, σ=asℝ₊))
tP = TransformedLogDensity(trans, p)
∇P = ADgradient(:ForwardDiff, tP);

results = mcmc_with_warmup(Random.GLOBAL_RNG, ∇P, 1000)

The difference is that in the second one I wrapped up β₀, β₁, σ into a struct called DynamicParameters and changed the relevant functions. The reason I want to do this is to make my functions somewhat cleaner by passing a struct in my real application. Please forgive me if this looks silly.

Thanks a lot for your help,
Best regards,
Jack

@tpapp
Copy link
Owner

tpapp commented Apr 18, 2020

There are multiple things going on here:

  1. as for code organization, I would generally define the posterior in a single callable. If that is not a good option for you (eg you reuse the same prior for multiple models), then there is no harm passing around the θ as a NamedTuple to both the prior and the posterior, and just have them extract what they need. So I would @unpack where I need the parameters, and not before.

  2. I would not define a struct for parameters, but you can spline in a NamedTuple for keywords, using the f(; θ...) syntax. That said, it may complain about unused parameters. So instead of using the @with_kw macro, you could define your own constructor that takes keywords and just tag a _... at the end to collect unused ones.

  3. If you are using ForwardDiff, everything needs to work with Dual types, too. See http://www.juliadiff.org/ForwardDiff.jl/latest/dev/how_it_works.html .

Most of these issues are not specific to DynamicHMC BTW, just general Julia code organization questions (I am happy to help though). Are you just learning Julia?

@y1my1
Copy link
Author

y1my1 commented Apr 18, 2020

Thanks a lot for your suggestions. That's extremely helpful. I will follow your advice to solve the problem.

To be honest, I just started programming in Julia for a few weeks. Not to seek for excuses, as you know for a large portion of economics people, I used to program in Stata (if that's programming) and some Python, but I like Julia a lot and plan to use it as my main programming language looking forward.

I thought it may be a good way to be familiar with Julia by solving my estimation problem through learning by doing. And I mostly learn from others' programs and search google when encountering some new to me. I guess I need to study more thoroughly the Julia manual. But maybe the Julia manual is a little bit too detailed for the beginner. Do you have something to recommend to enhance my Julia programming skills? Thanks.

Best regards,

@tpapp
Copy link
Owner

tpapp commented Apr 18, 2020

I think you are actually doing fine; as you have figured out a lot of things. Just give it time.

IMO the best source is the manual, and reading code from others; eg hanging around on Discourse and reading/contributing.

@y1my1
Copy link
Author

y1my1 commented Apr 18, 2020

Great. I found Discourse a great place for constructive discussions but I generally read.

Thanks. Have a great weekend.

@y1my1 y1my1 closed this as completed Apr 18, 2020
@y1my1
Copy link
Author

y1my1 commented Apr 18, 2020

By the way, I also want to contribute someway, please feel free to let me know if there is something I could help here related to DynamicHMC or somewhere else generally in Julia community.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

2 participants