stan_init does not work if written in Stanmodel function

[aakhmetz]

When I write stan_init as a part of Stanmodel function, it looks like it is finally not assigned:

days_inc = vcat(fill.(df_inc[!,:day],df_inc[!,:n])...)
stan_data = Dict("J" => length(days_inc), "onset_day_reported" => days_inc)
stan_init = Dict("mu" => 2.5, "sigma" => 1.0, "incubation" => days_inc .+ 0.5)

stan_model = Stanmodel(
    model = stan_code, 
    init = [stan_init]);

_, stan_chns, _ = stan(stan_model, stan_data, summary = false);

As the log-file shows later:

data
  file = noname_3.data.R
init = 2 (Default)

Writing something like w/o brackets:

stan_model = Stanmodel(
    model = stan_code, 
    init = stan_init);

resolves in error:

MethodError: Cannot `convert` an object of type Dict{String,Any} to an object of type Array{Dict{String,Any},1}
Closest candidates are:

The working example is

_, stan_chns, _ = stan(stan_model, stan_data, init=stan_init, summary = false);

but as I understood from the manual the former variant should also work.

[nb] however, I tested this in Julia 1.3 [nb2] probably related to #28

[itsdfish]

Hi @aakhmetz. Without a MWE, my best guess is that your entry called incubation in the stan_init variable needs to be a Float64 rather than an array. You might be able to do something like this to populate your dictionary:

days_inc = rand(10)
d = Dict("incubation[$i]"=>v for (i,v) in enumerate(days_inc))

Result:

Dict{String,Float64} with 10 entries:
  "incubation[4]" => 0.526845
  "incubation[6]" => 0.868194
  "incubation[2]" => 0.645691
  "incubation[5]" => 0.972755
  "incubation[8]" => 0.52399
  "incubation[7]" => 0.172707
  "incubation[9]" => 0.557837
  ⋮               => ⋮

[aakhmetz]

@itsdfish Thank you for checking. I tried but no success with a new version:

MethodError: Cannot `convert` an object of type Dict{String,Float64} to an object of type Array{Dict{String,Any},1}
Closest candidates are:
  convert(::Type{Array{T,N}}, !Matched::FillArrays.Zeros{V,N,Axes} where Axes) where {T, V, N} at /home/aakhmetz/.julia/packages/FillArrays/LYCik/src/FillArrays.jl:320
  convert(::Type{Array{T,N}}, !Matched::FillArrays.Ones{V,N,Axes} where Axes) where {T, V, N} at /home/aakhmetz/.julia/packages/FillArrays/LYCik/src/FillArrays.jl:320
  convert(::Type{Array{T,N}}, !Matched::FillArrays.AbstractFill{V,N,Axes} where Axes) where {T, V, N} at /home/aakhmetz/.julia/packages/FillArrays/LYCik/src/FillArrays.jl:312

My thought was that init should be equally treated as data.

If you have some time, could you check the gist? https://gist.github.com/aakhmetz/a2a754a2d6f31209010cc7bd8eeeb64a It may be the MWE

[itsdfish]

Thanks for the MWE. I had to add the other parameters to your init_stan dictionary. I was able to replicate your issue: the init worked when called from stan(), but not when set in the Stanmodel object. My guess is that there is a bug or the feature was removed and the docs were not updated. Rob would have a better idea.

Here is a working version version:

Note that Float64() is not necessary and I merged the two dictionaries so that it works.

using CmdStan, Test, Statistics

tempdir = pwd()

nwarmup, nsamples, nchains = 1000, 1000, 4

stan_code = "
data {
    int<lower = 1> J; // number of cases
    real<lower = 0> onset_day_reported[J]; //day of illness onset
}

parameters {
    real<lower = 0> incubation[J];

    real mu;
    real<lower = 0> sigma;
}

model {
    mu ~ normal(0, 5);
    sigma ~ cauchy(0, 5);

    incubation ~ lognormal(mu,sigma);

    for (j in 1:J) {
        target += normal_lpdf(incubation[j] | onset_day_reported[j]+0.5, 0.5);
    }

}

generated quantities {
    real incubation_mean = exp(mu+sigma^2/2.0);
    real incubation_sd = sqrt((exp(sigma^2)-1.0)*exp(2.0*mu+sigma^2));
}
"

df_inc = DataFrame(Dict("n" => [1, 5, 7, 9, 17, 18, 16, 6, 8, 10, 7, 1, 2, 4, 4, 1],
    "day" => [6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 25]))
days_inc = vcat(fill.(df_inc[!,:day],df_inc[!,:n])...)

stan_data = Dict("J" => length(days_inc), "onset_day_reported" => days_inc)

# stan_init = Dict("mu" => 2.5, "sigma" => 0.2, "incubation" => days_inc .+ 0.5)
stan_init = Dict("incubation[$i]"=>v+.5 for (i,v) in enumerate(days_inc))
temp_init = Dict("mu" => 2.5, "sigma" => 1.0)
merge!(stan_init, temp_init)

stan_model = Stanmodel(
    model = stan_code, 
    nchains = nchains,
    num_warmup = nwarmup,
    num_samples = nsamples);

