Sesión 13

Sesión 13#

Aprendizaje de estructura en redes bayesianas#

Objetivos:

  • Comprender el problema de aprendizaje de estructura en redes bayesianas.

  • Entender la puntuación de máxima verosimilitud y sus limitaciones.

  • Comparar las puntuaciones AIC, BIC y BD, y sus propiedades.

Referencias:

  • Probabilistic Graphical Models: Principles and Techniques, By Daphne Koller and Nir Friedman. Cap. 18.

  • Mastering Probabilistic Graphical Models Using Python, By Ankur Ankan and Abinash Panda. Cap. 5.

  • Probabilistic Graphical Models Specialization, Coursera. Prof. Daphne Koller.

  • CS228 Lecture Notes — Structure Learning. Ermon Group, Stanford. https://ermongroup.github.io/cs228-notes/learning/structure/

1. Motivación#

En las sesiones anteriores aprendimos a estimar los parámetros de una red bayesiana, asumiendo que la estructura ya era conocida.
Sin embargo, en muchos problemas reales no contamos con esa estructura de antemano.

El aprendizaje de estructura busca responder:

Dado un conjunto de datos \(\mathcal{D}\), ¿cuál es la red bayesiana más adecuada para modelar las dependencias entre las variables?

1.1. ¿Por qué aprender la estructura?#

  • Modelado sin expertos: cuando el conocimiento del dominio no alcanza para definir la estructura manualmente.

  • Descubrimiento de conocimiento: cuando identificar las dependencias estadísticas es el objetivo mismo.

1.2. Importancia de la estructura correcta#

Supongamos que la estructura verdadera del modelo es:

\[X \rightarrow Y \rightarrow Z\]

¿Qué pasa si falta un arco?

  • El modelo codifica independencias incorrectas.

  • La distribución verdadera \(P^*\) no puede ser representada.

  • Puede generalizar mejor (menos parámetros, menos sobreajuste).

¿Qué pasa si hay un arco de más?

  • El modelo codifica dependencias espurias.

  • Sí puede representar \(P^*\) correctamente.

  • Más parámetros → mayor riesgo de sobreajuste y peor generalización.

1.3. Dos enfoques principales#

Existen dos familias de métodos para aprender la estructura:

Enfoque

Idea central

Ventajas

Desventajas

Basado en puntuación

Definir una función que evalúa cada DAG y buscar el máximo

Flexible, robusto con ruido

Espacio de búsqueda exponencial

Basado en restricciones

Hacer pruebas de independencia estadística e inferir la estructura

Interpretable, usa conocimiento previo

Requiere muchos datos para tener poder estadístico

En esta sesión nos concentramos principalmente en el enfoque basado en puntuación.

1.4. Enfoque basado en puntuación: forma general#

Toda puntuación que regula complejidad sigue la forma:

\[\mathrm{score}(\mathcal{G}: \mathcal{D}) = \underbrace{\ell(\mathcal{G}: \mathcal{D})}_{\text{ajuste a los datos}} - \underbrace{\varphi(|\mathcal{D}|) \cdot \|\mathcal{G}\|}_{\text{penalización por complejidad}}\]

donde:

  • \(\ell(\mathcal{G}: \mathcal{D})\) es la log-verosimilitud.

  • \(\|\mathcal{G}\|\) es el número de parámetros del modelo.

  • \(\varphi(|\mathcal{D}|)\) controla cuánto penalizamos la complejidad.

La elección de \(\varphi\) define los distintos criterios:

Criterio

\(\varphi(D)\)

Notas

Sin penalización (MLE)

\(0\)

Sobreajuste, no se usa

AIC

\(1\)

Penalización fija

BIC

\(\frac{\log M}{2}\)

Penalización crece con los datos

2. Puntuación de máxima verosimilitud#

2.1. Definición#

La idea es encontrar la estructura \(\mathcal{G}\) que maximiza la log-verosimilitud de los datos:

\[\mathrm{score}_L(\mathcal{G}: \mathcal{D}) = \log \mathcal{L}\big((\hat{\theta}, \mathcal{G}): \mathcal{D}\big)\]

