I've tried to reproduce the model from a PYMC3 and Stan comparison. But it seems to run slowly and when I look at @code_warntype there are some things -- K and N I think -- which the compiler seemingly calls Any.

I've tried adding types -- though I can't add types to turing_model's arguments and things are complicated within turing_model because it's using autodiff variables and not the usuals. I put all the code into the function do_it to avoid globals, because they say that globals can slow things down. (It actually seems slower, though.)

Any suggestions as to what's causing the problem? The turing_model code is what's iterating, so that should make the most difference.

using Turing, StatsPlots, Random


function scale(w0::Float64, w1::Array{Float64,1})
    scale = √(w0^2 + sum(w1 .^ 2))
    return w0 / scale, w1 ./ scale

function do_it(iterations::Int64)::Chains
    K = 10                         # predictor dimension
    N = 1000                       # number of data samples

    X = rand(N, K)                 # predictors    (1000, 10)

    w1 = rand(K)                   # weights       (10,)
    w0 = -median(X * w1)           # 50% of elements for each class (number)
    w0, w1 = scale(w0, w1)         # unit length (euclidean)
    w_true = [w0, w1...]

    y = (w0 .+ (X * w1)) .> 0.0    # labels
    y = [Float64(x) for x in y]

    σ = 5.0
    σm = [x == y ? σ : 0.0 for x in 1:K, y in 1:K]

    @model turing_model(X, y, σ, σm) = begin

        w0_pred ~ Normal(0.0, σ)
        w1_pred ~ MvNormal(σm)

        p = sigmoid.(w0_pred .+ (X * w1_pred))

        @inbounds for n in 1:length(y)
            y[n] ~ Bernoulli(p[n])

    @time chain = sample(turing_model(X, y, σ, σm), NUTS(iterations, 200, 0.65));

    # ϵ = 0.5
    # τ = 10

    # @time chain = sample(turing_model(X, y, σ), HMC(iterations, ϵ, τ));

    return (w_true=w_true, chains=chain::Chains)

chain = do_it(1000)
  • This doesn't run as written. Can you add enough code to let the example be run? It makes it much easier to diagnose performance problems. Thanks – Oscar Smith May 30 at 3:40
  • Sorry, I'd forgotten the two supporting functions which are used in setting up the data. Should run now. – Wayne May 30 at 11:58
  • The type-instability does not appear about N or K. N and K is set to an Int64 and you never change its binding. Although you set y to a BitArray and then a Float64 array, Julia gets rid of the instability issue during compilation. The type-stability is due to the model. I am not sure how Turing.jl works but probably this type-instability is expected. I would try asking this question on Discourse and give a link to here. – hckr May 31 at 19:56

Your Answer


Browse other questions tagged or ask your own question.