search

TuringLang / Turing.jl

Bayesian inference with probabilistic programming.
https://turinglang.org
MIT License
2.05k stars 219 forks source link

Bug in Gibbs sampler #99

Closed yebai closed 7 years ago

yebai commented 7 years ago 

using Turing, Distributions
using PyPlot, PyCall

K = 10
N = 51

T = [[0.272545164 6.5637631e-37 7.93970182e-6 8.61092244e-5 0.165656847 0.287081146 0.000218659173 4.06806091e-13 0.274399041 5.09466694e-6]       
 [0.0835696383 7.05396766e-6 0.206100213 1.52600546e-11 1.53908757e-8 6.5269808e-11 0.673013596 0.0193819808 0.0179275025 3.2272789e-18]       
 [0.886728867 8.63212428e-7 3.94046562e-11 0.110467781 1.11735715e-6 2.12946427e-18 5.16880247e-8 0.000236340361 0.00111600188 0.00144897701]  
 [8.16557898e-11 0.0168992485 6.27227283e-14 0.010767298 0.426756301 1.78050405e-9 6.41316248e-5 0.321159984 0.0369047216 0.187448314]       
 [0.130652141 4.70465408e-6 0.000473601393 0.0378509164 9.06618543e-10 0.000778622816 0.00029383779 0.829914494 8.69688804e-26 3.16817317e-5]  
 [0.223054659 0.288152163 7.35806925e-19 0.0185562602 1.73073908e-8 0.400936069 1.17437994e-12 0.0443974641 0.0249033671 7.2202285e-18]        
 [4.78064507e-16 0.320444079 0.00385904296 1.26156421e-9 0.00688364264 0.00447186979 0.156660567 0.169796226 0.333780163 0.00410440864]        
 [4.84470444e-13 2.5025163e-14 2.78748146e-7 0.00245132866 3.03033036e-12 0.00284425237 8.49830551e-7 6.60111797e-8 0.994702713 5.11768273e-7] 
 [0.16839645 1.80280379e-7 5.68958062e-11 0.134838199 0.00020810431 0.0861188042 0.0517409105 0.361825373 3.31239961e-11 0.196871978]          
 [0.0474764005 1.16126593e-6 5.96036112e-8 0.00128470373 1.30134792e-6 0.0374283978 0.310068428 2.27075277e-19 0.00647484474 0.597264703]] 

obs = [ 0.0,   7.72711051,   2.76189162,   8.8216901 ,
        10.80174329,   8.87655587,   0.47685358,   9.51892527,
         7.82538035,   5.52629325,  10.75167786,   5.94925434,
        -0.96912603,   1.65160838,   1.65005965,  -0.99642713,
         7.37803004,   5.40821392,   9.44046498,   8.51761132,
         9.76981763,   5.980154  ,   9.19558142,   5.33965621,
         6.2388448 ,   2.77755879,   6.67731151,   8.52411613,
        11.31057577,   8.11554144,   6.64705471,   8.02025435,
         9.84003587,   3.03943679,  -2.93966727,   2.04372567,
        -0.93734763,   3.66943525,   6.12876571,  -2.07758649,
         1.10420963,  -0.23197037,   3.64908206,  14.14671815,
         6.96651114,   7.28554932,   9.06049355,   6.54246834,
        11.22672275,   7.41962631,   8.45635411 ]

#means   = zeros(Float64,K)
initial = ones(Float64, K)/K

@model big_hmm(y) = begin
    states = tzeros(Int,N)
    means  = tzeros(Float64, K)

    for i = 1:K
        means[i] ~ Normal(0, 10)
    end
    states[1] ~ Categorical(initial)
    for i = 2:N
        states[i] ~ Categorical(vec(T[states[i-1],:]))
        y[i] ~ Normal(means[states[i]], 4)
    end
    res = means[states[:]]
    return(res, means)
end

g = Gibbs(100, HMC(2, 0.2, 5, :means), PG(10,2, :states))
@time chain = @sample(big_hmm(obs), g);

s = MambaChains(chain)

Running this example produces the following error:

KeyError: key :states not found
 in getindex at dict.jl:688 [inlined]
 in getindex(::Turing.Sample, ::Symbol) at io.jl:36
 in update(::Turing.VarInfo, ::Turing.Sample, ::Set{Symbol}) at gibbs_helper.jl:5
 in run(::Function, ::Dict{Any,Any}, ::Turing.GibbsSampler{Turing.Gibbs}) at gibbs.jl:67
 in sample(::Function, ::Dict{Any,Any}, ::Turing.Gibbs) at gibbs.jl:81
 in macro expansion; at util.jl:184 [inlined]
 in anonymous at <missing>:?
 in include_string(::String, ::String) at loading.jl:441
 in include_string(::String, ::String, ::Int64) at eval.jl:28
 in include_string(::Module, ::String, ::String, ::Int64, ::Vararg{Int64,N}) at eval.jl:32
 in (::Atom.##53#56{String,Int64,String})() at eval.jl:40
 in withpath(::Atom.##53#56{String,Int64,String}, ::String) at utils.jl:30
 in withpath(::Function, ::String) at eval.jl:46
 in macro expansion at eval.jl:57 [inlined]
 in (::Atom.##52#55{Dict{String,Any}})() at task.jl:60
yebai commented 7 years ago

UPDATE: after modifying the model to return states, I get the following error:

MethodError: Cannot `convert` an object of type Array{Int64,1} to an object of type ForwardDiff.Dual{N,T<:Real}
This may have arisen from a call to the constructor ForwardDiff.Dual{N,T<:Real}(...),
since type constructors fall back to convert methods.
 in copy!(::Base.LinearFast, ::Array{ForwardDiff.Dual,1}, ::Base.LinearFast, ::Array{Array{Int64,1},1}) at abstractarray.jl:559
 in vectorize(::Distributions.Categorical{Float64}, ::Array{Int64,1}) at hmc_helper.jl:28
 in update(::Turing.VarInfo, ::Turing.Sample, ::Set{Symbol}) at gibbs_helper.jl:7
 in run(::Function, ::Dict{Any,Any}, ::Turing.GibbsSampler{Turing.Gibbs}) at gibbs.jl:67
 in sample(::Function, ::Dict{Any,Any}, ::Turing.Gibbs) at gibbs.jl:81
 in macro expansion; at util.jl:184 [inlined]
 in anonymous at <missing>:?
 in include_string(::String, ::String) at loading.jl:441
 in include_string(::String, ::String, ::Int64) at eval.jl:28
 in include_string(::Module, ::String, ::String, ::Int64, ::Vararg{Int64,N}) at eval.jl:32
 in (::Atom.##53#56{String,Int64,String})() at eval.jl:40
 in withpath(::Atom.##53#56{String,Int64,String}, ::String) at utils.jl:30
 in withpath(::Function, ::String) at eval.jl:46
 in macro expansion at eval.jl:57 [inlined]
 in (::Atom.##52#55{Dict{String,Any}})() at task.jl:60

This error also occur in another mixture model example. It seems that the Gibbs sampler try to differentiate through a discrete variable.

xukai92 commented 7 years ago

Yeah. I saw the same error as well. I'm working on fix this.

102

yebai commented 7 years ago

Fixed by #103

yebai commented 7 years ago

After updating the model with a prior over transitioning matrix, I run into another error

using Turing, Distributions
using ForwardDiff: Dual
using PyPlot, PyCall

K = 10
N = 51

obs = [ 0.0,   7.72711051,   2.76189162,   8.8216901 ,
        10.80174329,   8.87655587,   0.47685358,   9.51892527,
         7.82538035,   5.52629325,  10.75167786,   5.94925434,
        -0.96912603,   1.65160838,   1.65005965,  -0.99642713,
         7.37803004,   5.40821392,   9.44046498,   8.51761132,
         9.76981763,   5.980154  ,   9.19558142,   5.33965621,
         6.2388448 ,   2.77755879,   6.67731151,   8.52411613,
        11.31057577,   8.11554144,   6.64705471,   8.02025435,
         9.84003587,   3.03943679,  -2.93966727,   2.04372567,
        -0.93734763,   3.66943525,   6.12876571,  -2.07758649,
         1.10420963,  -0.23197037,   3.64908206,  14.14671815,
         6.96651114,   7.28554932,   9.06049355,   6.54246834,
        11.22672275,   7.41962631,   8.45635411 ];

@model bayeshmm(y) = begin
    s = tzeros(Int64, N)
    m = tzeros(Dual, K)
    T = Array{Array}(K)
    for i = 1:K
        m[i] ~ Normal(0, 10)
        T[i] ~ Dirichlet(ones(K))
    end
    s[1] ~ Categorical(ones(Float64, K)/K)
    for i = 2:N
        s[i] ~ Categorical(vec(T[s[i-1]]))
        y[i] ~ Normal(m[s[i-1]], 4)
    end
    res = m[s[:]]
    return(res, s, m)
end

g = Gibbs(300, HMC(2, 0.2, 5, :m, :T), PG(100,2, :s))
c = @sample(bayeshmm(obs), g);

ERROR

MethodError: Cannot `convert` an object of type Array{ForwardDiff.Dual{0,Float64},1} to an object of type ForwardDiff.Dual{0,Float64}
This may have arisen from a call to the constructor ForwardDiff.Dual{0,Float64}(...),
since type constructors fall back to convert methods.

 in push!(::Array{ForwardDiff.Dual{0,Float64},1}, ::Array{ForwardDiff.Dual{0,Float64},1}) at ./array.jl:479
 in varInfo2samples(::Turing.VarInfo) at /Users/yebai/.julia/v0.5/Turing/src/samplers/support/gibbs_helper.jl:41
 in run(::Function, ::Dict{Any,Any}, ::Turing.GibbsSampler{Turing.Gibbs}) at /Users/yebai/.julia/v0.5/Turing/src/samplers/gibbs.jl:71
 in sample(::Function, ::Dict{Any,Any}, ::Turing.Gibbs) at /Users/yebai/.julia/v0.5/Turing/src/samplers/gibbs.jl:81
 in anonymous at ./<missing>:?