donde \(\hat{\theta} = \theta_{\text{MLE}}\) son los parámetros de máxima verosimilitud dado \(\mathcal{G}\) y \(\mathcal{D}\).

2.2. Ejemplo: dos variables \(X\) e \(Y\)#

Consideremos dos estructuras posibles:

Estructura \(\mathcal{G}_0\): \(X \quad Y\) (sin arco)

\[\mathrm{score}_L(\mathcal{G}_0: \mathcal{D}) = \sum_{d=1}^{M} \left(\log\hat{\theta}_{x[d]} + \log\hat{\theta}_{y[d]}\right)\]

Estructura \(\mathcal{G}_1\): \(X \rightarrow Y\) (con arco)

\[\mathrm{score}_L(\mathcal{G}_1: \mathcal{D}) = \sum_{d=1}^{M} \left(\log\hat{\theta}_{x[d]} + \log\hat{\theta}_{y[d]|x[d]}\right)\]

2.3. Diferencia entre estructuras e información mutua#

La diferencia entre ambas puntuaciones resulta ser:

\[\mathrm{score}_L(\mathcal{G}_1: \mathcal{D}) - \mathrm{score}_L(\mathcal{G}_0: \mathcal{D}) = M \cdot I_{\hat{P}}(X; Y)\]

donde \(I_{\hat{P}}(X; Y)\) es la información mutua empírica:

\[I_{\hat{P}}(X; Y) = \sum_{x,y} \hat{P}(x, y) \log \frac{\hat{P}(x,y)}{\hat{P}(x)\hat{P}(y)}\]

La información mutua mide qué tan lejos están \(X\) e \(Y\) de ser independientes.

2.4. Descomposición general#

Para cualquier red bayesiana con \(n\) variables:

\[\mathrm{score}_L(\mathcal{G}: \mathcal{D}) = M \sum_{i=1}^{n} I_{\hat{P}}(X_i, \mathrm{Pa}_{X_i}) - M \sum_{i=1}^{n} H_{\hat{P}}(X_i)\]

donde:

  • \(I_{\hat{P}}(X_i; \mathrm{Pa}_{X_i})\) es la información mutua entre el nodo y sus padres.

  • \(H_{\hat{P}}(X_i) = -\sum_{x} \hat{P}(x) \log \hat{P}(x)\) es la entropía del nodo.

La entropía no depende de la estructura \(\mathcal{G}\). Por tanto, maximizar la puntuación equivale a maximizar la información mutua entre cada nodo y sus padres.

2.5. Limitación: sobreajuste#

Dado que \(I_{\hat{P}}(X; Y) \geq 0\) siempre, y es cero solo si \(X\) e \(Y\) son independientes en \(\hat{P}\):

  • Agregar un arco nunca reduce la puntuación y casi siempre la aumenta.

  • Incluso si \(X \perp Y\) en \(P^*\), en datos finitos casi siempre \(I_{\hat{P}}(X;Y) > 0\).

La puntuación de máxima verosimilitud se maximiza con redes completamente conectadas → sobreajuste.

Por esta razón, esta puntuación no se usa en la práctica.

import numpy as np
import pandas as pd

M = 10000
np.random.seed(42)

# Datos
raw_data = np.random.randint(low=0, high=2, size=(M, 2))
data = pd.DataFrame(data=raw_data, columns=['X', 'Y'])
data['Z'] = (data['Y'] + np.random.randint(0, 2, M)) % 2

# Probabilidades empíricas
P_X  = data['X'].value_counts(normalize=True)
P_Y  = data['Y'].value_counts(normalize=True)
P_Z  = data['Z'].value_counts(normalize=True)
P_XY = data.groupby(['X','Y']).size() / M
P_YZ = data.groupby(['Y','Z']).size() / M
P_XZ = data.groupby(['X','Z']).size() / M

