Need some help converting JAGS code into Turing Model

Specific Domains Probabilistic programming
1

I am trying to implement this JAGS model in Turing.

Screen Shot 2021-07-06 at 3.08.26 PM
Screen Shot 2021-07-06 at 3.08.26 PM1191×430 26.2 KB

JAGS code from:

Zhan, P., Jiao, H., Man, K., & Wang, L. (2019). Using JAGS for Bayesian cognitive diagnosis modeling: A tutorial. Journal of Educational and Behavioral Statistics , 44 (4), 473-503.

This is the code I tried to implement. Kind of stuck as unable to debug from the errors Turing is throwing (Newbie in Turing).

@model function DINA(X_Matrix, Q_Matrix)
    I = size(X_Matrix)[1]
    J,K = size(Q_Matrix)
    
    s ~ filldist(Beta(1,1),J)
    g ~ Truncated.(Beta(1,1), 0, 1 .- s)
    
    delta = fill(1,K)
    C = all_patterns([0,1],I)
    α = generate_α(C,K)
    η = compute_η(α, Q_Matrix)
    
    p = Matrix{Float64}(undef, I,J)
    for i in 1:I
        for j in 1:J
            p[i,j] = g[j] + (1 - s[j] - g[j]) * η[i,j]
        end
    end
    
    for i in 1:I
        for j in 1:J
            X_Matrix[i,j] ~ Bernouli(p[i,j])
        end
    end
end


function all_patterns(C,N)
    Matrix(transpose(hcat(collect.(reverse.(Iterators.product(fill(C,N)...))[:])...)))
end


function compute_η(A,Q)
    I = size(A)[1]
    J, K = size(Q)
    
    eta = typeof(A)(undef, I,J)
    for i in 1:I
        for j in 1:J 
            eta[i,j] = prod(A[i,:].^Q[j,:])
        end
    end
    return eta
end 

function generate_α(A, K)
    A[sample(begin:end, K), :]
end

I also tried the g parameter with simple Beta Distribution instead of truncated as part of debugging and it didn’t help either.

2 
using Turing
using LazyArrays
using StatsPlots
using StatisticalRethinking;
using DataFrames
using Optim
using StatsBase
using RCall

@rlibrary(CDM)
R"""
data(data.ecpe, package="CDM")
X <- data.ecpe$data[,-1]
Q <- data.ecpe$q.matrix
""";
@rget X;
@rget Q;

dina_model =DINA(Matrix(X), Matrix(Q))
sample(dina_model, MH(), MCMCThreads(),3000,3)

Error:

