SMC sampling with Julia 1.1 produces samples with 0 variance for discrete RVs!

[trappmartin]

I observed this behaviour a few times already. At first I thought I couldn't reproduce the issue but now I figured out the corner case and here is a MWE.

using Turing
using LinearAlgebra
using StatsFuns: logistic

x = vcat(randn(100,3), randn(100,3).+2)
y = vcat(zeros(100), ones(100))

@model lr_nuts(x, y, σ) = begin

    N,D = size(x)

    α ~ Normal(0, σ)
    β ~ MvNormal(zeros(D), ones(D)*σ)

    p = Vector{Float64}(undef, D)
    s = Vector{Bool}(undef, D)
    for d = 1:D
        p[d] ~ Beta(1/2, 1/2)
        s[d] ~ Bernoulli(p[d])
    end

    for n = 1:N
        y[n] ~ Bernoulli(logistic(dot(x[n,s], β[s]) + α))
    end
end

chn = sample(lr_nuts(x,y,sqrt(10)), SMC(1000));

This code produces a chain with variance in the discrete RV's s if I'm using Julia 1.0.4. If I'm using Julia 1.1 instead, the variance for all discrete RV's (s) is zero.

[yebai]

This is likely not a bug - it appears that SMC is suffering from serious degeneracy issues here: s[1] and s[2] almost always collapsed into one value, while s[3] occasionally produce different values. Using a slightly more robust PG sampler would produce the following results

julia> chn = sample(lr_nuts(x,y,sqrt(10)), PG(100,1000))
[PG] Sampling...100% Time: 0:07:40
[ Info: [PG] Finished with
[ Info:   Running time    = 460.4168884250001;
Object of type Chains, with data of type 1000×12×1 Array{Union{Missing, Float64},3}

Log evidence      = 1.347086266753746e-15
Iterations        = 1:1000
Thinning interval = 1
Chains            = 1
Samples per chain = 1000
internals         = lp
parameters        = p[2], s[3], elapsed, β[3], s[1], s[2], p[1], α, β[2], p[3], β[1]

2-element Array{ChainDataFrame,1}

Summary Statistics

│ Row │ parameters │ mean     │ std       │ naive_se   │ mcse       │ ess     │ r_hat   │
│     │ Symbol     │ Float64  │ Float64   │ Float64    │ Float64    │ Any     │ Any     │
├─────┼────────────┼──────────┼───────────┼────────────┼────────────┼─────────┼─────────┤
│ 1   │ elapsed    │ 0.460417 │ 0.0381922 │ 0.00120774 │ 0.00472707 │ 30.3797 │ 1.08612 │
│ 2   │ p[1]       │ 0.521195 │ 0.230015  │ 0.00727372 │ 0.0731973  │ 4.60085 │ 1.42517 │
│ 3   │ p[2]       │ 0.911269 │ 0.0607229 │ 0.00192023 │ 0.0144352  │ 13.275  │ 1.03824 │
│ 4   │ p[3]       │ 0.340721 │ 0.0684537 │ 0.0021647  │ 0.0170574  │ 12.0584 │ 1.21026 │
│ 5   │ s[1]       │ 0.944    │ 0.230037  │ 0.0072744  │ 0.0548979  │ 25.8544 │ 1.06013 │
│ 6   │ s[2]       │ 0.998    │ 0.044699  │ 0.00141351 │ 0.002      │ 7.13199 │ 1.001   │
│ 7   │ s[3]       │ 0.994    │ 0.0772656 │ 0.00244335 │ 0.006      │ 6.49518 │ 1.00505 │
│ 8   │ α          │ -5.45379 │ 1.30607   │ 0.0413017  │ 0.423328   │ 5.06959 │ 1.31852 │
│ 9   │ β[1]       │ 1.00529  │ 0.311868  │ 0.00986212 │ 0.0660955  │ 18.3526 │ 1.10466 │
│ 10  │ β[2]       │ 1.95927  │ 0.298234  │ 0.00943099 │ 0.0680354  │ 8.18564 │ 1.19784 │
│ 11  │ β[3]       │ 1.9802   │ 0.478192  │ 0.0151217  │ 0.114609   │ 5.27897 │ 1.23839 │

Or running SMC with a larger number of particles:

julia> chn = sample(lr_nuts(x,y,sqrt(10)), SMC(10000))
Object of type Chains, with data of type 10000×12×1 Array{Union{Missing, Float64},3}

Log evidence      = -41.36666799787091
Iterations        = 1:10000
Thinning interval = 1
Chains            = 1
Samples per chain = 10000
internals         = lp
parameters        = p[2], s[3], β[3], s[1], s[2], p[1], α, le, β[2], p[3], β[1]

2-element Array{ChainDataFrame,1}

Summary Statistics

│ Row │ parameters │ mean     │ std         │ naive_se    │ mcse       │ ess     │ r_hat   │
│     │ Symbol     │ Float64  │ Float64     │ Float64     │ Float64    │ Any     │ Any     │
├─────┼────────────┼──────────┼─────────────┼─────────────┼────────────┼─────────┼─────────┤
│ 1   │ le         │ -41.3667 │ 7.10578e-15 │ 7.10578e-17 │ 0.0        │ 40.1606 │ 0.9999  │
│ 2   │ p[1]       │ 0.852187 │ 0.294901    │ 0.00294901  │ 0.0275693  │ 40.1606 │ 1.25    │
│ 3   │ p[2]       │ 0.710609 │ 0.115764    │ 0.00115764  │ 0.00894093 │ 45.1968 │ 1.07854 │
│ 4   │ p[3]       │ 0.728462 │ 0.243582    │ 0.00243582  │ 0.0216744  │ 40.1606 │ 1.15641 │
│ 5   │ s[1]       │ 0.8698   │ 0.33654     │ 0.0033654   │ 0.0313053  │ 40.5063 │ 1.16267 │
│ 6   │ s[2]       │ 0.9883   │ 0.107537    │ 0.00107537  │ 0.00588965 │ 351.587 │ 1.01181 │
│ 7   │ s[3]       │ 0.8169   │ 0.386768    │ 0.00386768  │ 0.0355487  │ 40.1606 │ 1.25597 │
│ 8   │ α          │ -7.21101 │ 0.949255    │ 0.00949255  │ 0.0865559  │ 40.1606 │ 1.22792 │
│ 9   │ β[1]       │ 2.04657  │ 1.61176     │ 0.0161176   │ 0.138484   │ 40.1606 │ 1.18215 │
│ 10  │ β[2]       │ 1.91447  │ 0.566133    │ 0.00566133  │ 0.0395129  │ 102.47  │ 1.03663 │
│ 11  │ β[3]       │ 1.20322  │ 1.84933     │ 0.0184933   │ 0.161411   │ 44.093  │ 1.08028 │

[trappmartin]

This is interesting! Thanks. I should have realised by myself.

However, it’s still strange that I got different behaviour depending on the Julia version. I’ll check if this was only a side effect.

[yebai]

Hi @trappmartin, is this issue resolved now?

[trappmartin]

I don't think so. But I give it another try.

[yebai]

Libtask is now tested on Julia 1.0, 1.1, 1.2 and nightly. So I don't think it's a task-copying related issue.

See https://travis-ci.org/TuringLang/Libtask.jl

[trappmartin]

Looks good, I'll still want to have a look at it before closing the issue.

[yebai]

I think this is resolved. Pls, reopen if there is a need.