What is the canon way to randomly draw one (or more) parameter sets from a sampled distribution?

Specific Domains Probabilistic Programming
1

Using turing I have done

n_steps = 1000
n_chains = 4
chain = sample(model, NUTS(), MCMCThreads(), n_steps, n_chains)

I can plot the output using

plot(chain)

image
image532×398 94.7 KB

However, if I want to draw a random parameter set to do something with, is there a canon way to do this? I.e. here I have done inference for parameters γ, ν, σI and want to draw a single value from that joint distribution. I practise, I think I could use

n_ps = 3
collect(chain.value[rand(1:n_steps), :, rand(1:n_chains)])[1:n_ps]

however, I also suspect that there might exist a more direct interface.

2

With the Array constructor you can at least avoid accessing the value field, e.g. Array(chain, append_chains=false)[rand(1:n_chains)][rand(1:n_steps), :] (Getting started · MCMCChains.jl). If the chains can be pooled together, it is Array(chain)[rand(1:n_steps), :].

1 Like
3

you can do:

sample(chain[:γ], 10)
1 Like
4

Thanks, works great! If I want to sample all three parameters, is there some version of

sample(chain[:γ, :ν, :σI], 1)[1] # errors

that I can do (which is a bit nicer than

[sample(chain[:γ], 1)[1], sample(chain[:ν], 1)[1] , sample(chain[:σI], 1)[1] ]

?

5

you could do sample.((chain[:γ], chain[:ν]), 10)

1 Like
6

That works, thanks a lot!

7

Just checking, when I do this, the different parameter values are correctly correlated, i.e. I am not sampling γ and ν values independently of each other, right?

8

You didn’t share your model, but I don’t see how γ and ν can be correlated?

1 Like
9

For reference, this is the full workflow:

using Catalyst
sir = @reaction_network begin
    γ, S + I --> 2I
    ν, I --> R
end

using OrdinaryDiffEqDefault, Plots
u0 = [:S => 999.0, :I => 1.0, :R => 0.0]
p_true = [:γ => 0.0005, :ν => 0.1]
oprob_true = ODEProblem(sir, u0, 100.0, p_true)
sol_true = solve(oprob_true)
plot(sol_true; label = ["S (true)" "I (true)" "R (true)"], lw = 4)

using Distributions
σI = 10.0
t_measurement = 1.0:2:100.0
I_observed = sol_true(t_measurement; idxs = :I)
I_observed = [rand(Normal(I, σI)) for I in I_observed]
plot!(t_measurement, I_observed; label = "I (measured)", color = 2)

using Turing, SymbolicIndexingInterface
setp_oop = SymbolicIndexingInterface.setp_oop(oprob_true, [:γ, :ν])
@model function sir_likelihood(I_observed, oprob_base, setp_oop, saveat)
    # Defines the parameters we wish to estimate and their prior distributions.
    γ ~ LogUniform(0.00001, 0.001)
    ν ~ LogUniform(0.01, 1.0)
    σI ~ LogUniform(0.1, 100.0)

    # Simulate the model for parameter values γ, ν. Saves the solution at the measurement times.
    p = setp_oop(oprob_base, [γ, ν])
    oprob_base = remake(oprob_base; p)
    sol = solve(oprob_base; saveat, verbose = false, maxiters = 10000)

    # If simulation was unsuccessful, the likelihood is -Inf.
    if !SciMLBase.successful_retcode(sol)
        Turing.@addlogprob! -Inf
        return nothing
    end

    # Computes the likelihood of the observations.
    for idx in eachindex(I_observed)
        I_observed[idx] ~ Normal(sol[:I][idx], σI)
    end
end

n_steps = 1000
n_chains = 4
sir_model = sir_likelihood(I_observed, oprob_true, setp_oop, t_measurement)
chain = sample(sir_model, NUTS(), MCMCThreads(), n_steps, n_chains; progress = false)

The result should be a 3d posterior over the estimated parameters γ, ν, and σI, respectively(?). Imagine there is a correlation between γ and ν in this distribution. In this case, wouldn’t sampling a value of γ and one of ν independently yield a different result from sampling a single point in my posterior distribution and extracting its γ and ν values? I.e. there there are two different interpretations of the sampled point

first.(sample.((chain[:γ], chain[:ν]), 1))

where I am uncertain which of them it is.

1 Like