ForneyLab.step! must be explicitly imported

[erlebach]

I am trying to duplicate results in Winn and Bishop's book "Model-based machine learning" and feel that ForneyLab can accomplish this task. Here is some code:

# Calculate the posterior of p_skills[i], i = {1,2}
g = FactorGraph()
T = [0.9 0.1; 0.2 0.8] |> transpose

nb_skills = 2
nb_questions = 3
p_skills = Vector{Variable}(undef, nb_skills) # one-hot coding
is_question_correct = Vector{Variable}(undef, nb_questions)

a = Vector{Int32}(undef, nb_skills)
b = Vector{Int32}(undef, nb_skills)
p4 = Vector{Variable}(undef, 2)

for i in 1:nb_skills
    @RV p_skills[i] ~ Beta(a[i], b[i])  # p_skills prior
end

# The first argument to Transition must be a vector
# Is my use of Transition allowed?
@RV is_question_correct[1] ~ Transition([p_skills[1], 1-p_skills[1]], T)
@RV is_question_correct[2] ~ Transition([p_skills[2], 1-p_skills[2]], T)
@RV p4[1] = p_skills[1] * p_skills[2]    # Deterministic function
@RV p4[2] = 1 - p4[1]

# p4 is a 2-vector
@RV is_question_correct[3] ~ Transition(p4, T)    # Requires that p4 have type Vector{Variable}

# Observed variables
placeholder(is_question_correct[1], :is_correct1, index=1, datatype=Int)
placeholder(is_question_correct[2], :is_correct2, index=2, datatype=Int)
placeholder(is_question_correct[3], :is_correct3, index=3, datatype=Int)

There are no errors so far. I am working in Jupiter-lab 3.3 and Julia 7.2 and the latest version of ForneyLab 0.11.4 on the Mac Air M1. Now I create some source code:

algo = messagePassingAlgorithm(p_skills[1])
source_code = algorithmSourceCode(algo)
source_expr = Meta.parse(source_code) # Parse and load the algorithm in scope
println(source_code) # Inspect the algorithm codej

I get the message:

KeyError: key nothing not found

