.~ seems to give incorrect answers

[xukai92]

using Turing

@model naive_bayes(image, label, D, N, C, ::Type{T}=Float64) where {T<:Real} = begin
    m = Matrix{T}(undef, D, C)
    for c = 1:C
        m[:,c] ~ MvNormal(fill(0, D), 10)
    end
    # m .~ Normal(0, 10) # this gives incorrect answer

    for n = 1:N
        image[:,n] ~ MvNormal(m[:,label[n]], 1)
    end
end

Checking answer via

using MLDatasets

image_raw = reshape(MNIST.traintensor(Float64), 28*28, :)
label = MNIST.trainlabels() .+ 1

# Pre-processing

using MultivariateStats

D_pca = 40
pca = fit(PCA, image_raw; maxoutdim=D_pca)

image = transform(pca, image_raw)

@info "Peformed PCA to reduce the dimension to $D_pca"

# Data function

get_data(n=100) = Dict(
    "C" => 10, 
    "D" => D_pca, 
    "N" => n, 
    "image" => image[:,1:n], 
    "label" => label[1:n]
)
data = get_data()
model = naive_bayes(data["image"], data["label"], data["D"], data["N"], data["C"])

alg = HMC(0.1, 4)
n_samples = 2_000

chain = sample(model, alg, n_samples)

m_data = Array(group(chain, :m))

m_bayes = mean(
    map(
        i -> reconstruct(pca, Matrix{Float64}(reshape(m_data[i,:], D_pca, 10))), 
        1_000:100:2_000
    )
)

using PyPlot

function make_imggrid(imgs, n_rows, n_cols; gap::Int=1)
    n = length(imgs)
    d_row, d_col = size(first(imgs))
    imggrid = 0.5 * ones(n_rows * (d_row + gap) + gap, n_cols * (d_col + gap) + gap)
    i = 1
    for row = 1:n_rows, col = 1:n_cols
        if i <= n
            i_row = (row - 1) * (d_row + gap) + 1
            i_col = (col - 1) * (d_col + gap) + 1
            imggrid[i_row+1:i_row+d_row,i_col+1:i_col+d_col] .= imgs[i]
        else
            break
        end
        i += 1
    end
    return imggrid
end
    
imgs = [reshape(m_bayes[:,i], 28, 28)' for i in 1:size(m_bayes, 2)]

plt.figure(figsize=(5, 2))
plt.imshow(make_imggrid(imgs, 2, 5), cmap="gray")

It should give something like below:

Edit (05/2021): Updated example

[torfjelde]

I think this is still an issue btw. Though we no longer give the wrong answer, we now just throw an error if the RHS of ~ is not a Vector (in your example it's a Matrix):

DynamicPPL.jl/src/compiler.jl

Lines 266 to 268 in aaab9bb

	top = [:($tmpright = $right),
	:($tmpright isa Union{$Distribution,AbstractVector{<:$Distribution}}
	|| throw(ArgumentError($DISTMSG)))]

Ran into this myself trying to help someone out on Discourse: https://discourse.julialang.org/t/bayesian-logistic-regression-with-turing-jl/60105/9

Is there a particular reason why we don't allow broadcasting over arbitrary arrays? Maybe @devmotion or @mohamed82008 remembers why we decided to enforce Vector only?

[mohamed82008]

I think at some point we could broadcast over arrays iirc

[mhauru]

Fixed by #804, which will be live in the next breaking release.