_, stan_chns, _ = stan(stan_model, stan_data, summary = false, init = stan_init);

[aakhmetz]

@itsdfish Thank you! That explains my confusion, because I could not understand the error message coming from julia when I used init writting in the Stanmodel. Good!

[goedman]

Thanks Andrei for bringing this up. And thanks (again) Chris for jumping in. This works (replacing the latter part of Chris' MWE version:

stan_model = Stanmodel(
    model = stan_code,
    nchains = nchains,
    init = [stan_init],
    num_warmup = nwarmup,
    num_samples = nsamples);

_, stan_chns, _ = stan(stan_model, stan_data, summary = false);

chns = set_section(stan_chns, Dict(
  :parameters => ["mu", "sigma", "incubation_mean", "incubation_sd"],
  :incubations => ["incubation.$i" for i in 1:116],
  :internals => ["lp__", "accept_stat__", "stepsize__", "treedepth__", "n_leapfrog__",
    "divergent__", "energy__"]
  )
)

show(chns)

The bug is that I never updated init in Stanmodel to accept a single Dict, it expects an DataDict[]:

?Stanmodel

...
  * `monitors::String[] `        : Variables saved for post-processing
  * `data::DataDict[]`           : Observed input data
  * `random::Random`             : Random seed settings
  * `init::DataDict[]`           : Initial values for parameters
  * `output::Output`             : File output options
  * `printsummary=true`          : Show computed stan summary
...

I could turn DataDict into Union{ Dict{Symbol, Any}, Vector{Dict{Symbol, Any}} } and then handle it in stancode.jl. Vaguely I remember I used to have a warning for this issue in the docs (that it had to be a Vector{Dict} in Stanmodel). I later on preferred data and init to be passed in to stan(), what do you think?

[itsdfish]

No problem.

If you always want to process a vector, another option is to make an alternative method and pass it to the current method. For example,

my_fun(vector{Float64}::x) = ...
my_fun(Float64::x) = my_fun([x])

[goedman]

Great suggestion. Most folks provide a single Dict, so that is kind of the default.

[goedman]

Unfortunately, as I probably also found out several years ago, given that these arguments are all keyword arguments they don't participate in dispatch. I'm sure there is a way to handle this, but it is not trivial.

I'll release a version with the note restored in the on-line documentation about setting up init in Stanmodel().

I actually would prefer to remove this option altogether, but that is a breaking change and might not be worth it.

[itsdfish]

Ah. I understand. If you do not want to make a breaking change, you can dispatch on keywords with a simple hack:

my_fun(;x) = my_fun(x)
my_fun(x::Vector{Float64}) = println(typeof(x))
my_fun(x::Float64) = my_fun([x])

Result:

julia> my_fun(x=.3)
Array{Float64,1}

[goedman]

Hi Chris,

With multiple keyword arguments I have never succeeded.

Something like below might work. This is one of my many mental blocks in Julia. I mean the issue of contravariance.

datatype = Union{
  Dict{String, T} where {T <: Any},
  Vector{Dict{String, T}} where {T <: Any}
}

mutable struct Mymodel
  data::datatype
  init::datatype
end

function Mymodel(;data=Nothing, init=Nothing)

  # Converts required if not defined????
  mydata = data == Nothing ? Dict{String, Any}() : data
  myinit = init == Nothing ? Dict{String, Any}() : init

  Mymodel(mydata, myinit)
end

dd = Dict("theta" => 1, "alpha" => 3.0)
id = Dict("sigma" => 1.0)

m1 = Mymodel(data=dd, init=id)

@show m1
println()

dd2 = [Dict("theta" => 1, "alpha" => 3.0),
  Dict("theta" => 1, "alpha" => 3.0)]

id2 = [Dict("sigma" => 1.0)]

m2 = Mymodel(data=dd2, init=id2)

@show m2
println()

m3 = Mymodel(data=dd2)

@show m3

But I can't figure out why I can't just create a Mymodel with Nothings in there (this is what stan() expects. I tried extending the datatype Union, but clearly not in the right way.

[itsdfish]

Hmmm. One minor change you could make is:

function Mymodel(;data= Dict{String, Any}(), init= Dict{String, Any}())
  Mymodel(data, init)
end

Can you elaborate on creating an instance of Mymodel with Nothing?

[goedman]

Yes, you are right. That was leftover of including Nothing in the datadict Union.

Anyway, this feature was confusing and not consistently implemented (didn't work for data in Stanmodel() and in case of init it took precedence over init data passed in using stan().

So I decided to remove both fields from the Stanmodel.

Updated the docs, README (with credits to both of you) and will release it as v5.4.0 as soon as Travis is done.

Thanks for your help!

[goedman]

v5.4.0 has been merged.

[aakhmetz]

Thank you, @itsdfish and @goedman, for your feedback, clarifying the situation, and getting it fixed! I think the issue can be closed now. I presume I will need to update my StanJulia though.

[goedman]

Your welcome. Just Pkg.update("CmdStan") should do the trick. Or Pkg.update() to update all your packages.