Problems with deterministic distribution

[mateuszbaran]

I'm trying to use deterministic distribution in my model as outlined here: https://discourse.julialang.org/t/is-there-a-turing-alternative-to-pm-deterministic-from-pymc3/38667/2 and so far my code fails:

Distributions.rand(rng::AbstractRNG, d::Determin) = d.val
Distributions.logpdf(d::Determin, x::T) where T<:Real = zero(x)

Bijectors.bijector(d::Determin) = Identity{0}() # <= `0` indicates 0-dimensional, i.e. univariate

@model model1(a, ::Type{T} = Float64) where {T} = begin
    ma ~ Normal(0.0, 0.5)
    for i in 1:length(a)
        a[i] ~ Normal(ma, 0.05)
    end

    b ~ Determin(-exp(-2*ma))
    return ma, b
end

function test_model1()
    # sampler = HMC(0.001, 10)
    # sampler = HMCDA(1000, 0.75, 0.05)
    sampler = NUTS(0.65)
    c1 = sample(model1(a1), sampler, 2000)

    println(c1)
end

Typical results look like this:

julia> test_model1()
┌ Info: Found initial step size
└   ϵ = 0.025
Object of type Chains, with data of type 1000×14×1 Array{Float64,3}

Iterations        = 1:1000
Thinning interval = 1
Chains            = 1
Samples per chain = 1000
internals         = acceptance_rate, hamiltonian_energy, hamiltonian_energy_error, is_accept, log_density, lp, max_hamiltonian_energy_error, n_steps, nom_step_size, numerical_error, step_size, tree_depth
parameters        = b, ma

2-element Array{ChainDataFrame,1}

Summary Statistics
  parameters                       mean                       std                 naive_se                      mcse       ess   r_hat
  ──────────  ─────────────────────────  ────────────────────────  ───────────────────────  ────────────────────────  ────────  ──────
           b  13816403112283807744.0000  7133491167269512192.0000  225580797572648256.0000  2172001283591027712.0000    4.8969  1.3982
          ma                     0.0690                    0.0179                   0.0006                    0.0007  830.2122  0.9998

Quantiles
  parameters                       2.5%                      25.0%                      50.0%                      75.0%                      97.5%
  ──────────  ─────────────────────────  ─────────────────────────  ─────────────────────────  ─────────────────────────  ─────────────────────────
           b  -5026632495806943232.0000  12090906044107700224.0000  15653120955249981440.0000  18317459011682740224.0000  22252367610115112960.0000
          ma                     0.0349                     0.0565                     0.0694                     0.0811                     0.1047

HMC and HMCDA samplers also give similarly bad results bud SMC actually works. Did I make something stupid here or is this a bug? A similar model in PyMC3 works with the NUTS sampler:

a1 = [0.05, 0.03, 0.02, 0.03, 0.04, 0.09, 0.13, 0.16]

with pm.Model() as model3:
    mu_a = pm.Normal('mu_a', 0.0, 0.5)
    aim = pm.Normal('aim', mu=mu_a, sigma=0.05, observed=a1)
    b = pm.Deterministic('b', -np.exp(-2*mu_a))

    tr = pm.sample(2000, tune=700)
    pm.traceplot(tr)
    display(pm.summary(tr))

  | mean | sd | hdi_3% | hdi_97% | mcse_mean | mcse_sd | ess_mean | ess_sd | ess_bulk | ess_tail | r_hat
-- | -- | -- | -- | -- | -- | -- | -- | -- | -- | -- | --
0.069 | 0.017 | 0.037 | 0.101 | 0.000 | 0.0 | 1860.0 | 1860.0 | 1851.0 | 2777.0 | 1.0
-0.872 | 0.030 | -0.923 | -0.811 | 0.001 | 0.0 | 1860.0 | 1860.0 | 1851.0 | 2777.0 | 1.0

[devmotion]

Isn't the logpdf implementation incorrect and should rather be something like

Distributions.logpdf(d::Determin, x) = x == d.T ? zero(float(x)) : oftype(float(x), -Inf)

? I haven't performed any debugging yet but I assume that the incorrect logpdf will lead HMC to return samples that do not adhere to the restriction that you want to impose by the Determin distribution. I'm wondering as well if maybe a DiscreteNonParametric distribution with support only at a single point would be more appropriate. Alternatively, you could try to use a Normal distribution with zero variance and check what happens instead of implementing your own distributions and bijectors.