# Productos marginales 
P_XP_Y = pd.Series({(x,y): P_X[x] * P_Y[y] for (x,y) in P_XY.index})
P_YP_Z = pd.Series({(y,z): P_Y[y] * P_Z[z] for (y,z) in P_YZ.index})
P_XP_Z = pd.Series({(x,z): P_X[x] * P_Z[z] for (x,z) in P_XZ.index})

# Scores
score_G0 = M * ((P_X * np.log(P_X)).sum() +
                (P_Y * np.log(P_Y)).sum() +
                (P_Z * np.log(P_Z)).sum())

score_G1 = score_G0 + M * (P_XY * np.log(P_XY / P_XP_Y)).sum()
score_G2 = score_G1 + M * (P_YZ * np.log(P_YZ / P_YP_Z)).sum()
score_G3 = score_G2 + M * (P_XZ * np.log(P_XZ / P_XP_Z)).sum()

print(f"score_L(G0) = {score_G0:.4f}  [X    Y    Z]")
print(f"score_L(G1) = {score_G1:.4f}  [X -> Y    Z]")
print(f"score_L(G2) = {score_G2:.4f}  [X -> Y -> Z]")
print(f"score_L(G3) = {score_G3:.4f}  [X -> Y -> Z, X -> Z]")

score_L(G0) = -20794.3192  [X    Y    Z]
score_L(G1) = -20793.8192  [X -> Y    Z]
score_L(G2) = -20793.5459  [X -> Y -> Z]
score_L(G3) = -20793.1041  [X -> Y -> Z, X -> Z]

3. Puntuaciones con penalización por complejidad#

3.1. AIC: Akaike Information Criterion#

\[\mathrm{score}_{AIC}(\mathcal{G}: \mathcal{D}) = \ell(\mathcal{G}: \mathcal{D}) - \|\mathcal{G}\|\]

La penalización es fija: cada parámetro extra cuesta exactamente 1 en log-verosimilitud.

  • No es consistente: incluso con infinitos datos puede seleccionar un modelo más grande que el verdadero.

  • Es útil cuando el objetivo es predicción (no recuperar la estructura exacta).

BIC penaliza el arco extra: selecciona la estructura más simple y correcta.

def aic_score(log_likelihood, n_params):
    return log_likelihood - n_params

# Número de parámetros libres
# G0: P(X=1) + P(Y=1) = 2 parámetros
# G1: P(X=1) + P(Y=1|X=0) + P(Y=1|X=1) = 3 parámetros
params_G0 = 2
params_G1 = 3

aic_G0 = aic_score(score_G0, params_G0)
aic_G1 = aic_score(score_G1, params_G1)

print("── AIC ──")
print(f"score_L(G0) = {score_G0:.4f}   params = {params_G0}")
print(f"score_L(G1) = {score_G1:.4f}   params = {params_G1}")
print()
print(f"AIC(G0) = {score_G0:.4f} - {params_G0} = {aic_G0:.4f}")
print(f"AIC(G1) = {score_G1:.4f} - {params_G1} = {aic_G1:.4f}")

── AIC ──
score_L(G0) = -20794.3192   params = 2
score_L(G1) = -20793.8192   params = 3

AIC(G0) = -20794.3192 - 2 = -20796.3192
AIC(G1) = -20793.8192 - 3 = -20796.8192

3.2. BIC: Bayesian Information Criterion#

\[\mathrm{score}_{BIC}(\mathcal{G}: \mathcal{D}) = \ell(\mathcal{G}: \mathcal{D}) - \frac{\log M}{2} \cdot \|\mathcal{G}\|\]

La penalización crece con el tamaño del dataset: a más datos, mayor castigo por complejidad.

  • Asintóticamente consistente: conforme \(M \to \infty\), la información mutua (que crece con \(M\)) domina sobre el log del número de parámetros (que crece con \(\log M\)), recuperando \(\mathcal{G}^*\).

  • Es el criterio más usado en la práctica.

def bic_score(log_likelihood, n_params, M):
    return log_likelihood - (np.log(M) / 2) * n_params

bic_G0 = bic_score(score_G0, params_G0, M)
bic_G1 = bic_score(score_G1, params_G1, M)

phi_bic = np.log(M) / 2