Stacktrace:
  [1] getindex(h::Dict{Interface, Type}, key::Nothing)
    @ Base ./dict.jl:481
  [2] collectInboundTypes(entry::ScheduleEntry, #unused#::Type{SumProductRule{Addition}}, inferred_outbound_types::Dict{Interface, Type})
...

I looked everywhere for clues on how to solve this, unsuccessfully so far. I would appreciate any help on this. Thanks!

[albertpod]

Hi @erlebach!

Could you please point out which model from the book you are referring to? This will help to better understand the model.

[erlebach]

The model is shown in figure 2.14, page 74, in Chapter 2. I removed the 4th question.

[albertpod]

Hi @erlebach!

I see. ForneyLab.jl can run inference in this model, however, it won't rely on loopy-belief propagation. (FL doesn't support this type of inference algorithm as opposed to ReactiveMP.jl.

I think that for the AND operator you may want to use the Contingency node that implements cross-tabulation.

We will provide you with the solution by the end of the week. At the moment all developers of ForneyLab.jl are busy.

[erlebach]

Sounds great. In the meantime, I will check out the Contingency lab. Note that the network in my example has no loops. I guess another solution would be to reformulate the problem with continuous distributions? Many Complex systems are cyclic aren't they? I guess loopy propagation is hard to automate.

[albertpod]

I see. I will have a closer look. Reformulating the problem in terms of continuous distributions will ease the life of ForneyLab indeed.

As for real-life scenarios you are right: cycles are usually unavoidable, however, there are ways to circumvent this, e.g. variational message passing. As for loopy BP, ReactiveMP.jl has this algorithm, however, let's focus on the ForneyLab's solution at the moment.

We'll get back to you soon.

[erlebach]

Thanks.

[semihakbayrak]

Hi @erlebach!

For the time being, let me provide you with an implementation that executes approximate inference with importance sampling on the model Fig. 2.13. I guess this is the model you want to implement.

using ForneyLab

graph = FactorGraph()

f1(x) = x*0.9 + (1-x)*0.1
f2(x) = x*0.2 + (1-x)*0.8
f3(x,z) = 0.9999*x*z + 0.0001*(1-x*z) #soften the AND function a bit to avoid numerical problems
@RV csharp ~ Bernoulli(0.5)
@RV sql ~ Bernoulli(0.5)
@RV s1 ~ Nonlinear{Sampling}(csharp,g=f1)
@RV s2 ~ Nonlinear{Sampling}(sql,g=f2)
@RV s3 ~ Nonlinear{Sampling}(csharp,sql,g=f3)
@RV isCorrect1 ~ Bernoulli(s1)
@RV isCorrect2 ~ Bernoulli(s2)
@RV isCorrect3 ~ Bernoulli(s3)
placeholder(isCorrect1, :isCorrect1)
placeholder(isCorrect2, :isCorrect2)
placeholder(isCorrect3, :isCorrect3)

algo = messagePassingAlgorithm([csharp,sql])
source_code = algorithmSourceCode(algo)
eval(Meta.parse(source_code))
println(source_code)

data = Dict(:isCorrect1 => 1., :isCorrect2 => 0., :isCorrect3 => 0.) # example query
marginals = step!(data)

After executing the inference, you can check the mean values of the approximate posteriors with the followings

mean(marginals[:sql]) mean(marginals[:csharp]) You'll see that the former is approximately 0.26 while the letter one is approximately 0.66 with these configurations.

Please note that Bernoulli has one parameter in ForneyLab instead of two, as one parameter is sufficient to handle Binary cases. So please be careful with that also while using ReactiveMP package.

[erlebach]

That is so kind of you! I will study this carefully and learn from it!

I have simplified the code with the function I believe is the correct one: addNoise. Here is simplified code and a new error, which I do not understand.

using ForneyLab

graph = FactorGraph()

function addNoise(has_skill)
    if has_skill == 1.0
        res = sample(ProbabilityDistribution(Bernoulli, p=0.9))
        # f = x*0.9 + (1-x)*0.1
    else
        res = sample(ProbabilityDistribution(Bernoulli, p=0.2))
        # f = x*0.2 + (1-x)*0.8
    end
    return res
end

@RV csharp ~ Bernoulli(0.5)
@RV isCorrect1 ~ Nonlinear{Sampling}(csharp, g=addNoise, n_samples=50)
placeholder(isCorrect1, :isCorrect1)

algo = messagePassingAlgorithm([csharp,sql])
source_code = algorithmSourceCode(algo)
eval(Meta.parse(source_code))
println(source_code)

The error message:

MethodError: no method matching logPdf(::ProbabilityDistribution{Univariate, PointMass}, ::Float64)
Closest candidates are:
  logPdf(::ProbabilityDistribution{<:ForneyLab.VariateType, Function}, ::Any) at ~/src/2022/ForneyLab.jl/src/probability_distribution.jl:132
  logPdf(::ProbabilityDistribution{MatrixVariate, Dirichlet}, ::Any) at ~/src/2022/ForneyLab.jl/src/factor_nodes/dirichlet.jl:71
  logPdf(::ProbabilityDistribution{Multivariate, Dirichlet}, ::Any) at ~/src/2022/ForneyLab.jl/src/factor_nodes/dirichlet.jl:70
  ...

Stacktrace:
  [1] (::ForneyLab.var"#322#323"{typeof(addNoise), Message{PointMass, Univariate}})(z::Float64)
    @ ForneyLab ~/src/2022/ForneyLab.jl/src/engines/julia/update_rules/nonlinear_sampling.jl:47
  [2] logPdf(dist::ProbabilityDistribution{Univariate, Function}, x::Float64)
    @ ForneyLab ~/src/2022/ForneyLab.jl/src/probability_distribution.jl:132
  [3] _broadcast_getindex_evalf
    @ ./broadcast.jl:670 [inlined]
  [4] _broadcast_getindex
    @ ./broadcast.jl:643 [inlined]
  [5] getindex
    @ ./broadcast.jl:597 [inlined]
  [6] copy
    @ ./broadcast.jl:899 [inlined]
  [7] materialize(bc::Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{1}, Nothing, typeof(logPdf), Tuple{Vector{ProbabilityDistribution{Univariate, Function}}, Vector{Float64}}})
    @ Base.Broadcast ./broadcast.jl:860
  [8] sampleWeightsAndEntropy(x::ProbabilityDistribution{Univariate, Bernoulli}, y::ProbabilityDistribution{Univariate, Function})
    @ ForneyLab ~/src/2022/ForneyLab.jl/src/factor_nodes/sample_list.jl:207
  [9] prod!(x::ProbabilityDistribution{Univariate, Bernoulli}, y::ProbabilityDistribution{Univariate, Function}, z::ProbabilityDistribution{Univariate, SampleList}) (repeats 2 times)
    @ ForneyLab ~/src/2022/ForneyLab.jl/src/factor_nodes/sample_list.jl:232
 [10] *(x::ProbabilityDistribution{Univariate, Bernoulli}, y::ProbabilityDistribution{Univariate, Function})
    @ ForneyLab ~/src/2022/ForneyLab.jl/src/ForneyLab.jl:132
 [11] step!(data::Dict{Symbol, Float64}, marginals::Dict{Any, Any}, messages::Vector{Message})
    @ Main ./none:8
 [12] step! (repeats 2 times)
    @ ./none:5 [inlined]
 [13] top-level scope
    @ In[19]:2
 [14] eval
    @ ./boot.jl:373 [inlined]
 [15] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)
    @ Base ./loading.jl:1196

