Contigency factor node

[mateusjoffily]

Hi all,

Is there a Contingency factor node in RxInfer equivalent to the one in ForneyLab?

Thanks,
Mateus

[wouterwln]

Hi Mateus, there is not exactly the Contingency factor node. However, there is the general DiscreteTransition node, which can operate like a Contingency distribution. The DiscreteTransition models p(y | X) when y is categorical, and X is a collection of categorical random variables. We did this because the interface of always having normalized distributions coming out of the node is more convenient for RxInfer (if you are looking for "given some of my categoricals, what is the distribution of this other categorical, the DiscreteTransition node is your friend). What is the exact setting you are looking at? Maybe we can cast your Contingency problem to a DiscreteTransition problem. (It should be noted that DiscreteTransition is a renormalized version of a Contingency distribution, so you should be fine :)

[mateusjoffily]

Hi Wouter,

Thank you very much for your help. In my problem, y is binary and X is a collection of two binary variables. Please, see below the original ForneyLab code:

@rv x ~ Dirichlet(placeholder(:α, dims=(2,)))
y = Vector{Variable}(undef, 2)
for k = 1:2
@rv y[k] ~ Categorical(x)
end
Contingency(y[2], y[1], placeholder(:Y, dims=(2,2)))

Where can I find documentation on the DiscreteTransition node?

Thanks,
Mateus

[wouterwln]

I am not really sure what you want to achieve with the Contingency node here. From your explanation, it seems that y is actually a vector of 2 binary random variables. Do you want to learn the dependencies between the two different y's? I don't really see X back in your forneylab code, only x (which is Dirichlet and not a collection of two binary variables). The DiscreteTransition node is explained a bit in the HMM example and in the POMDP example. The exact documentation is not there yet, since the node is quite new and I haven't had time to document it.

What do you think of this model?

@model function transition_demo(y)
    A ~ DirichletCollection(ones(2,2))
    x ~ Dirichlet(ones(2))
    for i in 1:2
        y[i] ~ Categorical(x)
    end
    y[1] ~ DiscreteTransition(y[2], A)
end

transition_constraints = @constraints begin
       q(y,A) = q(y)q(A)
end

infer(model = transition_demo(), data = (y = [[0.0, 1.0], [1.0, 0.0]],), constraints = transition_constraints)

[mateusjoffily]

Sorry for not having been clearer, I should have renamed the variables in my code to match those in your first message. You're absolutely right about the description of my variables. If I had renamed them in my code, I would probably have realized that I misunderstood the use of the DiscreteTransition node.

In the ForneyLab code, I used the Contingency node to specify prior beliefs about the dependencies between the two binary variables, y[1] and y[2]. I know the parameters of the Contingency matrix, I don't need to learn them. I want to estimate the posterior beliefs about x given those prior beliefs about y (before observation). I would expect to get it using the command: result.posteriors[:x].alpha, after running infer().

I don't think this is what your proposed model does, but is there any way of implementing it in RxInfer?

It worked fine in ForneyLab. I write below the full code that I used in ForneyLab:

fg = FactorGraph()

@rv x ~ Dirichlet(placeholder(:α, dims=(2,)))
y = Vector{Variable}(undef, 2)
for k = 1:2
@rv y[k] ~ Categorical(x)
end
Contingency(y[2], y[1], placeholder(:Y, dims=(2,2)))

q = PosteriorFactorization(x, y, ids=[:Xm, :Ym])
algo = messagePassingAlgorithm(q, free_energy=true)
code = algorithmSourceCode(algo, free_energy=true)
eval(Meta.parse(code))

estimate conditional beliefs about x

data = Dict(:α => [1, 1],
:Y => [0.5 0; 0 0.5])
marginals = Dict{Symbol, ProbabilityDistribution}()
marginals[:x] = ProbabilityDistribution(Dirichlet, a=[1,1])
for i=1:n_its
stepYm!(data, marginals)
stepXm!(data, marginals)
end

update hyperparameters

α_p = marginals[:x].params[:a]

[wouterwln]

I kind of fail to see where the actual data you have about either x or y enters the model, but the following code puts a prior on y that models a joint dependency between the two (they are always equal):

@model function transition_demo(y_mat, α)
    x ~ Dirichlet(α)
    y[1] ~ Categorical(x)
    y[2] ~ Categorical(x)
    y[1] ~ DiscreteTransition(y[2], y_mat)
end

init = @initialization begin
    q(x) = Dirichlet([1.0, 1.0])
end

constraints = @constraints begin
    q(x, y) = q(x)q(y)
end

result = infer(model=transition_demo(α=[1, 1]), data=(y_mat=[1.0 0.0; 0.0 1.0],), iterations=10, initialization=init, constraints=constraints)

Adding an additional prior on y[2] allows you to fully model any contingency table:

@model function transition_demo(y_mat, α)
    x ~ Dirichlet(α)
    y[1] ~ Categorical(x)
    y[2] ~ Categorical(x)
    y[1] ~ DiscreteTransition(y[2], y_mat)
    y[2] ~ Categorical([0.8, 0.2])
end

init = @initialization begin
    q(x) = Dirichlet([1.0, 1.0])
end

constraints = @constraints begin
    q(x, y) = q(x)q(y)
end

result = infer(model=transition_demo(α=[1, 1]), data=(y_mat=[1.0 0.0; 0.0 1.0],), iterations=10, initialization=init, constraints=constraints)

Is this what you are looking for?

[mateusjoffily]

It worked perfectly! Thank you very much!

Regarding where the actual data enter my model, the data will be in y. After estimating the posterior beliefs about x conditional on the Contingency table, I replace the model with the one below:

@model function transition_demo_new(y, α)
x ~ Dirichlet(α)
y[1] ~ Categorical(x)
y[2] ~ Categorical(x)
end

result = infer(model=transition_demo_new(α=result.posteriors[:x][end].alpha), data=(y=[[0.0, 1.0], [1.0, 0.0]],))

where α is result.posteriors[:x][end].alpha) calculated with the previous model. In the new model I didn't impose the constraints. Perhaps there is a better way of achieving this?

