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.

1 Like
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
11

Yeah, just use the whole posterior, in general the posterior could be any shape and MCMC is run to infer it without too many assumptions. In any case, pairwise scatterplots are a handy check to go with the check on the mixing and convergence of chains.

using PairPlots, CairoMakie
pairplot(chain)

image
image1200×1276 76 KB

1 Like
12

The method I shared samples directly from the marginals. So yeah in that case the correlation is lost. I think in your case you want to sample directly from the chain. But maybe @penelopeysm or someone else from turing can better help with that (there is also a turing channel on slack)

1 Like