Graphical and Probabilistic Models

What Are Probabilistic Graphical Models?

Probabilistic graphical models (PGMs) represent joint distributions over many random variables using graphs. Nodes are random variables; edges encode conditional (in)dependence. PGMs offer:

• Compact representation: factor the joint distribution exploiting independence
• Efficient inference: compute marginals and conditionals without full enumeration
• Learning: estimate parameters and structure from data
• Interpretability: explicit representation of variable relationships

Two main families:

Bayesian Networks

Structure

A DAG where each node X has a Conditional Probability Distribution (CPD): P(X | Parents(X)).

The joint distribution factorizes as:

P(X₁, …, Xₙ) = Πᵢ P(Xᵢ | Parents(Xᵢ))

Example — Simple diagnostic network:

Smoking → Lung Cancer → Dyspnea
                     ↘
Tuberculosis ————————→ XRay

Each node stores a CPT (Conditional Probability Table): - P(Smoking) — marginal (no parents) - P(LungCancer | Smoking) — 2×2 table - P(Dyspnea | LungCancer, Tuberculosis) — 4×2 table

d-Separation

Three patterns determine conditional independence via d-separation:

Two nodes are d-separated if every path between them is blocked. D-separation ↔ conditional independence.

Inference in Bayesian Networks

Variable elimination: sum out unwanted variables in a specific order:

P(Y | E=e) ∝ Σ_{X \ {Y}} P(X₁,…,Xₙ, E=e)

Apply chain rule, then eliminate variables by summing over them left to right.

Belief propagation: pass messages between nodes on tree-structured graphs. Exact inference for trees; approximate (loopy BP) for general graphs.

from pgmpy.models import BayesianNetwork
from pgmpy.factors.discrete import TabularCPD
from pgmpy.inference import VariableElimination

# Define structure
model = BayesianNetwork([("Smoking","LungCancer"),
                          ("LungCancer","Dyspnea"),
                          ("Tuberculosis","Dyspnea")])

# Define CPTs
cpd_smoke = TabularCPD("Smoking", 2, [[0.7],[0.3]])
cpd_lung  = TabularCPD("LungCancer", 2,
                        [[0.99, 0.90],[0.01, 0.10]],
                        evidence=["Smoking"], evidence_card=[2])
cpd_tb    = TabularCPD("Tuberculosis", 2, [[0.99],[0.01]])
cpd_dys   = TabularCPD("Dyspnea", 2,
                        [[0.9,0.3,0.2,0.1],[0.1,0.7,0.8,0.9]],
                        evidence=["LungCancer","Tuberculosis"], evidence_card=[2,2])

model.add_cpds(cpd_smoke, cpd_lung, cpd_tb, cpd_dys)
model.check_model()

infer = VariableElimination(model)
result = infer.query(["LungCancer"], evidence={"Dyspnea": 1, "Smoking": 1})
print(result)

Learning Bayesian Networks

Parameter learning (structure known): count-based MLE or Bayesian estimation (Dirichlet prior).

Structure learning: find the DAG that best explains data. - Score-based: search DAG space to maximize BIC/BDe score - Constraint-based: run independence tests (PC algorithm) - Hybrid: combine both (MMHC)

from pgmpy.estimators import MaximumLikelihoodEstimator, BDeuScore, HillClimbSearch

# Parameter learning
model.fit(data, estimator=MaximumLikelihoodEstimator)

# Structure learning
hc = HillClimbSearch(data)
best_model = hc.estimate(scoring_method=BDeuScore(data))

Markov Random Fields (MRFs)

Structure

Undirected graph where edges indicate direct dependency. The joint distribution factorizes over cliques (fully-connected subgraphs):

P(X₁, …, Xₙ) = (1/Z) Πc ψc(Xc)

where ψc are potential functions and Z is the partition function (normalization constant — computing Z is often intractable).

Applications: image segmentation (pixels connected to neighbors), social networks, physics (Ising model), error-correcting codes.

Ising Model (Binary MRF on a Grid)

P(x) ∝ exp(β Σ_{(i,j)∈E} xᵢxⱼ + h Σᵢ xᵢ)

where xᵢ ∈ {−1, +1}. β controls neighbor coupling (ferromagnetic when β>0); h is external field.

Conditional Random Fields (CRF)

CRFs are discriminative MRFs that model P(Y | X) directly — avoids modeling P(X).

Standard linear-chain CRF for sequence labeling:

P(y | x) = (1/Z(x)) exp(Σt Σk θk fk(yₜ₋₁, yₜ, x, t))

Feature functions fk capture label-label transitions and label-observation interactions.

CRFs are the classical model for NER, POS tagging, and semantic segmentation (dense CRF for pixel labeling).

import sklearn_crfsuite