Regarding your second model, I'm not sure I completely understood why an additional prior on y[2] would allow to fully model any contingency table. Could you give me any hints on this?

Thank you very much for your help :)

[wouterwln]

Hey Mateus, I see I missed your last reply, but I see that Thijs (author of forneylab) has modeled the Contingency node for you. I hope this is sufficient and I hope you will have fun using RxInfer!

[ThijsvdLaar]

@mateusjoffily I've ported the ForneyLab implementation of the Contingency node to RxInfer. In the example below the Contingency node directly models a joint distribution $p(y_1, y_2|B)$, normalized over all entries in $B$. In contrast, the DiscreteTransition node models a conditional distribution $p(y_1|y_2, B)$, normalized over the columns of $B$.

The node definition and rules are given by,

using RxInfer

using Base: BroadcastFunction
using ReactiveMP: normalize!
using LinearAlgebra: Diagonal, tr


@node Contingency Stochastic [ out1, out2, p ]

@average_energy Contingency (q_out1_out2::Contingency, q_p::Any) = begin
    return -tr(components(q_out1_out2)' * mean(BroadcastFunction(log), q_p))
end


#------------
# Sum-Product
#------------

@rule Contingency(:out1, Marginalisation) (m_out2::Categorical, m_p::PointMass) = begin
    p = mean(m_p) * probvec(m_out2)
    normalize!(p, 1)
    return Categorical(p)
end

@rule Contingency(:out2, Marginalisation) (m_out1::Categorical, m_p::PointMass) = begin
    p = mean(m_p)' * probvec(m_out1)
    normalize!(p, 1)
    return Categorical(p)
end


#------------
# Variational
#------------

@rule Contingency(:out1, Marginalisation) (m_out2::Categorical, q_p::Any) = begin
    p = clamp.(exp.(mean(BroadcastFunction(log), q_p) * probvec(m_out2)), tiny, Inf)
    normalize!(p, 1)
    return Categorical(p)
end

@rule Contingency(:out2, Marginalisation) (m_out1::Categorical, q_p::Any) = begin
    p = clamp.(exp.(mean(BroadcastFunction(log), q_p)' * probvec(m_out1)), tiny, Inf)
    normalize!(p, 1)
    return Categorical(p)
end


#----------
# Marginals
#----------

@marginalrule Contingency(:out1_out2) (m_out1::Categorical, m_out2::Categorical, q_p::Any) = begin
    D = map(e -> clamp(exp(e), tiny, huge), mean(BroadcastFunction(log), q_p))
    B = Diagonal(probvec(m_out1)) * D * Diagonal(probvec(m_out2))
    P = map!(Base.Fix2(/, sum(B)), B, B) # In-place normalization
    return Contingency(P, Val(false))
end

And a simple example with sumproduct message passing,

# An example for the contingency node with sum-product message passing.
# In this example, the contingency node models a joint distribution p(y_1, y_2 | B),
# while a discrete transition node would model the conditional p(y_1 | y_2, B).

@model function contingency_sp(B_hat)
    y[1] ~ Categorical([0.5, 0.5])
    y[2] ~ Categorical([0.5, 0.5])
    y[1] ~ Contingency(y[2], B_hat)
end

result = infer(model       = contingency_sp(), 
               data        = (B_hat=[0.4 0.1; 0.1 0.4],), # All entries in the matrix normalize
               free_energy = true)