Easiest way to programatically set initial parameters based on the mean of the priors

[DominiqueMakowski]

This is somewhat of a follow-up on #2286

I would like to set my initial parameters to the mean (or the mode) of the prior distributions.

using Turing, StatsFuns, DataFrames, Random


@model function mod(y, x1, x2)

    μ_intercept ~ Normal(0, 0.5)
    μ_x1 ~ Normal(1, 0.5)
    μ_x2 ~ Normal(2, 0.5)
    σ ~ truncated(Normal(0.0, 1), lower=0)

    μ = μ_intercept .+ μ_x1 .* x1 .+ μ_x2 .* x2
    y .~ Normal.(μ, σ)
end

The most optimal way in theory, as far as I can see, would be to have a prior extractor function that extracts the distributions used from a model. This would be wonderful for many other applications, including visualization of models, but in my case it would allow to then compute the initial parameters analytically:

# Ideal solution: extract priors as a collection of distributions
priors = (
    μ_intercept=Normal(0, 0.5),
    μ_x1=Normal(1, 0.5),
    μ_x2=Normal(2, 0.5),
    σ=truncated(Normal(0.0, 1), lower=0)
)  # Assuming this was obtained with something like extract_priors()
initial_parameters = mean.(collect(priors))

In the issue mentioned above, the alternative proposed was to first sample from the priors, and then compute the indices of interest. However, two approaches were suggested.

1. Using rand()

# Existing solution
fit = mod(rand(10), rand(10), rand(10))
priors_samples = hcat([rand(Vector, fit) for _ in 1:100]...)
initial_parameters = mean(priors_samples, dims=2)

2. Using sample()

# Alternative
fit = mod(rand(10), rand(10), rand(10))
chains = sample(fit, Prior(), 100)
initial_parameters = mean(Array(chains), dims=1)

Both of these solutions have some overhead as one must fit the model first. Additionally, using sample() feels slower than rand() (?)

My question is, what would be the fastest way and the most "Turingy" approach? Thanks

[torfjelde]

We actually have a method for extracting the priors: extract_priors:)

So I think the following should do the trick:

mapreduce(DynamicPPL.tovec ∘ mean, vcat, values(extract_priors(model)))

This does the following:

1. Extract the priors (into an OrderedDict)
2. Compute the mean.
3. Convert that to a vec using the tovec internal method.
4. Reduce using vcat

[DominiqueMakowski]

Once again I spoke too fast, the number of parameters that I get is not what the model is expecting... will dig further

[torfjelde]

What version of Turing are you on? The tovec approach works on my end:) I'm on v0.35.2

[DominiqueMakowski]

julia>         mapreduce(DynamicPPL.tovec ∘ mean, vcat, values(extract_priors(x)))
ERROR: UndefVarError: `tovec` not defined in `DynamicPPL`
Suggestion: check for spelling errors or missing imports.
Stacktrace:
 [1] getproperty(x::Module, f::Symbol)
   @ Base .\Base.jl:42
 [2] top-level scope
   @ c:\Users\domma\Dropbox\RECHERCHE\Studies\DoggoNogoValidation\pilot\analysis\2_models_sample.jl:44

(analysis) pkg> status
Status `C:\Users\domma\Dropbox\RECHERCHE\Studies\DoggoNogoValidation\pilot\analysis\Project.toml`
  [336ed68f] CSV v0.10.15
⌃ [13f3f980] CairoMakie v0.12.15
  [a93c6f00] DataFrames v1.7.0
⌅ [366bfd00] DynamicPPL v0.28.6
⌃ [da2b9cff] Mooncake v0.4.1
  [a68b5a21] ParetoSmooth v0.7.11
⌃ [0e71a2a6] SequentialSamplingModels v0.11.12
  [4c63d2b9] StatsFuns v1.3.2
⌃ [fce5fe82] Turing v0.34.1

[torfjelde]

Ah yes, so you have to bump Turing to 0.35 (which should also update DynamicPPL.jl to the most recent version). Then this should work 👍

[DominiqueMakowski]

works! cheers