crf = sklearn_crfsuite.CRF(algorithm="lbfgs", c1=0.1, c2=0.1, max_iterations=100)
crf.fit(X_train, y_train)  # X_train: list of sentences; each sentence = list of feature dicts

Gaussian Graphical Models

Special case: all variables are jointly Gaussian. The precision matrix Ω = Σ⁻¹ encodes the conditional independence structure:

Ωᵢⱼ = 0 ↔ Xᵢ ⊥ Xⱼ | all other variables

Graphical Lasso (glasso): learn sparse Ω by L1-regularized MLE:

max log det Ω − tr(SΩ) − λ‖Ω‖₁

from sklearn.covariance import GraphicalLasso

model = GraphicalLasso(alpha=0.1, max_iter=100)
model.fit(X)
precision_matrix = model.precision_   # zero entries = no edge
covariance_matrix = model.covariance_

Applications: gene regulatory networks, financial correlations, neuroimaging connectivity.

Hidden Markov Models (HMM)

A generative model for sequences with hidden states:

- Hidden state sequence: z₁, z₂, …, zₜ (Markov chain) P(zₜ | zₜ₋₁) — transition matrix A - Observation sequence: x₁, x₂, …, xₜ P(xₜ | zₜ) — emission matrix B - Initial distribution: π (prior over z₁)

z₁ → z₂ → z₃ → … (hidden)
↓    ↓    ↓
x₁   x₂   x₃    (observed)

Three fundamental HMM problems:

from hmmlearn import hmm
import numpy as np

# Discrete HMM (Multinomial emissions)
model = hmm.MultinomialHMM(n_components=3, n_iter=100, tol=0.01)
model.fit(X_sequences, lengths=sequence_lengths)

# Decode most likely states
states = model.predict(X_new, lengths=[len(X_new)])

# Continuous (Gaussian) HMM
gauss_hmm = hmm.GaussianHMM(n_components=4, covariance_type="full", n_iter=100)
gauss_hmm.fit(X_continuous)

Applications: speech recognition, genomic sequence analysis, financial regime detection, gesture recognition.

Latent Dirichlet Allocation (LDA)

A hierarchical Bayesian model for topic modeling in text documents:

• Each document is a mixture of topics: θd ~ Dirichlet(α)
• Each topic is a distribution over words: φk ~ Dirichlet(β)
• Each word is drawn by: choose topic z ~ θd, then word w ~ φz

from sklearn.decomposition import LatentDirichletAllocation
from sklearn.feature_extraction.text import CountVectorizer

vectorizer = CountVectorizer(max_features=5000, min_df=5, stop_words="english")
X_counts = vectorizer.fit_transform(documents)

lda = LatentDirichletAllocation(n_components=20, max_iter=20,
                                  learning_method="online",
                                  random_state=42)
lda.fit(X_counts)

# Show top words per topic
feature_names = vectorizer.get_feature_names_out()
for topic_idx, topic in enumerate(lda.components_):
    top_words = [feature_names[i] for i in topic.argsort()[-10:]]
    print(f"Topic {topic_idx}: {' '.join(top_words)}")

# Document-topic distribution
doc_topics = lda.transform(X_counts)  # shape (n_docs, n_topics)

Variational Inference

When exact inference (computing P(latent | data)) is intractable, approximate with a simpler distribution:

ELBO = E_q[log P(x, z)] − E_q[log q(z)] = log P(x) − KL(q(z) ‖ P(z|x))

Maximize ELBO to find the best approximation q(z) to the true posterior P(z|x). This is the foundation of: - Variational Autoencoders (VAEs) - Mean-field variational inference - Black-box VI (BBVI)

MCMC (Markov Chain Monte Carlo)

Sample from complex posteriors by constructing a Markov chain that has the target distribution as its stationary distribution.

Metropolis-Hastings: 1. Propose a new state x' ~ q(x' | x) 2. Accept with probability min(1, P(x')q(x|x') / P(x)q(x'|x))

Gibbs Sampling: Cycle through variables, sampling each from its conditional P(xᵢ | x_{-i}).

Hamiltonian Monte Carlo (HMC) / NUTS: Uses gradient information to make efficient proposals in high dimensions. Used in Stan, PyMC.

import pymc as pm
import numpy as np

with pm.Model() as bayesian_regression:
    # Priors
    alpha = pm.Normal("alpha", mu=0, sigma=10)
    beta  = pm.Normal("beta",  mu=0, sigma=10, shape=n_features)
    sigma = pm.HalfNormal("sigma", sigma=1)

    # Likelihood
    mu = alpha + pm.math.dot(X_train, beta)
    y_obs = pm.Normal("y_obs", mu=mu, sigma=sigma, observed=y_train)

    # Sample
    trace = pm.sample(2000, tune=1000, chains=4, return_inferencedata=True)

pm.plot_posterior(trace, var_names=["alpha", "beta"])

Summary: When to Use Which Model