Probabilistic Inference

Overview

The structures of sds.probabilistic store a model: its variables, its graph and its factors. They answer no question about it beyond the probability of one complete assignment (joint()). This subpackage holds the inference algorithms that do:

Function

Model

Answers

Exact?

variable_elimination()

Bayesian network, Markov random field

\(P(Q \mid e)\) for one or more query variables

yes

belief_propagation()

Bayesian network, Markov random field

\(P(X \mid e)\) for every variable at once

on tree-shaped models

forward()

hidden Markov model

\(\log P(o_{1:T})\) and \(P(S_t \mid o_{1:t})\)

yes

viterbi()

hidden Markov model

\(\arg\max_{s_{1:T}} P(s_{1:T}, o_{1:T})\)

yes

Note

Factor offers no product or marginalization. Each algorithm copies the model’s factors into private working tables through scope and value(), computes on those copies, and returns plain Factor objects. The model is never modified.

Throughout this page, the running examples are two classic models: Pearl’s burglary alarm network and Russell and Norvig’s umbrella world.

        graph TD
    B((Burglary)) --> A((Alarm))
    E((Earthquake)) --> A
    A --> J((JohnCalls))
    A --> M((MaryCalls))

Variable Elimination

The posterior of a query given evidence is a sum of the factor product over every hidden variable. Written naively, that sum enumerates every joint assignment, exponentially many. Variable elimination distributes each sum over the product and eliminates the hidden variables one at a time: multiply the tables that mention the variable, sum it out, put the result back.

For the alarm network and the query \(P(B \mid j, m)\), eliminating \(E\) then \(A\) gives

\[P(B \mid j, m) \;\propto\; P(B) \sum_{a} P(j \mid a)\, P(m \mid a) \sum_{e} P(e)\, P(a \mid B, e)\]

and \(P(B = \mathrm{true} \mid j, m) \approx 0.284\): even when both neighbours call, a burglary remains less likely than a false alarm.

The cost is exponential in the largest intermediate table, which depends on the elimination order. By default the order is chosen greedily, always eliminating the variable that builds the smallest table.

from sds.algorithms.probabilistic_algorithms import variable_elimination

posterior = variable_elimination(alarm, [burglary],
                                 {john_calls: "true", mary_calls: "true"})
posterior.value({burglary: "true"})      # 0.284...
sds.algorithms.probabilistic_algorithms.variable_elimination(model, query, evidence=None, order=None)[source]

Compute the exact distribution P(query | evidence).

Parameters:

  • model (AbstractGraphicalModel) – A Bayesian network or a Markov random field. Its factors are read through factors() and never modified.

  • query (sequence of RandomVariable) – One or more variables of the model whose joint distribution is wanted. They must not be observed in evidence.

  • evidence (mapping of RandomVariable to str, optional) – Observed states. Defaults to no evidence.

  • order (sequence of RandomVariable, optional) – Order in which to eliminate the hidden variables (every variable that is neither queried nor observed). It must list each hidden variable exactly once. By default a greedy min-size order is used: at each step, eliminate the variable whose elimination builds the smallest table, ties broken by the model’s variable order.

Returns:

A factor over query (in the given order) whose values are normalized to sum to 1.

Return type:

Factor

Raises:

InvalidParameterError – If query is empty, contains duplicates, a variable outside the model, or an observed variable; if the evidence names an unknown variable or state; if order does not list exactly the hidden variables; or if the evidence has probability zero.

Examples

>>> from sds.probabilistic import BayesianNetwork, Factor, RandomVariable
>>> rain = RandomVariable("Rain", ("yes", "no"))
>>> wet = RandomVariable("Wet", ("yes", "no"))
>>> bn = BayesianNetwork()
>>> for v in (rain, wet):
...     bn.add_variable(v)
>>> bn.add_edge(rain, wet)
>>> bn.set_cpt(rain, Factor((rain,), {("yes",): 0.2, ("no",): 0.8}))
>>> bn.set_cpt(wet, Factor((wet, rain), {
...     ("yes", "yes"): 0.9, ("no", "yes"): 0.1,
...     ("yes", "no"): 0.1, ("no", "no"): 0.9}))
>>> posterior = variable_elimination(bn, [rain], {wet: "yes"})
>>> round(posterior.value({rain: "yes"}), 4)
0.6923

Notes

Time and space complexity: O(n · d^w), with n the number of variables, d the largest domain size and w the size of the largest table built during elimination (the induced width of the chosen order, plus one). On a tree-structured model w is at most the largest factor scope.