Debugging is quite difficult!

[semihakbayrak]

Hi @erlebach!

I have a better understanding of the model now. I'll put the code below with explanations and step by step we will build the model in Figure 2.5 by ForneyLab. Please see comments in the code and let me know if there is anything not clear.

using ForneyLab

graph = FactorGraph()

# Let us go step by step
# Page 53 Assumption 2 postulates that there is an equal chance that this person is knowledgable in csharp and sql
# We insert this prior knowledge through Bernoulli distributions with parameters 0.5 as they do in the book

@RV csharp ~ Bernoulli(0.5)
@RV sql ~ Bernoulli(0.5)

# So far, we created the Figure 2.2. Now, let us create the AddNoise function.
# If the person has that skill than %90 correct answer chance, if not %20 correct answer chance.
# Note that this case refers to setting the parameter of Bernoulli likelihood node to 0.9 and 0.2 respectively.
# Let's do that below
AddNoise(hasSkill) = if hasSkill == 1 return 0.9 else return 0.2 end

# Now we are ready to build Figure 2.3
@RV s1 ~ Nonlinear{Sampling}(csharp,g=AddNoise,n_samples=50)
@RV s2 ~ Nonlinear{Sampling}(sql,g=AddNoise,n_samples=50)

@RV isCorrect1 ~ Bernoulli(s1)
@RV isCorrect2 ~ Bernoulli(s2)

# We are done with Figure 2.3. Our model will do exactly the same thing they describe in the book. We just adjust the model
# specification in such a way that ForneyLab can execute inference.

# Let us implement Figure 2.4, as well. By starting implementing the AND function
ANDfunc(csharp,sql) = csharp*sql
# Now compose AND and AddNoise functions so as to implement Figure 2.5
Compositefunc(csharp,sql) = AddNoise(ANDfunc(csharp,sql))
# Let us implement a mediator variable that will define the parameter of Bernoulli likelihood
@RV s3 ~ Nonlinear{Sampling}(csharp,sql,g=Compositefunc,n_samples=50)
# Connecting this to a Bernoulli likelihood node completes the definition for question 3
@RV isCorrect3 ~ Bernoulli(s3)

