Latent Variables

New to Julia
,
1

I am struggling with this error for modeling latent variables using Turing

MethodError: no method matching Normal(::Array{Float64,1}, ::Float64)

My data generation code is

using Distributions
using Statistics, LinearAlgebra
#using VegaDatasets
using Turing, TuringModels

function SEM_b(N_ID, NY)
    #loading matrix 2*6
    ld2 = [1 .71 .69 0 0 0;
    0 0 0 1 .35 1.5]

    cov2 = Matrix(1.0I,2,2)

    eta2 = MvNormal([0; 0],cov2)
    eta2b = rand(eta2, N_ID) #2*N matrix

    eps2 = MvNormal(fill(0, NY),diagm(ones(NY)))

    eps2b = rand(eps2, N_ID) #NY*N_ID matrix

    Y2 = transpose(eta2b)*ld2 + transpose(eps2b)

    transpose(eta2b), transpose(eps2b), Y2

end

and my Turing code is

@model MM_SEM(Y, N_eta) = begin

    ll,rr = size(Y)
    μ = Array{Float64,2}(undef,ll,rr)

    λ ~ filldist(Normal(0,5),4)
    η_σ ~ filldist(truncated(Normal(0,5),0,Inf),N_eta)

    σ ~ filldist(truncated(Normal(0,5),0,Inf),rr)

    η_ρ ~ LKJ(N_eta,1.)

    L = η_σ .* η_ρ .* η_σ'

    L = (L' + L)/2

    eta ~ filldist(MvNormal(fill(0,N_eta), L),ll)

    μ[:,1] = eta[1,:]
    μ[:,2] = λ[1]*eta[1,:]
    μ[:,3] = λ[2]*eta[1,:]

    μ[:,4] = eta[2,:]
    μ[:,5] = λ[3]*eta[2,:]
    μ[:,6] = λ[4]*eta[2,:]

    for j in 1:rr
        Y[:,j] ~ Normal(μ[:,j],σ[j])
    end
    #for i ∈ 1:ll

        #μ[i,1] = eta[1,i]
        #μ[i,2] = λ[1]*eta[1,i]
        #μ[i,3] = λ[2]*eta[1,i]

        #μ[i,4] = eta[2,i]
        #μ[i,5] = λ[3]*eta[2,i]
        #μ[i,6] = λ[4]*eta[2,i]

        #for j in 1:rr
        #    Y[i,j] ~ Normal(μ[i,j],σ[j])
        #end

    #end
end

and

chain = sample(MM_SEM(a[3],2), NUTS(0.95), 1000);
2

Here you calling Normal with an array as the first argument, and it seems there is no such method, only Normal with two scalars.

3

Thank you very much for your reply. Sorry, I am a bit confused at this moment. If I change it to

     for i ∈ 1:ll

        μ[i,1] = eta[1,i]
        μ[i,2] = λ[1]*eta[1,i]
        μ[i,3] = λ[2]*eta[1,i]

        μ[i,4] = eta[2,i]
        μ[i,5] = λ[3]*eta[2,i]
        μ[i,6] = λ[4]*eta[2,i]

        for j in 1:rr
            Y[i,j] ~ Normal(μ[i,j],σ[j])
        end

    end

I end up with this error also

TypeError: in typeassert, expected Float64, got a value of type ForwardDiff.Dual{Nothing,Float64,10}
in include_string at base\loading.jl:1088
in top-level scope at test.jl:78
in sample at Turing\RzDvB\src\inference\Inference.jl:154
in #sample#1 at Turing\RzDvB\src\inference\Inference.jl:154 
in sample at Turing\RzDvB\src\inference\Inference.jl:164 
in #sample#2 at Turing\RzDvB\src\inference\Inference.jl:164
in Sampler at Turing\RzDvB\src\inference\hmc.jl:376 
in DynamicPPL.Sampler at Turing\RzDvB\src\inference\hmc.jl:384
in Turing.Inference.HMCState at Turing\RzDvB\src\inference\hmc.jl:605
in #HMCState#58 at Turing\RzDvB\src\inference\hmc.jl:624
in find_good_stepsize at AdvancedHMC\P9wqk\src\trajectory.jl:778 
in #find_good_stepsize#13 at AdvancedHMC\P9wqk\src\trajectory.jl:778 
in find_good_stepsize##kw at AdvancedHMC\P9wqk\src\trajectory.jl:713 
in #find_good_stepsize#12 at AdvancedHMC\P9wqk\src\trajectory.jl:718
in phasepoint at AdvancedHMC\P9wqk\src\hamiltonian.jl:69 
in ∂H∂θ at AdvancedHMC\P9wqk\src\hamiltonian.jl:31 
in ∂logπ∂θ at Turing\RzDvB\src\inference\hmc.jl:474 
in gradient_logp at Turing\RzDvB\src\core\ad.jl:84 
in gradient_logp at Turing\RzDvB\src\core\ad.jl:84 
in gradient_logp at Turing\RzDvB\src\core\ad.jl:121
in gradient! at ForwardDiff\sdToQ\src\gradient.jl:33
in gradient! at ForwardDiff\sdToQ\src\gradient.jl:37 
in chunk_mode_gradient! at ForwardDiff\sdToQ\src\gradient.jl:140
in f at Turing\RzDvB\src\core\ad.jl:111 
in Model at DynamicPPL\uRBQJ\src\model.jl:84 
in Model at DynamicPPL\uRBQJ\src\model.jl:96 
in evaluate_threadsafe at DynamicPPL\uRBQJ\src\model.jl:135
in _evaluate at DynamicPPL\uRBQJ\src\model.jl:145 
in macro expansion at DynamicPPL\uRBQJ\src\model.jl 
in  at base\none
in #11 at test.jl:62 
in setindex! at base\array.jl:849
4

https://github.com/TuringLang/Turing.jl/issues/716

(Googling the error message helps many times)

5

Can you provide a MWE? I don’t know how to get a from your code.

Besides, it seems like you are adapting code from:
https://github.com/StatisticalRethinkingJulia/TuringModels.jl/blob/master/scripts/13/m13.1.jl

So do you try to write code like original form: wait .~ Normal.(mu , sigma) instead of an explicit looping?

6

Yes. I did adapt the part about using the LKJ correlation matrix and covariance matrix from the Julia version of Statistical Rethinking.
Sorry

a = SEM_b(300,6)

This is my first program in Julia, so I am not good at expressing my code in a condensed form.

Just to add, I am finally able to get it to work even though it is not efficient

@model function MM_SEM(Y, N_eta, ::Type{T} = Matrix{Float64}) where T

    ll,rr = size(Y)

    μ = T(undef, ll, rr)

    λ ~ filldist(Normal(0,5),4)
    η_σ ~ filldist(truncated(Normal(0,5),0,Inf),N_eta)

    σ ~ filldist(truncated(Normal(0,5),0,Inf),rr)

    η_ρ ~ LKJ(N_eta,1.)

    L = η_σ .* η_ρ .* η_σ'

    L = (L' + L)/2

    eta ~ filldist(MvNormal(fill(0,N_eta), L),ll)

    for i = 1:ll

        μ[i,1] = eta[1,i]
        μ[i,2] = λ[1]*eta[1,i]
        μ[i,3] = λ[2]*eta[1,i]

        μ[i,4] = eta[2,i]
        μ[i,5] = λ[3]*eta[2,i]
        μ[i,6] = λ[4]*eta[2,i]

        for j = 1:rr
            Y[i,j] ~ Normal(μ[i,j],σ[j])
        end

    end

    return λ, σ, η_σ, η_ρ
end

#need to check on the LKJ corr
chain = sample(MM_SEM(a[3],2), NUTS(0.95), 1000);