How to properly calculate inference for Bayesian net?

[ferencbartok]

I'm trying to do inference like this:

#model is BayesianNetwork
inference = VariableElimination(model)
        for node in model.nodes():
            if node in observations:
                query = inference.query(variables=[node])
                infer_probabilities[node] = query.values
            else:
                query = inference.query(variables=[node], evidence=observations)
                infer_probabilities[node] = query.values

Is there an easier/better way?
If I'd like to calculate marginals without any evidence I use this:
query = inference.query(variables=model.nodes())
Then I get: Unable to allocate 32.0 GiB for an array with shape...

The problem with the first one is that it is pretty slow compared to the existing (smile) solution and I don't understand why. (We'd like to migrate to open-source from SMILE.)
If I use BeliefPropagation then I get: invalid value encountered in scalar divide S_by_C = {key: S[key] / C[key] for key in S}

I'd appreciate any pointers, thank you!

[ankurankan]

@ferencbartok I am not exactly sure what you are trying to compute in the example above. Could you elaborate a bit please? If you are looking to compute the (posterior) marginal distributions over all variables given some evidence you could do something like:

In [13]: from pgmpy.utils import get_example_model
In [14]: model = get_example_model('alarm')
In [15]: from pgmpy.inference import VariableElimination
In [16]: infer = VariableElimination(model)
In [18]: observation = {'CVP': 'LOW', 'PAP': 'HIGH'}
In [19]: for node in model.nodes():
    ...:     if node not in observation:
    ...:         infer.query([node], evidence=observation)

By default pgmpy tries to compute the joint disribution over all the variables in the variables argument, that's why it throws an error when you try to do query = inference.query(variables=model.nodes()) because the joint distribution over all the variables would be too big to fit in memory. You could pass joint=False to the query method to get marginal distributions. However, if you pass too many variables, it can still result in memory error because it first computes the joint distribution.

I imagine that SMILE is faster in this case because it might be using an approximate Belief Propagation algorithm (and is written in C++), which computes the marginals over all the nodes together. Variable Elimination isn't designed to do that and our Belief Propagation algorithm currently doesn't have the functionality to access these marginals right now. If you could explain a bit what you are trying to compute, I can try to modify the Belief Propagation algorithm to do that.

[ferencbartok]

Hm, I used this:

inference = ApproxInference(model)
result = inference.query(variables=list(model.nodes()), joint=False, n_samples=500)
for node in model.nodes():
    infer_probabilities[node] = result[node].values

But if I try as you suggested I still get the wrong numbers. Also, it's much slower compared to querying all the nodes at once.

inference_approx = ApproxInference(model)
for node in model.nodes():
    infer_probabilities[node] = inference_approx.query([node], n_samples=100).values