# Lastly, let us inform ForneyLab that we will run queries by the following placeholders
placeholder(isCorrect1, :isCorrect1)
placeholder(isCorrect2, :isCorrect2)
placeholder(isCorrect3, :isCorrect3)

# Algorithm generation
algo = messagePassingAlgorithm([csharp,sql])
source_code = algorithmSourceCode(algo)
eval(Meta.parse(source_code))
println(source_code)

# Let us run an example query
data = Dict(:isCorrect1 => 1., :isCorrect2 => 0., :isCorrect3 => 1.) # example query
marginals = step!(data)

# Now we can fetch the weighted samples for csharp and sql variables as the following
marginals[:sql]
marginals[:csharp]

# It is possible to calculate statistics like mean
mean(marginals[:sql])
mean(marginals[:csharp])

So this completes the model specification. Again, this model specification I think same with what they describe in the book. The way they define likelihoods is through tabular probability table consisting of [0.9 0.1;0.2 0.8]. We achieve the same specification by directly putting Bernoulli likelihoods, parameters of which are defined by AddNoise function.

Note that you do not need to explicitly sample from distributions as you did in your model specification. Setting the parameters of Bernoulli is enough and ForneyLab will handle sampling for us.

[erlebach]

Thanks! It all makes sense. The question I have is why did my code fail? It looks identical to your version except for the definition of the AddNoise function. Is it that ForneyLab could not handle me AddNoise function? The error messages are quite cryptic as a result of meta-programming I guess. I will work backwards from your version. Thanks!

[semihakbayrak]

In my opinion ForneyLab can handle your AddNoise function, too. I see two main reasons why it didn't work for you. First of all you try to infer sql variable that you didn't define in the model: algo = messagePassingAlgorithm([csharp,sql]). Secondly, and more importantly, your likelihood node is directly a deterministic function (Nonlinear node). In ForneyLab, we do not allow this yet. Likelihood nodes need to be soft factors (probability distributions other than kronecker delta, e.g., GaussianMeanVariance, Poisson, Beta, Bernoulli, etc.). In my code, you'll see that I use Bernoulli distributions as likelihood nodes. That's why I define the model on the book a little bit differently in syntax but same in effect: to make ForneyLab compatible with the model specification.

[erlebach]

Hi @semihakbayrak !

I did not notice your reply until after I sent this message. My questions are partially answered, especially the most important one: likelihood functions must be soft functions. But that does not explain why your approach works, and obviously, it does. You applied your Nonlinear function to a mixture of p=0.9 and p=0.2; you applied a Bernouilli to the result. It is not clear to me that the result of the AddNoise function is equivalent to a Bernoulli. It just so happens that the result is close enough. On the other hand, I believe my AddNoise function is more accurate, but then I must find a way to transform the Nonlinear sampling into a soft function.

After some checking of the demos, I think I understand. If one follows a called to Nonlinear{sample} with a call to any distribution, whether Gaussian or Bernoulli, Forneylab will fit the data to whatever distribution that is specified. In the case of a Gaussian, a Laplace approximation will be used. In the case of a Bernoulli, some other approximation will be adopted. Did I get this right? Thanks again for all your help!

Below is the message that I sent before seeing your message.

I ran the code and generated the same results as in the Chapter. Now, I will generalize it. But I have a few comments and questions.

