missing values OK for simulation but not inference

[baggepinnen]

I am just starting out with some experiments using Soss, and came up with the following simple model of a review process with articles and reviewers. The model is happy about simulating forward when provided data with missing values, but the call to dynamicHMC fails with

ERROR: MethodError: no method matching logpdf(::Normal{ForwardDiff.Dual{ForwardDiff.Tag{LogDensityProblems.var"#34#35"{LogDensityProblems.Transform ...

nr = 5 # Number of reviewers
na = 5 # Number of articles
reviewer_bias = rand(Normal(0,1), nr)
article_score = rand(Normal(0,2), na)
R = clamp.([rand(Normal(r+a, 0.1)) for r in reviewer_bias, a in article_score], -5, 5)
Rmask = rand(Bool, size(R))
R = replace(Rmask, 0=>missing) .* R

m = @model begin
    rα ~ Normal(0,1) |> iid(nr)
    aα ~ Normal(0,2) |> iid(na)
    rσ ~ Gamma(0.2)
    R ~ For(1:na) do j
            For(1:nr) do i
                Normal(rα[i] + aα[j], rσ)
            end
    end
end;


Rc = [R[:,i] for i in 1:na] # Took me a while to realize I could not pass a matrix, I had to pass the columns. eachcol is not working btw.

s = rand(m(R=Rc))
s.R # no error but strange result

@time post = dynamicHMC(m(), (R=Rc,), 10);

[cscherrer]

Hi Fredrik,

Thanks for the issue. We can add missing value support to Soss (because it dispatches on data type), but right now it's mostly... missing ;) I would probably write this model using (i,j,R) triples rather than a dense matrix. You lose locality, but that makes it easy to separate observed vs unobserved data.

BTW you can build the matrix version of the generative model you describe like this:

m = @model begin
    rα ~ Normal(0,1) |> iid(nr)
    aα ~ Normal(0,2) |> iid(na)
    rσ ~ Gamma(0.2)
    R ~ For(nr, na) do i,j
        Normal(rα[i] + aα[j], rσ)
    end
end;

[baggepinnen]

Thanks for the input :) I'll try to rewrite the model to the triplets in the meantime, I'm not interested in inferring the missing R anyways.

With the matrix version of the model, if I do the following (without any missing)

s = rand(m(R=R))
s.R

shouldn't s.R == R? I get a completely different R back

julia> norm(s.R - R)
17.58012425152081

Full example:

nr = 5 # Number of reviewers
na = 5 # Number of articles
reviewer_bias = rand(Normal(0,1), nr)
article_score = rand(Normal(0,2), na)
R = clamp.([rand(Normal(r+a, 0.1)) for r in reviewer_bias, a in article_score], -5, 5)
# Rmask = rand(Bool, size(R))
# R = replace(Rmask, 0=>missing) .* R
m = @model begin
    reviewer_bias ~ Normal(0, 1) |> iid(nr)
    reviewer_gain ~ Normal(1, 0.1) |> iid(nr)
    true_article_score ~ Normal(0,2) |> iid(na)
    rσ ~ Gamma(0.2)
    R ~ For(nr, na) do i,j
        Normal(reviewer_bias[i] + true_article_score[j] + reviewer_gain[i]*true_article_score[j], rσ)
    end
end;
s = rand(m(R=R))
norm(s.R - R)

[cscherrer]

Good question. Currently (should we change this?) m(R=R) returns

julia> m(R=R)
Joint Distribution
    Bound arguments: [R]
    Variables: [rσ, rα, aα, R]

@model begin
        rσ ~ Gamma(0.2)
        rα ~ Normal(0, 1) |> iid(nr)
        aα ~ Normal(0, 2) |> iid(na)
        R ~ For(nr, na) do i, j
                Normal(rα[i] + aα[j], rσ)
            end
    end

So you can think of it as R being assigned and then overwritten.

That syntax is intended for assigning values to arguments (Soss models are function-like), but I can see how it could be confusing.

It might also help to see what code it's generating, which you can do like this:

julia> Soss.sourceRand()(m)
quote
    rσ = rand(Gamma(0.2))
    rα = rand(iid(nr, Normal(0, 1)))
    aα = rand(iid(na, Normal(0, 2)))
    R = rand(For(((i, j)->begin
                        #= REPL[2]:6 =#
                        Normal(rα[i] + aα[j], rσ)
                    end), (nr, na)))
    (rα = rα, aα = aα, rσ = rσ, R = R)
end

BTW for the sparse representation of your model I think you'll need to take something like an array of pairs of an argument. Let me know if you get stuck working through it, I'd be happy to have a closer look

[baggepinnen]

One of the reasons I expected R to be fixed to what I sent in was the example in the readme where that appears to be the case for X:

julia> X = randn(6,2)
6×2 Array{Float64,2}:
  1.19156    0.100793  
 -2.51973   -0.00197414
  2.07481    1.00879   
 -0.97325    0.844223  
 -0.101607   1.15807   
 -1.54251   -0.475159  
julia> truth = rand(m(X=X));

julia> pairs(truth)
pairs(::NamedTuple) with 3 entries:
  :X => [1.19156 0.100793; -2.51973 -0.00197414; … ; -0.101607 1.15807; -1.5425…

I am actually not sure what m(R=R) does if it does not give the conditional distribution given R, it seems passing R like that has no effect in this case? Would the result be different if R was used on the RHS in an expression before it was assigned?

On another note related to usability. I was thinking about the use of missing. In turing missing indicates that the value is not there and it will be inferred. Using indices and values separately I could ignore the missing values completely, which in my application was okay since I did not care about them. What about if I indicated them with a singleton ignore::Ignore which would tell Soss that those can be completely left out of the inference procedure? I.e:

struct Ignored end
const ignored = Ignored()

R = ... # Array with missing values
replace!(R, missing => ignored)
# use R in model

and that would be rewritten automatically to using indices and values. The benefit would be that there would be fewer variables to sample that if the missing values would be inferred and would keep the array on the original structure as seen by the user.

[cscherrer]

Hi @baggepinnen,

I was just going back through issues and saw this one, sorry to leave you hanging.

I am actually not sure what m(R=R) does if it does not give the conditional distribution given R, it seems passing R like that has no effect in this case?

Currently when you pass in a named tuple argument, Soss takes what it needs (only values listed as arguments in the model) and ignores the rest. We could consider changing this, the big concerns are

• What would give the most expected and widely useful behavior
• What would lead to fewer tricky bugs

For example, we could change things so that

m(;kwargs...) = Do(m; kwargs...)

or

m(;kwargs...) = predict(m; kwargs...)

Or we could throw an error, and insist that the arguments match. I can imagine any of these being surprising to some, so I'd love to talk through the possibilities and find the best option.

On the missing values thing, we can support this, just need to write the code. Soss generates code based on types, so it works out great that Missing is a different type.

OTOH, there are different kinds of missingness. Maybe this should be something like a model for the data, and another model for a mask. Then we could have a default mask for the simple MCAR case.

I haven't thought much about this, but I think there are some interesting things we could do.

[cscherrer]

I've started MaskArrays.jl to help with this: https://github.com/cscherrer/MaskArrays.jl