(I'm not using it in your example as it of course should return the same but my own dataset.)

As it turns out the numbers seem okay except for a "few", there it seems like the order is reversed for the approximated inference vs. exact inference probabilities. But probably I'm doing something wrong. Unfortunately BeliefPropagation is much slower even with 100 samples only and not giving good enough results. I'll try and make better measurements for exact inference to see if pgmpy is a viable option regarding speed. Otherwise, it's a really nice library and fairly easy to code.
Thank you for your help!

[ferencbartok]

Just as a follow-up question:

Currently our most time-consuming calculation is this:

• do initial inference for all nodes/variables with 0 evidence (I use node and variable interchangeably)
• iterate over ~40-50-60% of our nodes (let's call them "observables"; size depends on the given network)
  • for every observable node set one of the states as evidence (do it for all possible states, which number can usually be around 2-3-4)
  • calculate inference for the other 40-50-60% of nodes*
  • calculate a value for the observable node which shows how much effect/impact changing it's state has on the rest of the variables (~comparing original inference with 0 evidence to the just now calculated probabilities)

This solution using any library is wasteful as most of the things we recalculate should be the same as before. I'm wondering if you could suggest some way to do this calculation faster using pgmpy?

[ankurankan]

I think the sampling based Approximate Inference could work in this case.

do initial inference for all nodes/variables with 0 evidence (I use node and variable interchangeably)

Using the sampling appraoch, the samples can be generated once (which is the most expensive step) and then the marginals can be computed using the same samples.

from pgmpy.utils import get_example_model
from pgmpy.inference import ApproxInference

model = get_example_model('alarm')
samples = model.simulate(int(1e5))

infer = ApproxInference(model)
for node in model.nodes():
    print(infer.query([node], samples=samples).values)

iterate over ~40-50-60% of our nodes (let's call them "observables"; size depends on the given network)

Similar to last case, data can be simulated under some evidence (i.e. observables) and that can be used to quickly calculate the marginals over the rest of the variables.

observable = {'BP': 'LOW', 'HISTORY': 'TRUE'}
samples = model.simulate(int(1e5), evidence=observable)

infer = ApproxInference(model)
for node in model.nodes():
    if node not in observable.keys():
        print(infer.query([node], samples=samples).values)

The other option as I had mentioned before would be to use Belief Propagation as that allows us to quickly compute the marginals but might need to wait till the refactor is done.

[ferencbartok]

Thank you for the reply!

I've tried ApproxInference again and I have the same strange problem as before, some values are reversed, some are not. What I mean by this:
query.values = array([0.75305, 0.24695]), then I check the state_name in the expected file and they are the ~same.
But for another I get array([0.71448, 0.28552]), but the expected should be the reversed like: array([0.28552, 0.71448]).
Unfortunately I have no idea why some are in good order and some not.

My code is:

infer_probabilities: InferProbabilities = {}

    if observations:
        samples = model.simulate(int(1e5), evidence=observations)
    else:
        samples = model.simulate(int(1e5))
    inference = ApproxInference(model)

    for node in model.nodes():
        if node in target_node_ids:
              query = inference.query(variables=[node], samples=samples)
              infer_probabilities[node] = query.values.tolist()

Do you have any idea why this could be happening?

Otherwise, I'm not sure why this approach would help in this use-case, we are not reusing the created samples between inferences. I do see the point that compared to earlier ApproxInference solution* now we just create samples once for each inference compared to once for each query. But this could have been done regardless of the use_case/algorithm, no?
(*) for node in... infer_probabilities[node] = inference_approx.query([node], n_samples=100).values

Maybe we could do that we calculate the samples only once at the beginning and then use that for all the queries in all the inferences +setting the extra evidence in the current iteration? So I mean something like this:

query = inference.query(variables=[node], samples=samples, evidence=the_current_1_node_1_state_evidence)

Unfortunately I can't test this because of the strange problem with ApproxInference.

[ankurankan]

@ferencbartok For the first issue of reversed values, I am fairly certain that it's happening because the CPDs have different state orders. You can verify that by printing the state order (query.state_names[query_var]) from both inferences. Alternatively, you can use the get_value method of the CPD to get the probability corresponding to a specific state, instead of comparing the values array.

In [12]: from pgmpy.utils import get_example_model

In [13]: model = get_example_model('alarm')
 
In [14]: from pgmpy.inference import VariableElimination

In [15]: infer = VariableElimination(model)

In [16]: query = infer.query(['HISTORY'])

In [17]: for state in model.states['HISTORY']:
    ...:     print(query.get_value(HISTORY=state))
    ...: 
tensor(0.0545, dtype=torch.float64)
tensor(0.9455, dtype=torch.float64)

In [18]: for state in model.states['HISTORY']:
    ...:     print(query.get_value(**{'HISTORY': state}))
    ...: 
tensor(0.0545, dtype=torch.float64)
tensor(0.9455, dtype=torch.float64)

Regarding the second point, when n_samples is specified in ApproxInference.query, it automatically generates those samples. But when samples is specified, it uses those already generated samples to compute the probabilities, but the samples should already be generated under evidence, if any. Therefore,

1. For computing marginals without evidence, we can generate samples once and compute all the marginals. Something like this would give the marginals over all variables:

samples = model.simulate(int(1e5))
infer.query(list(model.nodes()), samples=samples, joint=False)

2. When some evidence is specified, we would need to generate samples under that specified evidence. Something like this would give the marginal posterior distribution over all variables given some observation.

samples = model.simulate(int(1e5), evidence=observations)
infer.query(list(set(model.nodes()) - set(observations.keys()), samples=samples, joint=False)

Lastly, another potential way to do this is to generate samples just once and do rejection sampling to get samples under evidence. This approach only generates samples once but has the drawback that if the probability of a given evidence is low, the results can be quite inaccurate.

samples = model.simulate(int(1e6))
infer = ApproxInference(model)

# Compute the marginals as before
marginals = infer.query(list(model.nodes()), samples=samples, joint=False)

# Iterate over the given list of observations
observations = [{'HISTORY': 'FALSE', 'BP': 'LOW'}, {'HISTORY': 'TRUE', 'BP': 'HIGH'}, ... , ]
for observation in observations:
      samples_evid = samples.query(' and '.join([f'{k} == {repr(v)}' for k, v in observation.items()]))
      infer.query(list(set(model.nodes()) - set(observation.keys()), samples=samples_evid, joint=False)

[ferencbartok]

Thank you!

About the reversed issue: I've checked a few of my tests based on your code and the probabilities given the exact states are correct. But why is this happening when I set up everything in the same way, the only difference is that I use ApproxInference instead of VariableElimination? Everything else is the same.

About the second question:
Yes, I understand that now we are creating samples only for each inference compared to the earlier solution where samples were generated for every variable query. However, it would still be too slow because generating samples for each inference takes a lot of time. Maybe using smaller samplesize for the "only 1 evidence" use-cases could speed things up, I'll try that later.

As for the rejection sampling approach, this seems promising, but I have 2 problems:

• I can't reach/see the query.values, so the actual posterior probabilities
• I think there is something wrong in this case with variable names that contain '-' character:

observations = [{'e5aa1768-1bc1-45ca-bf21-9571d55d75f7': 0},
                    {'e5aa1768-1bc1-45ca-bf21-9571d55d75f7': 1}]
e5aa1768 -1 bc1 -45 ca -bf21 -9571 d55d75f7 ==0 
                ^^^
SyntaxError: invalid syntax

This works:
observations = [{'a': '0', 'b': '1'}, {'a': '1', 'b': '0'}]
Any ideas why this could be happening?