1) I substituted my AddNoise function (where I directly sample the Bernoulli's) and got the same result as yours. I do not understand why your approach works. Consider the following two lines from your version of the code:

AddNoise(hasSkill) = if hasSkill == 1 return 0.9 else return 0.2 end

# Now we are ready to build Figure 2.3
@RV s1 ~ Nonlinear{Sampling}(csharp,g=AddNoise,n_samples=50)
@RV isCorrect1 ~ Bernoulli(s1)

s1 is generated by sampling a distribution of just two numbers: 0.9 and 0.2. That is probably why you follow up with a Bernoulli distribution to generate isCorrect1. What is strange, is that if I replace your function by mine (where I sample within the function):

function addNoise(has_skill)
    if has_skill == 1.0
        res = sample(ProbabilityDistribution(Bernoulli, p=0.9))
        # f = x*0.9 + (1-x)*0.1
    else
        res = sample(ProbabilityDistribution(Bernoulli, p=0.2))
        # f = x*0.2 + (1-x)*0.8
    end
    return res
end

@RV s1 ~ Nonlinear{Sampling}(csharp,g=AddNoise,n_samples=50)
@RV isCorrect1 ~ Bernoulli(s1)

I get the same result as you do, even if I apply a Bernoulli to s1. I do not understand why this should be the case.

If I write instead:

@RV isCorrect1 ~ Nonlinear{Sampling}(csharp,g=AddNoise,n_samples=50)

I get the error:

ERROR: LoadError: MethodError: no method matching logPdf(::ProbabilityDistribution{Univariate, PointMass}, ::Float64)
Closest candidates are:
  logPdf(::ProbabilityDistribution{<:ForneyLab.VariateType, Function}, ::Any) at ~/.julia/packages/ForneyLab/y2b8v/src/probability_distribution.jl:122
  logPdf(::ProbabilityDistribution{Multivariate, GaussianMeanVariance}, ::Any) at ~/.julia/packages/ForneyLab/y2b8v/src/factor_nodes/gaussian_mean_variance.jl:70
  logPdf(::ProbabilityDistribution{Univariate, Gamma}, ::Any) at ~/.julia/packages/ForneyLab/y2b8v/src/factor_nodes/gamma.jl:57

My question is: why does this error occur? Why is your version the correct approach (applying a Bernoulli following sampling from AddNoise. ) I really do not understand the behavior of ForneyLab under these conditions. Thanks!

[erlebach]

I would like to learn more. In what module (file) is

@RV isCorrect1 ~ Bernoulli(s1)

computed from weighted samples? I cannot find anything relevant in Bernoulli.jl under FactorNodes folder. I just found the nonlinear_conjugate.jl module, from which I should be able to trace backward.

[semihakbayrak]

Hi @erlebach!

Let me start with your second paragraph. Yes, indeed, this is how Nonlinear{Sampling} works. In those cases where messages are non-Gaussian, ForneyLab runs importance sampling automatically. So in this model specification, ForneyLab take samples from prior Bernoulli distributions and evaluate them on backward messages coming from Nonlinear nodes to weigh samples. In the end, what you get as posterior distributions are weighted samples.

Regarding your first question, I'll do my best to explain why my model specification is the same with the one the authors provide in the book. Let us go through question 1. The generative model I defined works as the following: with probability 0.5 this person has the csharp skill. We reflect this information through a Bernoulli prior. Then, if this person possesses the skill, with 0.9 probability she will give the correct answer and with 0.1 probability she will give incorrect answer. This is equivalent to defining a Bernoulli likelihood with the parameter 0.9. Let us consider the second case, where this person does not have the csharp skill. In this case, she chooses an answer randomly and has a chance of 0.2 to get the correct answer. Then I can define this case with a Bernoulli distribution with parameter 0.2. So far, so good! Please note that what we have is a mixture of Bernoulli distributions. With 0.5 chance, the likelihood is Bernoulli(0.9) and with 0.5 chance the likelihood is Bernoulli(0.2). The way I define this mixture of Bernoulli distributions is through the deterministic node AddNoise. I define one Bernoulli likelihood Bernoulli(s1). And the value of s1 is defined through AddNoise function. This function says the following: if this person has the csharp skill, set s1 to 0.9, which means with 0.9 chance Bernoulli will generate the value 1 that stands for correct, and with 0.1 probability Bernoulli likelihood will generate 0 that stands for incorrect. On the contrary, if this person does not have the skill, then AddNoise will return 0.2 and set s1 to 0.2. This means that with 0.2 probability the observation will be 1. So if I set s1 to 0.2 and sample from Bernoulli(s1) let's say 10000 times, roughly 2000 of the times I'll get 1 that stands for the correct answer.

If you check figure 2.8 in the book, you'll see that the authors implicitly define these Bernoulli likelihood distributions. What we did is just to define these soft factors explicitly in our model specification and that's all.

I understand your way of thinking and motivation to use your AddNoise deterministic function as likelihood nodes. It makes total sense. I wish we could support this functionality but unfortunately we do not support it yet. Hope to address it in future versions.

[semihakbayrak]

I would like to learn more. In what module (file) is

@RV isCorrect1 ~ Bernoulli(s1)

computed from weighted samples? I cannot find anything relevant in Bernoulli.jl under FactorNodes folder. I just found the nonlinear_conjugate.jl module, from which I should be able to trace backward.

Bernoulli nodes send Beta backward messages to variables s1, s2, and s3. You can find it under /src/engines/julia/update_rules/bernoulli.jl. Importance sampling and hence weighted samples are executed around Nonlinear{Sampling} node. These rules are defined in /src/engines/julia/update_rules/nonlinear_sampling.jl and /src/factor_nodes/sample_list.jl.

I am not sure if it can be helpful but for visualization of automated importance sampling procedure of ForneyLab, you can check out Figure A.2 in this paper.

[erlebach]

Thank you for these details. This is all very interesting, and it also strengthens my Julia skills. I actually read about 50% of your paper a few days ago. I did not realize who I was communicating with! I might have more questions. Did you actually write the Nonlinear sampling section of ForneyLab?

I understand now. Your approach makes complete sense now! Thanks for the great explanation!

[semihakbayrak]

Great to hear that you find it interesting! Yes, Nonlinear{Sampling} node and related components like SampleList are coded by @ivan-bocharov @ThijsvdLaar and me. So if you have any further questions, please feel free to reach out.

[erlebach]

Hi @semihakbayrak !

I have refactored my code to use arrays, and now I am trying to implement figure 2.14, which has a cycle. You told me that ForneyLab does not handle loops, so what are my options? 1) use variational inference (which I have not yet done with ForneyLab), and 2) transform the graph to an acyclic graph. That is not a viable solution for large graphs. Here is my current code:

g = FactorGraph()
n_samples = 50

AddNoise(hasSkill) = if hasSkill == 1 return 0.9 else return 0.2 end
ANDfunc(skillx, skilly) = skillx * skilly
CompositeFunc(skillx, skilly) = AddNoise(ANDfunc(skillx, skilly))

nb_skills = 2
nb_questions = 4
p_skills = Vector{Variable}(undef, nb_skills) # one-hot coding
is_question_correct = Vector{Variable}(undef, nb_questions)
s = Vector{Variable}(undef, nb_questions)

for i in 1:nb_skills
    @RV p_skills[i] ~ Bernoulli(0.5)
end

@RV s[1] ~ Nonlinear{Sampling}(p_skills[1], g=AddNoise, n_samples=n_samples)
@RV s[2] ~ Nonlinear{Sampling}(p_skills[2], g=AddNoise, n_samples=n_samples)
@RV s[3] ~ Nonlinear{Sampling}(p_skills[1], p_skills[2], g=CompositeFunc, n_samples=n_samples)
@RV s[4] ~ Nonlinear{Sampling}(p_skills[1], p_skills[2], g=CompositeFunc, n_samples=n_samples)

for i in 1:nb_questions
    @RV is_question_correct[i] ~ Bernoulli(s[i])
    placeholder(is_question_correct[i], :isCorrect*i)  
end

# When does pfz have an effect? 
pfz = PosteriorFactorization(p_skills[1], p_skills[2], ids=[:C, :S]) 
dalgo = messagePassingAlgorithm([p_skills[1], p_skills[2]])
source_code = algorithmSourceCode(algo)
eval(Meta.parse(source_code));
print(source_code)

Is it just a matter of uncommenting the line pfz = ...? I doubt it. Would be too easy. :-)

As expected, my current error relates to the loop:

ArgumentError: The input graph contains a loop around Interface 3 (3) of Equality equ_p_skills_2_2
.

Stacktrace:
 [1] children(vertices::Vector{Interface}, graph::ForneyLab.DependencyGraph{Interface}; allow_cycles::Bool, breaker_sites::Set{Interface}, restrict_to::Set{Interface})
   @ ForneyLab ~/src/2022/ForneyLab.jl/src/dependency_graph.jl:151
 [2] summaryPropagationSchedule(variables::Vector{Variable}, clusters::Vector{ForneyLab.Cluster}; limit_set::Set{Edge}, target_sites::Vector{Interface}, breaker_sites::Vector{Interface})
   @ ForneyLab ~/src/2022/ForneyLab.jl/src/message_passing.jl:187
 [3] messagePassingSchedule(pf::PosteriorFactor)
   @ ForneyLab ~/src/2022/ForneyLab.jl/src/algorithms/posterior_factor.jl:55
 [4] messagePassingAlgorithm(target_variables::Vector{Variable}, pfz::PosteriorFactorization; ep_sites::Vector{Tuple}, id::Symbol, free_energy::Bool)
   @ ForneyLab ~/src/2022/ForneyLab.jl/src/algorithms/inference_algorithm.jl:82
 [5] messagePassingAlgorithm (repeats 2 times)
   @ ~/src/2022/ForneyLab.jl/src/algorithms/inference_algorithm.jl:58 [inlined]
 [6] top-level scope
   @ In[17]:30
 [7] eval
   @ ./boot.jl:373 [inlined]
 [8] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)
   @ Base ./loading.jl:1196

[semihakbayrak]

Hi @erlebach!

You're right that we have to break the loop and variational inference could be a good option. ForneyLab supports variational inference through message passing, however, it only allows us to break the loops around soft factors. In your example, the loop appears around deterministic factors, which is unfortunately something ForneyLab can't handle.

As @albertpod mentioned, the other probabilistic programming tool we are developing in our lab, ReactiveMP.jl supports LoopyBP. I'll discuss it with my friends to see whether we have the required implementations to execute importance sampling in ReactiveMP, and if this is the case then we can run it in a loopy manner to get an approximate solution.

I'll inform you after our discussion with my friends.

[erlebach]

Hi @semihakbayrak ,

Thank you, again. I wonder if there is a way to trick my model into having a loop around a soft factor? I do have one question: what do you mean by "Loop around a soft factor"? A loop contains both nodes and variables, and goes through them. From my perspective, a loop does not "go around" anything. So I must be misunderstanding something. Could you please explain? Also, is there an example of how to break a loop "around" a soft factor using ForneyLab? That would also be helpful. I can't imagine that a simple problem like mine cannot be solved by most frameworks, even though Infer.net (Microsoft), which I refuse to use, can handle this type of system from what I understand.

I tried ReactiveMP yesterday, have an error with the most simplistic model:

# AddNoise(hasSkill) = if hasSkill == 1 return 0.9 else return 0.2 end

@model function toy3(n)
    AddNoise(hasSkill) = if hasSkill == 1 return 0.9 else return 0.2 end

    y = datavar(Float64, n)  # 3 questions

    csharp ~ Bernoulli(0.5)
    sql ~ Bernoulli(0.5)

    isCorrect1 ~ Bernoulli(Addnoise(csharp))
    isCorrect2 ~ Bernoulli(Addnoise(sql))

    return y, isCorrect1, isCorrect2
end

[erlebach]

One last question: in order for loopy propatiin to work without breaking up the loops, all that would be needed is to provide initializer messages to start the iteration going. Is that currently possible? Also, ismit possible to trace the flow of information through ForneyLab? By that I mean getting a recordmof what routines get called and in what order? That would help with understanding and debugging. Macros are great but coming kicate the debugging process.

[semihakbayrak]

Hi @erlebach!