TaskFailedException

    nested task error: TaskFailedException
    Stacktrace:
     [1] wait
       @ .\task.jl:317 [inlined]
     [2] threading_run(func::Function)
       @ Base.Threads .\threadingconstructs.jl:34
     [3] macro expansion
       @ .\threadingconstructs.jl:93 [inlined]
     [4] macro expansion
       @ ~\.julia\packages\AbstractMCMC\oou1a\src\sample.jl:304 [inlined]
     [5] (::AbstractMCMC.var"#33#44"{Bool, Base.Iterators.Pairs{Symbol, UnionAll, Tuple{Symbol}, NamedTuple{(:chain_type,), Tuple{UnionAll}}}, Int64, Int64, Vector{Any}, Vector{UInt64}, Vector{DynamicPPL.Sampler{MH{(), NamedTuple{(), Tuple{}}}}}, Vector{DynamicPPL.Model{var"#23#24", (:X_Matrix, :Q_Matrix), (), (), Tuple{Matrix{Int64}, Matrix{Int64}}, Tuple{}}}, Vector{Random._GLOBAL_RNG}})()
       @ AbstractMCMC .\task.jl:406
    
        nested task error: MethodError: no method matching assume(::Random._GLOBAL_RNG, ::DynamicPPL.SampleFromPrior, ::Vector{Truncated{Beta{Float64}, Continuous, Float64}}, ::AbstractPPL.VarName{:g, Tuple{}}, ::DynamicPPL.ThreadSafeVarInfo{DynamicPPL.UntypedVarInfo{DynamicPPL.Metadata{Dict{AbstractPPL.VarName, Int64}, Vector{Distribution}, Vector{AbstractPPL.VarName}, Vector{Real}, Vector{Set{DynamicPPL.Selector}}}, Float64}, Vector{Base.RefValue{Float64}}})
        Closest candidates are:
          assume(::Any, ::Union{DynamicPPL.SampleFromPrior, DynamicPPL.SampleFromUniform}, ::Distribution, ::AbstractPPL.VarName, ::Any) at C:\Users\axs162731\.julia\packages\DynamicPPL\IKTvc\src\context_implementations.jl:119
          assume(::Any, ::DynamicPPL.Sampler{var"#s158"} where var"#s158"<:Union{PG, SMC}, ::Distribution, ::AbstractPPL.VarName, ::Any) at C:\Users\axs162731\.julia\packages\Turing\Gntg0\src\inference\AdvancedSMC.jl:314
          assume(::Any, ::DynamicPPL.Sampler{var"#s158"} where var"#s158"<:IS, ::Distribution, ::AbstractPPL.VarName, ::Any) at C:\Users\axs162731\.julia\packages\Turing\Gntg0\src\inference\is.jl:62
          ...
        Stacktrace:
          [1] _tilde(rng::Random._GLOBAL_RNG, sampler::DynamicPPL.SampleFromPrior, right::Vector{Truncated{Beta{Float64}, Continuous, Float64}}, vn::AbstractPPL.VarName{:g, Tuple{}}, vi::DynamicPPL.ThreadSafeVarInfo{DynamicPPL.UntypedVarInfo{DynamicPPL.Metadata{Dict{AbstractPPL.VarName, Int64}, Vector{Distribution}, Vector{AbstractPPL.VarName}, Vector{Real}, Vector{Set{DynamicPPL.Selector}}}, Float64}, Vector{Base.RefValue{Float64}}})
            @ DynamicPPL ~\.julia\packages\DynamicPPL\IKTvc\src\context_implementations.jl:59
          [2] tilde(rng::Random._GLOBAL_RNG, ctx::DynamicPPL.DefaultContext, sampler::DynamicPPL.SampleFromPrior, right::Vector{Truncated{Beta{Float64}, Continuous, Float64}}, vn::AbstractPPL.VarName{:g, Tuple{}}, #unused#::Tuple{}, vi::DynamicPPL.ThreadSafeVarInfo{DynamicPPL.UntypedVarInfo{DynamicPPL.Metadata{Dict{AbstractPPL.VarName, Int64}, Vector{Distribution}, Vector{AbstractPPL.VarName}, Vector{Real}, Vector{Set{DynamicPPL.Selector}}}, Float64}, Vector{Base.RefValue{Float64}}})
            @ DynamicPPL ~\.julia\packages\DynamicPPL\IKTvc\src\context_implementations.jl:23
          [3] tilde_assume(rng::Random._GLOBAL_RNG, ctx::DynamicPPL.DefaultContext, sampler::DynamicPPL.SampleFromPrior, right::Vector{Truncated{Beta{Float64}, Continuous, Float64}}, vn::AbstractPPL.VarName{:g, Tuple{}}, inds::Tuple{}, vi::DynamicPPL.ThreadSafeVarInfo{DynamicPPL.UntypedVarInfo{DynamicPPL.Metadata{Dict{AbstractPPL.VarName, Int64}, Vector{Distribution}, Vector{AbstractPPL.VarName}, Vector{Real}, Vector{Set{DynamicPPL.Selector}}}, Float64}, Vector{Base.RefValue{Float64}}})
            @ DynamicPPL ~\.julia\packages\DynamicPPL\IKTvc\src\context_implementations.jl:52
          [4] #23
            @ .\In[281]:6 [inlined]
          [5] (::var"#23#24")(__rng__::Random._GLOBAL_RNG, __model__::DynamicPPL.Model{var"#23#24", (:X_Matrix, :Q_Matrix), (), (), Tuple{Matrix{Int64}, Matrix{Int64}}, Tuple{}}, __varinfo__::DynamicPPL.ThreadSafeVarInfo{DynamicPPL.UntypedVarInfo{DynamicPPL.Metadata{Dict{AbstractPPL.VarName, Int64}, Vector{Distribution}, Vector{AbstractPPL.VarName}, Vector{Real}, Vector{Set{DynamicPPL.Selector}}}, Float64}, Vector{Base.RefValue{Float64}}}, __sampler__::DynamicPPL.SampleFromPrior, __context__::DynamicPPL.DefaultContext, X_Matrix::Matrix{Int64}, Q_Matrix::Matrix{Int64})
            @ Main .\none:0
          [6] macro expansion
            @ ~\.julia\packages\DynamicPPL\IKTvc\src\model.jl:0 [inlined]
          [7] _evaluate
            @ ~\.julia\packages\DynamicPPL\IKTvc\src\model.jl:154 [inlined]
          [8] evaluate_threadsafe(rng::Random._GLOBAL_RNG, model::DynamicPPL.Model{var"#23#24", (:X_Matrix, :Q_Matrix), (), (), Tuple{Matrix{Int64}, Matrix{Int64}}, Tuple{}}, varinfo::DynamicPPL.UntypedVarInfo{DynamicPPL.Metadata{Dict{AbstractPPL.VarName, Int64}, Vector{Distribution}, Vector{AbstractPPL.VarName}, Vector{Real}, Vector{Set{DynamicPPL.Selector}}}, Float64}, sampler::DynamicPPL.SampleFromPrior, context::DynamicPPL.DefaultContext)
            @ DynamicPPL ~\.julia\packages\DynamicPPL\IKTvc\src\model.jl:144
          [9] Model
            @ ~\.julia\packages\DynamicPPL\IKTvc\src\model.jl:94 [inlined]
         [10] VarInfo
            @ ~\.julia\packages\DynamicPPL\IKTvc\src\varinfo.jl:126 [inlined]
         [11] VarInfo
            @ ~\.julia\packages\DynamicPPL\IKTvc\src\varinfo.jl:125 [inlined]
         [12] step(rng::Random._GLOBAL_RNG, model::DynamicPPL.Model{var"#23#24", (:X_Matrix, :Q_Matrix), (), (), Tuple{Matrix{Int64}, Matrix{Int64}}, Tuple{}}, spl::DynamicPPL.Sampler{MH{(), NamedTuple{(), Tuple{}}}}; resume_from::Nothing, kwargs::Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
            @ DynamicPPL ~\.julia\packages\DynamicPPL\IKTvc\src\sampler.jl:73
         [13] step
            @ ~\.julia\packages\DynamicPPL\IKTvc\src\sampler.jl:66 [inlined]
         [14] macro expansion
            @ ~\.julia\packages\AbstractMCMC\oou1a\src\sample.jl:97 [inlined]
         [15] macro expansion
            @ ~\.julia\packages\AbstractMCMC\oou1a\src\logging.jl:15 [inlined]
         [16] mcmcsample(rng::Random._GLOBAL_RNG, model::DynamicPPL.Model{var"#23#24", (:X_Matrix, :Q_Matrix), (), (), Tuple{Matrix{Int64}, Matrix{Int64}}, Tuple{}}, sampler::DynamicPPL.Sampler{MH{(), NamedTuple{(), Tuple{}}}}, N::Int64; progress::Bool, progressname::String, callback::Nothing, discard_initial::Int64, thinning::Int64, chain_type::Type, kwargs::Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
            @ AbstractMCMC ~\.julia\packages\AbstractMCMC\oou1a\src\sample.jl:88
         [17] #sample#3
            @ ~\.julia\packages\Turing\Gntg0\src\inference\Inference.jl:155 [inlined]
         [18] macro expansion
            @ ~\.julia\packages\AbstractMCMC\oou1a\src\sample.jl:313 [inlined]
         [19] (::AbstractMCMC.var"#976#threadsfor_fun#45"{UnitRange{Int64}, Bool, Base.Iterators.Pairs{Symbol, UnionAll, Tuple{Symbol}, NamedTuple{(:chain_type,), Tuple{UnionAll}}}, Int64, Vector{Any}, Vector{UInt64}, Vector{DynamicPPL.Sampler{MH{(), NamedTuple{(), Tuple{}}}}}, Vector{DynamicPPL.Model{var"#23#24", (:X_Matrix, :Q_Matrix), (), (), Tuple{Matrix{Int64}, Matrix{Int64}}, Tuple{}}}, Vector{Random._GLOBAL_RNG}})(onethread::Bool)
            @ AbstractMCMC .\threadingconstructs.jl:81
         [20] (::AbstractMCMC.var"#976#threadsfor_fun#45"{UnitRange{Int64}, Bool, Base.Iterators.Pairs{Symbol, UnionAll, Tuple{Symbol}, NamedTuple{(:chain_type,), Tuple{UnionAll}}}, Int64, Vector{Any}, Vector{UInt64}, Vector{DynamicPPL.Sampler{MH{(), NamedTuple{(), Tuple{}}}}}, Vector{DynamicPPL.Model{var"#23#24", (:X_Matrix, :Q_Matrix), (), (), Tuple{Matrix{Int64}, Matrix{Int64}}, Tuple{}}}, Vector{Random._GLOBAL_RNG}})()
            @ AbstractMCMC .\threadingconstructs.jl:48

Stacktrace:
  [1] sync_end(c::Channel{Any})
    @ Base .\task.jl:364
  [2] macro expansion
    @ .\task.jl:383 [inlined]
  [3] macro expansion
    @ ~\.julia\packages\AbstractMCMC\oou1a\src\sample.jl:282 [inlined]
  [4] macro expansion
    @ ~\.julia\packages\ProgressLogging\6KXlp\src\ProgressLogging.jl:328 [inlined]
  [5] (::AbstractMCMC.var"#31#42"{Bool, String, Base.Iterators.Pairs{Symbol, UnionAll, Tuple{Symbol}, NamedTuple{(:chain_type,), Tuple{UnionAll}}}, Int64, Int64, Vector{Any}, Vector{UInt64}, Vector{DynamicPPL.Sampler{MH{(), NamedTuple{(), Tuple{}}}}}, Vector{DynamicPPL.Model{var"#23#24", (:X_Matrix, :Q_Matrix), (), (), Tuple{Matrix{Int64}, Matrix{Int64}}, Tuple{}}}, Vector{Random._GLOBAL_RNG}, UnitRange{Int64}})()
    @ AbstractMCMC ~\.julia\packages\AbstractMCMC\oou1a\src\logging.jl:11
  [6] with_logstate(f::Function, logstate::Any)
    @ Base.CoreLogging .\logging.jl:491
  [7] with_logger(f::Function, logger::LoggingExtras.TeeLogger{Tuple{LoggingExtras.EarlyFilteredLogger{ConsoleProgressMonitor.ProgressLogger, AbstractMCMC.var"#1#3"{Module}}, LoggingExtras.EarlyFilteredLogger{Base.CoreLogging.SimpleLogger, AbstractMCMC.var"#2#4"{Module}}}})
    @ Base.CoreLogging .\logging.jl:603
  [8] with_progresslogger(f::Function, _module::Module, logger::Base.CoreLogging.SimpleLogger)
    @ AbstractMCMC ~\.julia\packages\AbstractMCMC\oou1a\src\logging.jl:34
  [9] macro expansion
    @ ~\.julia\packages\AbstractMCMC\oou1a\src\logging.jl:10 [inlined]
 [10] mcmcsample(rng::Random._GLOBAL_RNG, model::DynamicPPL.Model{var"#23#24", (:X_Matrix, :Q_Matrix), (), (), Tuple{Matrix{Int64}, Matrix{Int64}}, Tuple{}}, sampler::DynamicPPL.Sampler{MH{(), NamedTuple{(), Tuple{}}}}, ::MCMCThreads, N::Int64, nchains::Int64; progress::Bool, progressname::String, kwargs::Base.Iterators.Pairs{Symbol, UnionAll, Tuple{Symbol}, NamedTuple{(:chain_type,), Tuple{UnionAll}}})
    @ AbstractMCMC ~\.julia\packages\AbstractMCMC\oou1a\src\sample.jl:276
 [11] sample(rng::Random._GLOBAL_RNG, model::DynamicPPL.Model{var"#23#24", (:X_Matrix, :Q_Matrix), (), (), Tuple{Matrix{Int64}, Matrix{Int64}}, Tuple{}}, sampler::DynamicPPL.Sampler{MH{(), NamedTuple{(), Tuple{}}}}, parallel::MCMCThreads, N::Int64, n_chains::Int64; chain_type::Type, progress::Bool, kwargs::Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ Turing.Inference ~\.julia\packages\Turing\Gntg0\src\inference\Inference.jl:216
 [12] sample
    @ ~\.julia\packages\Turing\Gntg0\src\inference\Inference.jl:216 [inlined]
 [13] #sample#6
    @ ~\.julia\packages\Turing\Gntg0\src\inference\Inference.jl:201 [inlined]
 [14] sample
    @ ~\.julia\packages\Turing\Gntg0\src\inference\Inference.jl:201 [inlined]
 [15] #sample#5
    @ ~\.julia\packages\Turing\Gntg0\src\inference\Inference.jl:188 [inlined]
 [16] sample(model::DynamicPPL.Model{var"#23#24", (:X_Matrix, :Q_Matrix), (), (), Tuple{Matrix{Int64}, Matrix{Int64}}, Tuple{}}, alg::MH{(), NamedTuple{(), Tuple{}}}, parallel::MCMCThreads, N::Int64, n_chains::Int64)
    @ Turing.Inference ~\.julia\packages\Turing\Gntg0\src\inference\Inference.jl:188
 [17] top-level scope
    @ In[284]:1
 [18] eval
    @ .\boot.jl:360 [inlined]
 [19] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)
    @ Base .\loading.jl:1094
3

Try changing

g ~ Truncated.(Beta(1,1), 0, 1 .- s)

to

g ~ array_dist(Truncated(.Beta(1,1), 0, 1 .- s))

And see where it gets you. It’s also a little easier to debug this if you don’t use multithreading at the start – i.e. use sample(dina_model, MH(), 3000) just for debugging purposes.

1 Like
4

Hi. There are some wrong points in your Turing model.

I have tried to revise them.

@model function DINA(X_Matrix, Q_Matrix)
    N, I = size(X_Matrix)
    I, K = size(Q_Matrix)
    
    s ~ filldist(Beta(1, 1), I)
    g ~ arraydist(Truncated.(Beta(1, 1), 0, 1 .- s))
    
    aap = all_patterns([0, 1], K)
    C = size(aap, 1)
    delta = fill(1, C) # For Dirichlet
    _π ~ Dirichlet(delta)
    c ~ filldist(Categorical(_π), N)
    η = [prod(aap[c[n], :] .^ Q_Matrix[i,:]) for n in 1:N, i in 1:I]

    p = Matrix{Float64}(undef, N, I)
    for n in 1:N
        for i in 1:I
            p[n, i] = g[i] + (1 - s[i] - g[i]) * η[n, i]
        end
    end
    
    for n in 1:N
        for i in 1:I
            X_Matrix[n, i] ~ Bernoulli(p[n, i])
        end
    end
end