print("── BIC ──")
print(f"φ(M) = log({M})/2 = {phi_bic:.4f}")
print()
print(f"BIC(G0) = {score_G0:.4f} - {phi_bic:.4f} × {params_G0} = {bic_G0:.4f}")
print(f"BIC(G1) = {score_G1:.4f} - {phi_bic:.4f} × {params_G1} = {bic_G1:.4f}")

── BIC ──
φ(M) = log(10000)/2 = 4.6052

BIC(G0) = -20794.3192 - 4.6052 × 2 = -20803.5296
BIC(G1) = -20793.8192 - 4.6052 × 3 = -20807.6347

print("── Comparativa AIC y BIC ──")
print(f"{'':20} {'G0':>10} {'G1':>10} {'Gana':>8}")
print(f"{'score_L':20} {score_G0:>10.4f} {score_G1:>10.4f}")
print(f"{'AIC':20} {aic_G0:>10.4f} {aic_G1:>10.4f} {'G0' if aic_G0 > aic_G1 else 'G1':>8}")
print(f"{'BIC':20} {bic_G0:>10.4f} {bic_G1:>10.4f} {'G0' if bic_G0 > bic_G1 else 'G1':>8}")

── Comparativa AIC y BIC ──
                             G0         G1     Gana
score_L              -20794.3192 -20793.8192
AIC                  -20796.3192 -20796.8192       G0
BIC                  -20803.5296 -20807.6347       G0

3.3. BD Score: Bayesian Dirichlet#

La puntuación Bayesiana toma un enfoque diferente: en lugar de penalizar explícitamente la complejidad, integra sobre todos los posibles valores de los parámetros usando una prior Dirichlet:

\[\mathrm{score}_{BD}(\mathcal{G}: \mathcal{D}) = \log P(\mathcal{D} \mid \mathcal{G}) = \log \int P(\mathcal{D} \mid \theta, \mathcal{G}) \, P(\theta \mid \mathcal{G}) \, d\theta\]

Con una prior Dirichlet de hiperparámetros \(\alpha_{x,u}\), esta integral tiene forma cerrada:

\[\mathrm{score}_{BD}(\mathcal{G}: \mathcal{D}) = \sum_{i=1}^{n} \sum_{u \in \mathrm{Pa}_{X_i}} \left[ \log \frac{\Gamma(\alpha_{u})}{\Gamma(\alpha_{u} + N_u)} + \sum_{x} \log \frac{\Gamma(\alpha_{x,u} + N_{x,u})}{\Gamma(\alpha_{x,u})} \right]\]
from pgmpy.estimators import K2
from pgmpy.models import DiscreteBayesianNetwork

# Modelos 
model_G0 = DiscreteBayesianNetwork()
model_G0.add_nodes_from(['X', 'Y'])

model_G1 = DiscreteBayesianNetwork([('X', 'Y')])

# K2 score (BD)
k2 = K2(data)

k2_G0 = k2.score(model_G0)
k2_G1 = k2.score(model_G1)

print("── K2 (BD) ───")
print(f"K2(G0) = {k2_G0:.4f}")
print(f"K2(G1) = {k2_G1:.4f}")
print()
if k2_G0 > k2_G1:
    print("G0 gana")
else:
    print("G1 gana")

── K2 (BD) ───
K2(G0) = -13871.6225
K2(G1) = -13874.8091

G0 gana