Here is an example loop that is caused by soft factors and variational inference helps us break the loop. If the loop contains a deterministic node then ForneyLab throws error. This is what I mean with loop around deterministic node. ReactiveMP does not have this issue with loops but unfortunately it does not support importance sampling in the level we need to be able to execute inference in Figure 2.14.

I see your point with message initialization but inference scheduling algorithms are quite involved in ForneyLab so there is no simple way to make ForneyLab work in Figure 2.14. Even if you initiate messages, ForneyLab scheduling algorithms will detect loops and won't be able to come up with inference schedules.

println(source_code) shows the details of the inference schedule. There you can see in what order messages are calculated. But it is not a debugging tool because if you get error in algorithm generation phase, you won't be able to call println(source_code).

[erlebach]

Assume a deterministic function g(x) = x^2. If I change this function so that instead, g(x) returns a sample of N(x^2, 0.01), for example, could I then use VMP on a loopy graph? Thanks.

SInce you have not replied yet, here is another question: In ForneyLab, does one apply VMP on all variables or none of them? For instance, if I have a few clamped variables, it appears necessary to allow for variational factors for these variables. I am checking whether I can adapt the demo on Variational Mixtures to my problem within the framework of ForneyLab.

By the way, is this forum the best place for these discussions? Or do you guys hang out anywhere else for in-depth discussions?

using ForneyLab

g = FactorGraph()

# Specify generative model
@RV _pi ~ Beta(1.0, 1.0)
# @RV m_1 ~ GaussianMeanVariance(_pi, 100.0)
# @RV w_1 ~ Gamma(0.01, 0.01)
# @RV m_2 ~ GaussianMeanVariance(_pi, 100.0)
# @RV w_2 ~ Gamma(0.01, 0.01)

# _pi = 0.5
@RV m_1 ~ Clamp(0.9)
@RV m_2 ~ Clamp(0.2)
@RV w_1 ~ Clamp(100)
@RV w_2 ~ Clamp(100)

z = Vector{Variable}(undef, n)
y = Vector{Variable}(undef, n)
for i = 1:n
    @RV z[i] ~ Bernoulli(_pi)
    @RV y[i] ~ GaussianMixture(z[i], m_1, w_1, m_2, w_2)

    placeholder(y[i], :y, index=i)
end

# Build the algorithm
# if I remove one of the clamped variables from pfz, the function fails.

pfz = PosteriorFactorization(_pi, m_1, w_1, m_2, w_2, z, ids=[:PI, :M1, :W1, :M2, :W2, :Z])
algo = messagePassingAlgorithm(free_energy=true)

# Generate source code
source_code = algorithmSourceCode(algo, free_energy=true);

[albertpod]

@erlebach this kind of question fits better in the Discussions.

In short, you need to specify posterior factorization because GaussianMixture node supports only variational updates.

We are sorry for the long answers. At the moment, all developers are busy. We will get back to your questions as soon as we can.

[erlebach]

Thanks. Where can I find the Discussions? I do not see this options on Github.

[albertpod]

Hi @erlebach!

Here you go [button next to PR]

[ThijsvdLaar]

Will close in favor of discussions forum

[erlebach]

Hi,

The "Go" button leads to a 404 error in my browser. Why that does not happen to you is a mystery. Perhaps I need special permissions on Github? Any help is appreciated. Thanks.

Gordon

On Wed, May 11, 2022 at 4:53 AM Albert Podusenko @.***> wrote:

Hi @erlebach https://github.com/erlebach!

Here you go https://github.com/biaslab/ForneyLab.jl/discussions/landing [button next to PR]

— Reply to this email directly, view it on GitHub https://github.com/biaslab/ForneyLab.jl/issues/200#issuecomment-1123385687, or unsubscribe https://github.com/notifications/unsubscribe-auth/AACPIZEI4KFYHQPIUG2RZNDVJNYRJANCNFSM5UJYPVRA . You are receiving this because you were mentioned.Message ID: @.***>

[albertpod]

Hi @erlebach,

Indeed, the link was broken (fixed now). Check the attached screenshot.

Regards, Albert