Networks of discrete random variables

[nathanielvirgo]

This is more of a question than an issue. (If there's a better place to ask it please let me know.)

I'm hoping to do inference on networks of (mostly) discrete random variables. ReactiveMP sounded like it might be well suited to my applications, because message passing is likely to perform a lot better than the MCMC algorithms that other PPLs tend to use.

However, I can't get some simple examples to work, and I'm wondering whether ReactiveMP can be used for this sort of thing or if I'm just barking up the wrong tree.

Here's the simplest thing I tried, modifying the coin flip example from the getting started page:

@model function coin_model()
    y = datavar(Float64, 1)
    θ ~ Bernoulli(0.5)
    y[1] ~ Bernoulli(θ)
    return y, θ
end

result = inference(
    model = Model(coin_model),
    data  = (y = Vector([1.0]), )
)

Here the 'parameter' θ is drawn from a Bernoulli distribution and y is just a copy of θ. We condition on y=1, so we hope to find that the posterior for θ is also concentrated on 1. However, it gives this error message instead (plus a stack trace):

No analytical rule available to compute a product of distributions Bernoulli{Float64}(p=0.5) and Beta{Float64}(α=2.0, β=1.0).

I understand this error, though it seems that calculating the product of those distributions analytically should be straightforward (you just evaluate the beta distribution at each point in the support of the Bernoulli one). Does this mean this kind of thing isn't supported?

I also tried some variations on this slightly less silly model, where y is a noisy copy of θ instead of an exact copy:

@model function coin_model()
    y = datavar(Float64, 1)
    θ ~ Bernoulli(0.5)
    p ~ 0.2 + (θ*0.6)    
    y[1] ~ Bernoulli(p)
    return y, θ
end

result = inference(
    model = Model(coin_model),
    data  = (y = Vector([1.0]), )
)

but I couldn't get any version of this to work either. (This version doesn't seem to like the +.)

In short my question is whether I can use ReactiveMP to build up networks of discrete random variables and perform inference on them, or is this not supported?

[albertpod]

Hi @nathanielvirgo!

Thanks for trying out ReactiveMP.
In short, the answer to your question is no.

At the moment ReactiveMP works mainly with conjugate models.
In the first example, when you happen to multiply Beta and Bernoulli distributions the resulting marginal can't be expressed in a form of known analytical distribution.
To circumvent this, you can define your own multiplication function, i.e.:
function prod(::ProdAnalytical, left::Bernoulli, right::Beta)

As for the second model, there are two problems, i.e. ReactiveMP doesn't know how (a) to multiply Bernoulli distribution with a PointMass 0.6 and (b) to sum "scaled" Bernoulli distribution with another PointMass of 0.2.
In principle, you can specify your own rules that will tackle this problem as explained in Creating custom nodes and message computation rules

In the near future, ReactiveMP will support "black box" inference techniques where you won't have to specify custom rules and products.

In case of urgency, you can check out ForneyLab.jl in particular, the notebook on Variational Laplace and Sampling which provides approximate posteriors in a wider range of models (this features will be ported to ReactiveMP soon, so stay tuned)

I guess @bvdmitri and @bartvanerp will elaborate more on this issue/question.

[bvdmitri]

Hey @nathanielvirgo !

Thanks for your question. As Albert already pointed out, the speed and efficiency of message passing (and ReactiveMP.jl implementation) comes from the fact that we run inference in large parts of the graph analytically, hence requiring conjugate structure of the model. Regarding your first model, there is no analytical solution of product of Bernoulli and Beta distributions. What we mean here is that there is no known distribution that represent product of Bernoulli and Beta. I'm not sure your proposal (you just evaluate the beta distribution at each point in the support of the Bernoulli one) is correct because it does not represent any distribution and, as I understand, the result is not normalisable (support of Bernoulli is {0, 1} and if you evaluate pdf of Beta in these points it gives either 0.0 or Inf)

However, ReactiveMP.jl provides you a straightforward API to inject approximation methods and run inference on non-conjugate models. This is somewhat advanced usage of our API but it is definitely possible. We support the @constraints macro just for that purpose. You can find an extensive documentation here. I will give you an idea how to use for your first example (slightly modified):

@model function coin_model(n)
    y = datavar(Float64, n)
    θ ~ NormalMeanVariance(0.0, 10.0) # I changed prior here, also non conjugate
    for i in 1:n
        y[i] ~ Bernoulli(θ)
    end
    return y, θ
end

If you run this model without constraints, it will give the similar error:

No analytical rule available to compute a product of distributions NormalMeanVariance{Float64}(μ=0.0, v=10.0) and Beta{Float64}(α=1.6155273539839612, β=1.3844726460160388).

However, with @constraints macro, we can use, e.g. SampleList approximation (docs):

constraints = @constraints begin 
     q(θ) :: SampleList(1000, RightProposal())
end

n = 10

result = inference(
    model = Model(coin_model, n),
    constraints = constraints,
    returnvars = (θ = KeepLast(), ),
    data  = (y = rand(n), )
)

mean(result.posteriors[:θ]) # some value here, 0.450794357590217 in my case

If you really want to use Bernoulli distribution as a prior I would suggest you to look into custom functional form constraints. You may implement you own custom approximation method for Bernoulli * Beta product (which, e.g., returns either Bernoulli(1.0) or Bernoulli(0.0) if you want to run inference with discrete values). Here is the extensive documentation for that: https://biaslab.github.io/ReactiveMP.jl/stable/custom/custom-functional-form/#custom-functional-form .

As for the second model, I'm not sure what do you mean by scaling and shifting (0.2 + (θ*0.6)) variable that has Bernoulli distribution. That is something, that ReactiveMP.jl does not know how to do. There is nothing wrong, though, if you know what you are doing you can indeed extend the list of possible rules of + and * factor nodes, as Albert suggested.