The model’s factors are copied into working tables, and all products and sums happen on those copies: Factor itself offers no algebra (#84).

A variable that no factor mentions contributes a constant and is treated as uniform. For a Bayesian network, every CPT must therefore be set before inference.

Belief Propagation

Belief propagation computes every marginal at once by passing messages along the model’s factor graph, which links each factor to the variables of its scope. Each variable tells each neighbouring factor what the rest of the graph thinks of it, and each factor answers by summing its table against the other messages.

        graph LR
    B((B)) --- fA[P A given B,E]
    E((E)) --- fA
    fA --- A((A))
    A --- fJ[P J given A]
    A --- fM[P M given A]
    fJ --- J((J))
    fM --- M((M))
    B --- fB[P B]
    E --- fE[P E]

The alarm network’s factor graph above is a tree: messages settle after a few rounds and the marginals are exact, the same as variable elimination.

Warning

When the factor graph contains a cycle, the same procedure is loopy belief propagation. It frequently converges to useful approximations, but not to the exact answer. On the sprinkler network (Cloudy → Sprinkler, Cloudy → Rain, both → WetGrass), it estimates \(P(\mathrm{Rain} \mid \mathrm{WetGrass}) \approx 0.78\) where the exact value is about 0.71: the evidence travels around the cycle and is counted twice. If the messages do not settle within max_iterations, a RuntimeWarning is emitted.

sds.algorithms.probabilistic_algorithms.belief_propagation(model, evidence=None, *, max_iterations=100, tolerance=1e-09)[source]

Compute the marginal P(X | evidence) of every variable.

Parameters:

  • model (AbstractGraphicalModel) – A Bayesian network or a Markov random field, read through variables() and factors() only.

  • evidence (mapping of RandomVariable to str, optional) – Observed states. An observed variable’s marginal puts all its mass on the observed state.

  • max_iterations (int, optional) – Upper bound on the number of message-passing rounds (default 100).

  • tolerance (float, optional) – Convergence threshold on the largest change of any message between two rounds (default 1e-9).

Returns:

For every variable of the model, a one-variable factor holding its normalized marginal, in the model’s variable order.

Return type:

dict of RandomVariable to Factor

Raises:

InvalidParameterError – If the evidence names an unknown variable or state, if max_iterations is not positive, or if the evidence has probability zero.

Warns:

RuntimeWarning – If the messages have not converged after max_iterations rounds (possible only on a factor graph with cycles). The marginals returned are then those of the last round.

Examples

>>> from sds.probabilistic import BayesianNetwork, Factor, RandomVariable
>>> rain = RandomVariable("Rain", ("yes", "no"))
>>> wet = RandomVariable("Wet", ("yes", "no"))
>>> bn = BayesianNetwork()
>>> for v in (rain, wet):
...     bn.add_variable(v)
>>> bn.add_edge(rain, wet)
>>> bn.set_cpt(rain, Factor((rain,), {("yes",): 0.2, ("no",): 0.8}))
>>> bn.set_cpt(wet, Factor((wet, rain), {
...     ("yes", "yes"): 0.9, ("no", "yes"): 0.1,
...     ("yes", "no"): 0.1, ("no", "no"): 0.9}))
>>> marginals = belief_propagation(bn, {wet: "yes"})
>>> round(marginals[rain].value({rain: "yes"}), 4)
0.6923

Notes

Time complexity: \(O(I \sum_f |f| \cdot |\mathrm{scope}(f)|)\) for \(I\) rounds, with \(|f|\) the number of entries of factor \(f\). Space complexity: one message of the variable’s domain size per factor-variable pair, in each direction.

Evidence is applied by clamping: each observed variable receives an extra unary factor that is 1 on the observed state and 0 elsewhere. Messages are normalized after every update, which changes no marginal and keeps the numbers in a safe range.

Hidden Markov Models

        graph LR
    S1((S₁)) --> S2((S₂)) --> S3((S₃))
    S1 --> O1[O₁]
    S2 --> O2[O₂]
    S3 --> O3[O₃]

A hidden Markov model is a chain of hidden states, each emitting an observation. Two questions dominate its use: how likely is this observation sequence? (evaluation, the forward algorithm) and which state sequence best explains it? (decoding, the Viterbi algorithm). Both run in \(O(T\,|S|^2)\) instead of enumerating the \(|S|^T\) state sequences, by reusing, at each step, a quantity computed for the previous one.

Forward

The forward recursion carries, for each state, the probability of the observations so far ending in that state. It is rescaled at every step, which yields the filtered distributions \(P(S_t \mid o_{1:t})\) and keeps long sequences from underflowing; the likelihood is returned as a logarithm.

from sds.algorithms.probabilistic_algorithms import forward

log_p, filtered = forward(umbrella_hmm, ["yes", "yes"])
filtered[-1].value({weather: "rain"})    # 0.883
sds.algorithms.probabilistic_algorithms.forward(hmm, observations)[source]

Compute the log-likelihood of observations and the filtered states.

Parameters:

  • hmm (AbstractHiddenMarkovModel) – A fully configured hidden Markov model.

  • observations (sequence of str) – The observed symbols o_1 ... o_T, each in the observation domain.

Returns:

  • log_likelihood (float) – Natural logarithm of P(o_1, ..., o_T).

  • filtered (list of Factor) – For each time step t, a factor over the hidden-state variable holding P(S_t | o_1, ..., o_t).

Raises:

InvalidParameterError – If the sequence is empty, contains an unknown symbol, if the model is not fully configured, or if the sequence has probability zero under the model.

Return type:

Tuple[float, List[Factor]]

Examples

>>> from sds.probabilistic import Factor, HiddenMarkovModel, RandomVariable
>>> weather = RandomVariable("Weather", ("rain", "sun"))
>>> umbrella = RandomVariable("Umbrella", ("yes", "no"))
>>> hmm = HiddenMarkovModel()
>>> hmm.set_states(weather)
>>> hmm.set_observations(umbrella)
>>> prev = hmm.previous_states()
>>> hmm.set_initial_distribution(
...     Factor((weather,), {("rain",): 0.5, ("sun",): 0.5}))
>>> hmm.set_transition_model(Factor((weather, prev), {
...     ("rain", "rain"): 0.7, ("sun", "rain"): 0.3,
...     ("rain", "sun"): 0.3, ("sun", "sun"): 0.7}))
>>> hmm.set_emission_model(Factor((umbrella, weather), {
...     ("yes", "rain"): 0.9, ("no", "rain"): 0.1,
...     ("yes", "sun"): 0.2, ("no", "sun"): 0.8}))
>>> log_p, filtered = forward(hmm, ["yes", "yes"])
>>> round(filtered[-1].value({weather: "rain"}), 3)
0.883

Notes

Time complexity: \(O(T |S|^2)\). Space complexity: \(O(T |S|)\) for the filtered distributions.

The likelihood itself is math.exp(log_likelihood), which may underflow to 0.0 for long sequences even though the log stays finite.

Viterbi

Viterbi replaces the sum of the forward recursion by a maximum and keeps a back-pointer to the best predecessor of each state; following the pointers back from the best final state rebuilds the path. It works on logarithms, so a zero probability becomes \(-\infty\) and long sequences are safe.

from sds.algorithms.probabilistic_algorithms import viterbi

path, log_p = viterbi(umbrella_hmm, ["yes", "yes", "no", "yes", "yes"])
# ['rain', 'rain', 'sun', 'rain', 'rain']
sds.algorithms.probabilistic_algorithms.viterbi(hmm, observations)[source]

Return the most likely hidden-state path and its log-probability.

Parameters:

  • hmm (AbstractHiddenMarkovModel) – A fully configured hidden Markov model.

  • observations (sequence of str) – The observed symbols o_1 ... o_T, each in the observation domain.

Returns:

  • path (list of str) – The most likely states s_1 ... s_T. When several paths tie, the one preferring earlier states of the domain is returned.

  • log_probability (float) – Natural logarithm of the joint probability P(s_1 ... s_T, o_1 ... o_T) of that path with the observations, the quantity hmm.joint(path, observations) returns.

Raises:

InvalidParameterError – If the sequence is empty, contains an unknown symbol, if the model is not fully configured, or if every state path has probability zero.

Return type:

Tuple[List[str], float]

Examples

>>> from sds.probabilistic import Factor, HiddenMarkovModel, RandomVariable
>>> weather = RandomVariable("Weather", ("rain", "sun"))
>>> umbrella = RandomVariable("Umbrella", ("yes", "no"))
>>> hmm = HiddenMarkovModel()
>>> hmm.set_states(weather)
>>> hmm.set_observations(umbrella)
>>> prev = hmm.previous_states()
>>> hmm.set_initial_distribution(
...     Factor((weather,), {("rain",): 0.5, ("sun",): 0.5}))
>>> hmm.set_transition_model(Factor((weather, prev), {
...     ("rain", "rain"): 0.7, ("sun", "rain"): 0.3,
...     ("rain", "sun"): 0.3, ("sun", "sun"): 0.7}))
>>> hmm.set_emission_model(Factor((umbrella, weather), {
...     ("yes", "rain"): 0.9, ("no", "rain"): 0.1,
...     ("yes", "sun"): 0.2, ("no", "sun"): 0.8}))
>>> path, log_p = viterbi(hmm, ["yes", "yes", "no", "yes", "yes"])
>>> path
['rain', 'rain', 'sun', 'rain', 'rain']

Notes

Time complexity: \(O(T |S|^2)\). Space complexity: \(O(T |S|)\) for the back-pointers.

The most likely path is not the sequence of individually most likely states: filtering or smoothing each step separately can yield a sequence that is improbable, or even impossible, as a whole.

Best Practices

✅ Use variable_elimination for exact answers to a few specific queries; use belief_propagation when every marginal is needed and the model is tree-shaped.

✅ Set every CPT of a Bayesian network before inference: a variable no factor mentions is treated as uniform.

✅ Compare HMM sequences by log-likelihood; math.exp of a long sequence’s log-likelihood underflows to 0.0.

❌ Do not read a loopy belief propagation result as exact.

❌ Do not reconstruct the most likely path from the most likely state at each step: the individually best states can form an impossible sequence. Use viterbi.

References