array comprehension for random variables

[scheidan]

Defining an array of random variables (RV) with array comprehension leads to different error messages depending how it is used.

Model without array comprehension:

@model m1(y) = begin
    N = length(y)
    b = Array{Real}(undef, 4)
    for k in 1:4
        b ~ Normal(0, 0.2)
    end

    for i in 1:N
        y[1] ~ Normal(sum(b), 1.0)
    end
end

-> works

using array comprehension for unobserved RV:

@model m2(y) = begin
    N = length(y)

    b ~ [Normal(0, 0.2) for i in 1:4]

    for i in 1:N
        y[1] ~ Normal(sum(b), 1.0)
    end
end

-> ERROR: LoadError: UndefVarError: b not defined

using array comprehension for observed RV

@model m3(y) = begin
    N = length(y)
    b = Array{Real}(undef, 4)
    for k in 1:4
        b ~ Normal(0, 0.2)
    end

    y ~ [Normal(sum(b), 1.0) for i in 1:N]
end

-> ERROR: LoadError: AssertionError: Turing.observe only support vectorizing i.i.d distribution

tested on Turing v0.5.1+ #master

[trappmartin]

The following syntax for array comprehension for random variables is currently supported in Turing.

# Model specification.
@model m() = begin
    x = Vector{Float64}(undef, 10)
    x ~ [Normal(0, 1)]
    return x
end

# Create model function.
mf = m1()

# Draw sample from prior.
y = mf()

This syntax, however, seems inconvenient to me.

@xukai92, @yebai I suppose we have this syntax because x needs to be preallocated and the array comprehension was conveniently short. However, in favour or clarity, I would prefer to change the syntax to something more standard Julia, e.g.

# Model specification.
@model m() = begin
    N = 10
    x = Vector{Float64}(undef, N)
    x ~ [Normal(0, 1) for n in 1:N]
    return x
end

or

# Model specification.
@model m() = begin
    N = 10
    x = Vector{Float64}(undef, N)
    x[n] ~ [Normal(0, 1) for n in 1:N]
    return x
end

or is there a strong reason not to do that?

[mohamed82008]

Please use broadcasting now, .~, e.g:

@model m() = begin
    x = Vector{Float64}(undef, 10)
    x .~ Normal(0, 1)
    return x
end

Solved by #965.