/home/docs/checkouts/readthedocs.org/user_builds/modelos-graficos-probabilisticos/envs/v2026/lib/python3.13/site-packages/pgmpy/estimators/__init__.py:4: FutureWarning: `pgmpy.estimators.StructureScore` is deprecated and will be removed in v1.3.0. Use `pgmpy.structure_score` instead.
  from .StructureScore import (

4. Ejemplo práctico: datos de vinos#

Hill-Climb Search: ¿cómo funciona?#

Algoritmo greedy que construye la red paso a paso.

Pasos#

  1. Empieza con una red vacía: sin ningún arco.

  2. Prueba todas las operaciones posibles:

    • Añadir un arco

    • Quitar un arco

    • Invertir un arco

  3. Elige la operación que más suba el score.

  4. Repite hasta que ningún cambio mejore el score.

Ventaja#

Eficiente y fácil de combinar con cualquier score (AIC, BIC, K2).

Limitación#

Es greedy: solo ve el siguiente paso, no el panorama completo. Puede quedarse en un óptimo local, no necesariamente el global.

¿Por qué lo usamos aquí?#

Con 4 variables es suficiente para explorar cómo distintos scores producen estructuras distintas.

import os
import pandas as pd
import numpy as np
from pgmpy.models import DiscreteBayesianNetwork
from pgmpy.estimators import HillClimbSearch
import warnings
warnings.filterwarnings('ignore')

ruta = os.path.join('..', 'data', 'winequality-red.csv')
df = pd.read_csv(ruta, sep=';')
print(df.shape)
df.head()

(1599, 12)
fixed acidity volatile acidity citric acid residual sugar chlorides free sulfur dioxide total sulfur dioxide density pH sulphates alcohol quality
0 7.4 0.70 0.00 1.9 0.076 11.0 34.0 0.9978 3.51 0.56 9.4 5
1 7.8 0.88 0.00 2.6 0.098 25.0 67.0 0.9968 3.20 0.68 9.8 5
2 7.8 0.76 0.04 2.3 0.092 15.0 54.0 0.9970 3.26 0.65 9.8 5
3 11.2 0.28 0.56 1.9 0.075 17.0 60.0 0.9980 3.16 0.58 9.8 6
4 7.4 0.70 0.00 1.9 0.076 11.0 34.0 0.9978 3.51 0.56 9.4 5 
# alcohol, quality, sulphates, volatile acidity
vars_sel = ['alcohol', 'quality', 'sulphates', 'volatile acidity']
df_disc = df[vars_sel].copy()

# Discretizar por mediana: 0 = bajo, 1 = alto
for col in vars_sel:
    mediana = df_disc[col].median()
    df_disc[col] = (df_disc[col] > mediana).astype(int)

import logging
logging.getLogger('pgmpy').setLevel(logging.WARNING)

# Hill-Climb con cada score
for score in ['bic-d', 'aic-d', 'k2']:
    est = HillClimbSearch(df_disc)
    modelo = est.estimate(scoring_method=score, max_indegree=3)
    print(f"── {score} ───")
    print(f"Arcos: {list(modelo.edges())}\n")

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[12], line 7
      3 
      4 # Hill-Climb con cada score
      5 for score in ['bic-d', 'aic-d', 'k2']:
      6     est = HillClimbSearch(df_disc)
----> 7     modelo = est.estimate(scoring_method=score, max_indegree=3)
      8     print(f"── {score} ───")
      9     print(f"Arcos: {list(modelo.edges())}\n")

File ~/checkouts/readthedocs.org/user_builds/modelos-graficos-probabilisticos/envs/v2026/lib/python3.13/site-packages/pgmpy/estimators/HillClimbSearch.py:212, in HillClimbSearch.estimate(self, scoring_method, start_dag, tabu_length, max_indegree, expert_knowledge, epsilon, max_iter, show_progress)
    146 """
    147 Performs local hill climb search to estimates the `DAG` structure that
    148 has optimal score, according to the scoring method supplied. Starts at
   (...)    207 
    208 """
    210 # Step 1: Initial checks and setup for arguments
    211 # Step 1.1: Check scoring_method
--> 212 score, score_c = get_scoring_method(scoring_method, self.data, self.use_cache)
    213 score_fn = score_c.local_score
    215 # Step 1.2: Check the start_dag

File ~/checkouts/readthedocs.org/user_builds/modelos-graficos-probabilisticos/envs/v2026/lib/python3.13/site-packages/pgmpy/estimators/StructureScore.py:68, in get_scoring_method(scoring_method, data, use_cache, **kwargs)
     64         raise ValueError(
     65             "Unknown scoring method. Please refer documentation for a list of supported score metrics."
     66         )
     67     elif scoring_method.lower() not in list(supported_methods.keys()):
---> 68         raise ValueError(
     69             f"Incorrect scoring method for {var_type}, scoring_method should be one of"
     70             f"{list(supported_methods.keys())}, received {scoring_method}. {data.dtypes.unique()}"
     71         )
     72 elif isinstance(scoring_method, type(None)):
     73     # automatically determine scoring method, pick first one
     74     scoring_method = list(available_methods[var_type].keys())[0]

ValueError: Incorrect scoring method for continuous, scoring_method should be one of['bic-g', 'll-g', 'aic-g', 'bic-cg', 'll-cg', 'aic-cg'], received bic-d. [dtype('int64')]

BIC-D#

Arcos encontrados: 5

Origen

Destino

quality

alcohol

quality

sulphates

sulphates

alcohol

sulphates

volatile acidity

alcohol

volatile acidity

Estructura:

quality → alcohol → volatile acidity
quality → sulphates → volatile acidity
          sulphates → alcohol

from IPython.display import display, HTML

display(HTML("""
<script src="https://cdn.jsdelivr.net/npm/mermaid/dist/mermaid.min.js"></script>
<script>mermaid.initialize({startOnLoad:true, theme:'default'});</script>
<h3>BIC-D — 5 arcos</h3>
<div class="mermaid">
graph TD
    Q[quality] --> A[alcohol]
    Q --> S[sulphates]
    S --> A
    S --> V[volatile acidity]
    A --> V
</div>
<script>setTimeout(function() { mermaid.init(); }, 500);</script>
"""))

BIC-D — 5 arcos

graph TD Q[quality] --> A[alcohol] Q --> S[sulphates] S --> A S --> V[volatile acidity] A --> V

AIC-D#

Arcos encontrados: 6

Origen

Destino

quality

alcohol

quality

sulphates

quality

volatile acidity

sulphates

alcohol

sulphates

volatile acidity

alcohol

volatile acidity

Estructura:

quality → alcohol → volatile acidity
quality → sulphates → volatile acidity
quality → volatile acidity (directo)
          sulphates → alcohol

display(HTML("""
<h3>AIC-D — 6 arcos</h3>
<div class="mermaid">
graph TD
    Q[quality] --> A[alcohol]
    Q --> S[sulphates]
    Q --> V[volatile acidity]
    S --> A
    S --> V
    A --> V
</div>
<script>setTimeout(function() { mermaid.init(); }, 500);</script>
"""))

AIC-D — 6 arcos

graph TD Q[quality] --> A[alcohol] Q --> S[sulphates] Q --> V[volatile acidity] S --> A S --> V A --> V

K2#

Arcos encontrados: 6

Origen

Destino

quality

alcohol

quality

sulphates

quality

volatile acidity

alcohol

sulphates

alcohol

volatile acidity

sulphates

volatile acidity

Estructura:

quality → alcohol → sulphates → volatile acidity
quality → sulphates → volatile acidity
quality → volatile acidity (directo)
alcohol → volatile acidity (directo)

display(HTML("""
<h3>K2 — 6 arcos</h3>
<div class="mermaid">
graph TD
    Q[quality] --> A[alcohol]
    Q --> S[sulphates]
    Q --> V[volatile acidity]
    A --> S
    A --> V
    S --> V
</div>
<script>setTimeout(function() { mermaid.init(); }, 500);</script>
"""))

K2 — 6 arcos

graph TD Q[quality] --> A[alcohol] Q --> S[sulphates] Q --> V[volatile acidity] A --> S A --> V S --> V

Lectura recomendada#

https://pgmpy.org/examples/Structure_Learning.html

Koller & Friedman: Sección 18.3.2: Bayesian Score (página 794)

Introduce la puntuación Bayesiana completa con prior Dirichlet. Muestra cómo la integración sobre parámetros actúa como regularización natural, con mejores propiedades que BIC en muestras pequeñas.

CS228 Lecture Notes: Structure Learning (Ermon Group, Stanford) https://ermongroup.github.io/cs228-notes/learning/structure/

Complemento con comparación AIC/BIC/BD y métodos de búsqueda de estructura.