Ultimatively, I don't think any of these ideas will yield a perfect solution since it seems you need rather a solution to TuringLang/DynamicPPL.jl#94 (the general problem in your example is a duplicate of this issue IMO). Letting HMC (or any other sampler) deal with b doesn't seem to be correct or helpful in general, you would rather want to just perform the computation AND track it.

So a maybe better solution for now would be to just perform these computations afterwards based on your samples of mu_a, and not during inference.

[devmotion]

BTW regardless of your problem with Determin you can simplify (and probably speed up) your model by implementing it as

@model function model1(a)
    ma ~ Normal(0.0, 0.5)
	a ~ filldist(Normal(ma, 0.05), length(a)) # alternatively: a .~ Normal(ma, 0.05)
    b ~ Determin(-exp(-2*ma))
    return ma, b
end

T is never used and hence not needed, and filldist should be faster than the loop.

[mateuszbaran]

Thanks! I was just trying to replicate Deterministic from PyMC3. It's actually just a small part of my full model where T is actually used but I didn't know about filldist.

I tried your logpdf and it doesn't work either. For now I'll just keep using PyMC3 and I'll port my model when Turing is more mature.

[devmotion]

I quickly ran your model with

    b ~ Normal(-exp(-2*ma), 0)

instead of

    b ~ Determin(-exp(-2*ma))

and as I assumed basically all samples of the HMC sampler are rejected (which makes inference impossible) due to the fact that the Hamiltonian dynamics don't respect the deterministic constraint. HMC doesn't seem to be suitable for sampling b (which ideally you wouldn't do at all...) - PG seems to work fine.

So I guess for now performing these computations on the resulting chain instead of inside the model might be the easiest option, if one wants to use Turing + HMC.

[torfjelde]

And to make what @devmotion suggested easier, you can potentially make use of the generated_quantities method from here: https://github.com/cambridge-mlg/Covid19/blob/4d967d58eeada3f3399116fc815cb17fd4163c87/src/utils.jl#L79.

With that, your code would end up something like:

@model model1(a, ::Type{T} = Float64) where {T} = begin
    ma ~ Normal(0.0, 0.5)
    for i in 1:length(a)
        a[i] ~ Normal(ma, 0.05)
    end

    b = -exp(-2*ma)
    return (b = b, ) # or just `b` if you don't want a named tuple
end

chain = sample(model, ...)
results = generated_quantities(model, chain) # <= results in an array of named tuples

You can also use the arrtup2tuparr method in the file I linked above to convert results into a slightly more convenient format.

[torfjelde]

For anyone checking out this issue now, the generated_quantities method is dependent on TuringLang/DynamicPPL.jl#147 to be merged to work with the most recent versions of Turing.

[torfjelde]

Update:

function generated_quantities(m::Turing.Model, c::MCMCChains.Chains)
    # if `c` is multiple chains we pool them into a single chain
    chain = length(chains(c)) == 1 ? c : MCMCChains.pool_chain(c)

    varinfo = Turing.DynamicPPL.VarInfo(m)

    return map(1:length(chain)) do i
        Turing.DynamicPPL.setval!(varinfo, chain, i, 1)
        m(varinfo)
    end
end

should now work nicely for all cases in the most recent release 👍

[yebai]

@torfjelde Maybe add this to the Turing docs?

[torfjelde]

I'd prefer to just add it to one of the packages, but it's a bit unclear which package it should go in. Once we figure that out, we could add it to the docs 👍

[trappmartin]

Hi, I think it would be great if we could support this functionality in Turing. I suppose Turing.jl is the best place for it. Moreover, I would suggest adding functionality to obtain the log prob values of each datum at each iteration after inference. I am not aware that we have this atm. This is useful for model selection, see: https://github.com/TuringLang/MCMCChains.jl#deviance-information-criterion-dic, and would be a great feature especially if it can easily be used for other models too.

[devmotion]

Maybe I'm missing something here, but doesn't generated_quantities above provide this functionality? So you are wondering where to put it? IMO it could also be put into DynamicPPL (which does not depend on MCMCChains though) or MCMCChains (which does not depend on DynamicPPL though) since it is not Turing.jl-specific - it doesn't depend on any objects or functionality in Turing, but only on stuff that belongs to MCMCChains and DynamicPPL. Maybe it would be possible to define it for AbstractChains and put it into DynamicPPL (DynamicPPL already depends on AbstractMCMC)? It seems the main thing would be to replace pool_chain with something more generic.

[trappmartin]

Sounds good to me too if that is possible.
Regarding log prob, I’m aware that this should be possible I was merely asking to add an additional interface function for users that does that.