Sesión 10 extra 2

Contenido

Sesión 10 extra 2#

Inferencia aproximada#

Objetivos:

Estudiar inferencia aproximada y métodos de muestreo en modelos gráficos probabilísticos.

Recursos de consulta:

Probabilistic Graphical Models: Principles and Techniques, By Daphne Koller & Nir Friedman. Ch. 12.

Documentación pgmpy:

Pgmpy documentation

1.Introducción#

En inferencia exacta (por ejemplo, mediante Variable Elimination), el objetivo es calcular distribuciones marginales de una red bayesiana de forma analítica.

Sin embargo, este enfoque se vuelve computacionalmente costoso cuando el número de variables y conexiones crece.

Para resolver esto, la inferencia aproximada busca aproximar las distribuciones de probabilidad mediante métodos de muestreo (sampling).

1.1. Muestreo basado en partículas (Particle-based approximate inference)#

En lugar de representar la distribución de una red bayesiana mediante factores (como se hace en inferencia exacta), podemos aproximar la distribución conjunta mediante un conjunto finito de instanciaciones de las variables, llamadas partículas.

“In these methods, we approximate the joint distribution as a set of instantiations to all or some of the variables in the network.

These instantiations, often called particles, are designed to provide a good representation of the overall probability distribution.”

— Koller & Friedman (2009, Ch. 12)

En este enfoque, se representa la distribución conjunta $P(\mathcal{X})$ como un conjunto finito de instanciaciones aleatorias llamadas partículas:

\[ \mathcal{D} = \{ \xi[1], \xi[2], \dots, \xi[M] \} \]

Cada partícula $\xi[m]$ es una asignación de valores a las variables del modelo.

La idea central es que, si las partículas son generadas de forma adecuada, el conjunto $\mathcal{D}$ aproxima bien la distribución $P(\mathcal{X})$.

Muestreo por partículas#

Para generar partículas $\xi[m]$, se sigue el orden causal de la red (en este caso: $I, D, C, E, R$).

Cada partícula corresponde a una posible realidad simulada del modelo.

Ejemplo de muestreo:

Partícula	D (Dificultad)	I (Inteligencia)	C (Calificación)	E (Prueba)	R (Carta)
$\xi[1]$	Fácil ($d^0$)	Inteligente ($i^1$)	Alta ($c^2$)	Bueno ($e^1$)	Buena ($r^1$)
$\xi[2]$	Difícil ($d^1$)	No mucho ($i^0$)	Media ($c^1$)	Malo ($e^0$)	Mala ($r^0$)
$\xi[3]$	Fácil ($d^0$)	No mucho ($i^0$)	Baja ($c^0$)	Malo ($e^0$)	Mala ($r^0$)
…	…	…	…	…	…
$\xi[M]$	Difícil ($d^1$)	Inteligente ($i^1$)	Alta ($c^2$)	Bueno ($e^1$)	Buena ($r^1$)

Cada partícula $\xi[m]$ representa una instanciación completa del conjunto de variables, generada de acuerdo con las distribuciones condicionales del modelo.

1.2. Forward sampling en redes bayesianas#

Perfecto, y ¿cómo obtengo esas muestras de una red bayesiana?

Acá entra el forward sampling, método de simulación que nos permite generar partículas $xi[m]$ coherentes con la estructura causal del modelo.

El método de forward sampling consiste en generar partículas siguiendo el orden causal de la red bayesiana.

Algoritmo – Forward Sampling en una red bayesiana

Forward Sampling

Figura 1. Algoritmo de forward sampling para generar partículas en una red bayesiana. Koller & Friedman (2009, Ch. 12)

Paso	Descripción
(1) Orden topológico	Se ordenan las variables $X_1, \dots, X_n$ de modo que los padres aparezcan antes que sus hijos. Esto garantiza que cuando llegamos a una variable $X_i$, ya conocemos los valores de sus padres.
(2–3) Obtener valores de los padres	Para cada variable $X_i$, se obtienen los valores actuales $u_i$ de sus padres $\text{Pa}(X_i)$, que ya fueron muestreados.
(4) Muestreo condicional	Se genera un valor $x_i$ de acuerdo con la distribución condicional $P(X_i \mid u_i)$. En la práctica, se genera un número aleatorio $r \sim \text{Uniform}(0,1)$ y se selecciona el valor de $X_i$ correspondiente al intervalo acumulado que contiene $r$.
(5) Retorno de la partícula	El algoritmo devuelve la asignación completa $(x_1, \dots, x_n)$, que constituye una partícula del modelo.

💡 Ejemplo: Forward Sampling en la red del estudiante

Supongamos que queremos generar una partícula $\xi = (x_I,x_D,x_C,x_E,x_R)$ siguiendo el orden topológico:

\[ I \rightarrow D \rightarrow C \rightarrow E \rightarrow R \]

Paso 1️ – Muestreo de la inteligencia $I$

Padres: $\mathrm{Pa}(I)=\varnothing \Rightarrow u_I = \varnothing$

Distribución: $$ P(I) = \begin{cases} P(i^0) = 0.7 & \text{(no muy inteligente)} \\ P(i^1) = 0.3 & \text{(inteligente)} \end{cases} $$

Generamos un número aleatorio $r = 0.65$

Acumuladas: $[0,0.7)\Rightarrow i^0;\; [0.7,1]\Rightarrow i^1$

Como $r=0.65\in[0,0.7)$, $x_I=i^0$

Paso 2️ – Muestreo de la dificultad $D$

Padres: $\mathrm{Pa}(D)=\varnothing \Rightarrow u_D = \varnothing$

Distribución: $$ P(D) = \begin{cases} P(d^0) = 0.6 & \text{(fácil)} \\ P(d^1) = 0.4 & \text{(difícil)} \end{cases} $$

Generamos $r = 0.82$

Acumuladas: $[0,0.6)\Rightarrow d^0;\; [0.6,1]\Rightarrow d^1$

Como $r=0.82\in[0.6,1]$, $x_D=d^1$

Paso 3️ – Muestreo de la calificación $C \sim P(C \mid I, D)$

Padres: $\mathrm{Pa}(C)=\{I,D\}\Rightarrow u_C=(i^0,d^1)$

Dado que $I = i^0$ y $D = d^1$, usamos la tabla condicional:

Calificación $C$	$c^0$ (baja)	$c^1$ (media)	$c^2$ (alta)
Probabilidad	0.7	0.25	0.05

Generamos $r=0.31$

Acumuladas:
$[0,0.70)\Rightarrow c^0;\; [0.70,0.95)\Rightarrow c^1;\; [0.95,1]\Rightarrow c^2$.

Como $r=0.31\in[0,0.70)$, $x_C=c^0$

Paso 4️ – Muestreo de la prueba $E \sim P(E \mid I)$

Padres: $\mathrm{Pa}(E)=\{I\}\Rightarrow u_E=(i^0)$

Dado $I = i^0$, usamos:

Resultado $E$	$e^0$ (malo)	$e^1$ (bueno)
Probabilidad	0.95	0.05

Generamos $r = 0.12$

Acumuladas: $[0,0.95)\Rightarrow e^0;\; [0.95,1]\Rightarrow e^1$

Como $r=0.12\in[0,0.95)$, $x_E=e^0$

Paso 5️ – Muestreo de la carta $R \sim P(R \mid C)$

Padres: $\mathrm{Pa}(R)=\{C\}\Rightarrow u_R=(c^0)$

Dado $C = c^0$:

Carta $R$	$r^0$ (mala)	$r^1$ (buena)
Probabilidad	0.99	0.01

Generamos $r = 0.27$

Acumuladas: $[0,0.99)\Rightarrow r^0;\; [0.99,1]\Rightarrow r^1$

Como $r=0.27\in[0,0.99)$, $x_R=r^0$

Resultado final

Hemos generado la siguiente partícula:

\[ \xi = (x_I,x_D,x_C,x_E,x_R) = (i^0,\; d^1,\; c^0,\; e^0,\; r^0) \]

Cada vez que repitamos el proceso con diferentes números aleatorios, obtendremos nuevas partículas $\xi[m]$.

El conjunto de todas las partículas forma el conjunto de muestras:

\[ \mathcal{D} = \{ \xi[1], \xi[2], \dots, \xi[M] \} \]

El conjunto $\mathcal{D}$ aproxima la distribución conjunta $P(I, D, C, E, R)$ del modelo bayesiano.

Acumulación de probabilidades condicionales

Esta es una forma de entender cómo se acumulan las probabilidades condicionales para muestrear una variable $X_i$ dado los valores de sus padres $u_i$:

\[ \underbrace{0}_{\text{inicio}} \;\; \overbrace{P(X_i = x_i^1 \mid u_i)}^{\text{primer intervalo}} \;\; \overbrace{+\, P(X_i = x_i^2 \mid u_i)}^{\text{siguiente}} \;\; \overbrace{+\, P(X_i = x_i^3 \mid u_i)}^{\text{siguiente}} \;\; 1 \]

1.3. Estimación mediante Monte Carlo#

El promedio de las muestras (la frecuencia relativa) es una buena estimación de la probabilidad o expectativa verdadera en la red bayesiana.

El objetivo de forward sampling no es generar una sola partícula, sino un conjunto de partículas:

\[ \mathcal{D} = \{ \xi[1], \xi[2], \dots, \xi[M] \}, \]

que constituye una muestra del modelo bayesiano.

A partir de este conjunto $\mathcal{D}$, podemos estimar probabilidades o expectativas, por ejemplo:

\[ \boxed{P(Y = y)} \]

$P(R = r^1)$: la probabilidad de recibir una carta, o

\[ \boxed{P(Y =y \mid E=e)} \]

una probabilidad condicional con evidencia, por ejemplo $P(R = r^1 \mid E = e^0)$: la probabilidad de carta dado que el puntaje del examen fue malo.

Estas estimaciones se obtienen evaluando cuántas partículas cumplen la condición respecto al total generado, o más generalmente, promediando funciones sobre las partículas.

¿Cómo?

Para evaluar estas probabilidades, recuperamos la definición de la esperanza de una función $f(X)$ bajo la distribución $P(X)$:

\[ \begin{aligned} E_P[f(X)] = \sum_X f(X) P(X) \end{aligned} \tag{1} \]

Ecuación (1). Definición de la esperanza matemática: el promedio ponderado de los valores de $f(X)$ bajo la distribución $P(X)$.

y no podemos hacerlo de forma exacta, es decir, no podemos evaluar directamente $P(x)$, entonces aproximamos con las muestras:

\[ \begin{aligned} \hat{E}_\mathcal{D}[f] = \frac{1}{M} \sum_{m=1}^{M} f(\xi[m]) \end{aligned} \tag{2} \]

Ecuación (2). Aproximación mediante Monte Carlo. En lugar de recorrer todos los posibles $X$, tomamos $M$ muestras aleatorias $\xi[m]$ del modelo y calculamos el promedio de $f$ sobre ellas.

Así, reemplzamos la esperanza exacta (suma sobre todos los $X$) por una media muestral (suma sobre las partículas generadas).

Podemos expresar la probabilidad como una expectativa particular: la esperanza de una función indicadora que vale 1 cuando el evento ocurre y 0 en otro caso.

De esta forma, la fórmula general de Monte Carlo se convierte directamente en una estimación de probabilidad:

\[ \begin{aligned} \hat{P}_\mathcal{D}(y) = \frac{1}{M} \sum_{m=1}^{M} \mathbf{I}\{y[m] = y\} \end{aligned} \tag{3} \]

Ecuación (3). Finalmente, elegimos una función específica: $F(X) = \mathbb{1}_{Y = y}$, que nos permite estimar la probabilidad de que $Y = y$ contando cuántas partículas cumplen esa condición.

Al sustituir Ecuación (3) en Ecuación (2), la expectativa estimada se convierte en la frecuencia relativa del evento $Y=y$ en las partículas generadas; es decir, cuántas veces ocurre $Y=y$ dividido entre el total de muestras $M$.

Nota:
La función indicadora $\mathbf{I}\{Y = y\}$ es una función que vale 1 cuando el evento $Y = y$ ocurre y 0 en caso contrario:

\[\begin{split} > \mathbf{I}\{Y = y\} = > \begin{cases} > 1, & \text{si } Y = y \\ > 0, & \text{en otro caso} > \end{cases} > \end{split}\]

Al tomar la esperanza de esta función respecto a la distribución $P(Y)$:

\[ > \mathbb{E}_P[\mathbf{I}\{Y = y\}] = \sum_Y \mathbf{I}\{Y = y\} P(Y) > \]

todos los términos son cero excepto cuando $Y = y$, por lo que obtenemos directamente:

\[ > \mathbb{E}_P[\mathbf{I}\{Y = y\}] = P(Y = y) > \]

Es decir, la probabilidad de un evento es la esperanza de su función indicadora.

Por eso, al usar Monte Carlo y promediar los valores de $\mathbf{I}\{Y[m] = y\}$ en las muestras, estamos estimando $P(Y=y)$ mediante la frecuencia relativa de partículas en las que $Y = y$.

💡 Ejemplo: estimación de una probabilidad en la red del estudiante

Supongamos que queremos estimar la probabilidad de que la variable Carta ($R$) viable, es decir, $P(R = r^1)$, usando las partículas generadas por forward sampling.

Si obtenemos $M = 5$ muestras del modelo:

Muestra ($m$)	$R[m]$	$\mathbf{I}\{R[m] = r^1\}$
1	$r^0$	0
2	$r^1$	1
3	$r^1$	1
4	$r^0$	0
5	$r^1$	1

El promedio de estos valores nos da la estimación:

\[ \hat{P}_\mathcal{D}(R = r^1) = \frac{1}{5}(0 + 1 + 1 + 0 + 1) = 0.6 \]

Por lo tanto, con estas muestras estimamos que:

\[ P(R = r^1) \approx 0.6 \]

Interpretación:
El resultado $0.6$ significa que la carta es buena en el 60 % de las partículas generadas.

💡 Ejemplo: estimación de una probabilidad condicional con evidencia

Supongamos que ahora queremos estimar la probabilidad de que la variable Carta ($R$) sea buena, dado que la calificación fue alta, es decir:

\[ P(R = r^1 \mid C = c^2) \]

Generamos $M = 8$ partículas con forward sampling y luego filtramos las que cumplen la evidencia $C=c^2$:

Muestra ($m$)	$C[m]$	$R[m]$	$\mathbf{I}\{R[m] = r^1\}$
1	$c^2$	$r^1$	1
2	$c^0$	$r^0$	—
3	$c^2$	$r^1$	1
4	$c^1$	$r^0$	—
5	$c^2$	$r^1$	1
6	$c^0$	$r^0$	—
7	$c^2$	$r^1$	1
8	$c^2$	$r^0$	0

Solo consideramos las filas con $C=c^2$ (muestras 1, 3, 5, 7, 8).

En total, $N=5$ partículas cumplen la evidencia y 4 de ellas tienen $R=r^1$

La estimación Monte Carlo de la probabilidad condicional es: $$ \hat{P}_\mathcal{D}(R=r^1 \mid C=c^2) =\frac{\sum_{m=1}^M \mathbf{I}\{R[m]=r^1,\,C[m]=c^2\}} {\sum_{m=1}^M \mathbf{I}\{C[m]=c^2\}} =\frac{4}{5}=0.8 $$

1.4. Análisis del error#

Hasta aquí, hemos visto cómo estimar una probabilidad $\hat{P}_{\mathcal{D}}(Y=y)$ promediando las muestras generadas mediante forward sampling.

Cuando estimamos:

\[ \hat{P}(y)= \frac{\text{\# veces que ocurre} y}{M} \]

esta estimación varía cada vez que repetimos el muestreo.

¿Por qué? Porque las muestras son aleatorias, entonces hay dos preguntas clave:

¿qué tan lejos puede quedar mi estimación del valor real $P(y)$?

Cuántas muesstras $M$ necesito para estar segur@ de que el error será pequeño?

A continuación, analizamos cómo varía el error del estimador y cómo podemos acotarlo usando la desigualdad de Hoeffding.

1.4.1. Planteamiento matemático del error#

Queremos cuantificar la probabilidad de que el estimador difiera del valor verdadero más de un margen $\varepsilon$:

\[\begin{split} \begin{aligned} \underbrace{ P_\mathcal{D} \Big( \underbrace{ \big| \underbrace{\hat{P}_\mathcal{D}(y)}_{\text{estimación empírica}} - \underbrace{P(y)}_{\text{valor verdadero}} \big| }_{\text{error absoluto}} \ge \underbrace{\varepsilon}_{\substack{\text{tolerancia predefinida} \\ \text{(hiperparámetro de precisión)}}} \Big) }_{\text{probabilidad de que el error exceda el umbral } \varepsilon} \end{aligned} \end{split}\]

Esta expresión mide qué tan confiable es nuestra estimación.

¿Qué probabilidad hay de que mi estimación se equivoque más de lo que estoy dispuesta a tolerar?

Indica la probabilidad de que el error absoluto entre la estimación $\hat{P}_\mathcal{D}(y)$ y el valor verdadero $P(y)$ sea mayor o igual a una tolerancia predefinida $\varepsilon$.

1.4.2. Cota de Hoeffding#

Aplicando la desigualdad de Hoeffding, obtenemos:

\[ \underbrace{ P_\mathcal{D} \big( |\hat{P}_\mathcal{D}(y) - P(y)| \geq \varepsilon \big) }_{\text{probabilidad de que el error supere } \varepsilon} \;\; \leq \;\; \underbrace{ 2 e^{-2M \varepsilon^2} }_{\text{cota superior (Hoeffding)}} \tag{5} \]

La desigualdad de Hoeffding nos da una cota teórica sobre la probabilidad de cometer un error mayor que $\varepsilon$:

cuando más grande sea $M$, menor será el término $2 e^{-2M \varepsilon^2}$

Esto significa que la probabilidad de que el estimador se desvíe mucho del valor real, disminuye de forma exponencial con el número de partículas.

1.4.3. Determinación de la muestra#

En la práctica, podemos decidir:

Quiero que la probabilidad de equivocarme sea como mucho $\delta$ o equivalentemente, tener una confianza de $1 - \delta$ en mi estimación.

Entonces, imponemos que la cota de Hoeffding sea menor o igual a $\delta$:

\[ 2e^{-2M\varepsilon^2} \leq \delta \tag{6} \]

y despejamos el número mínimo de muestras necesarias:

\[ M \geq \frac{\ln(2/\delta)}{2\varepsilon^2} \tag{7} \]

Así, garantizamos que: con al menos $1 - \delta$ de confianza, el error será menor que $\varepsilon$.

1.5. Ejemplo en Código (forward sampling y análisis de error para $P(Y=y)$)#

from pgmpy.sampling import BayesianModelSampling
from pgmpy.inference import VariableElimination
import os
import numpy as np
import pickle

/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 (

Importamos el modelo bayesiano del estudiante.

ruta = os.path.join('..', 'data', 'student-model.pkl')

with open(ruta, "rb") as f:
    student_model = pickle.load(f)

print(type(student_model))

<class 'pgmpy.models.DiscreteBayesianNetwork.DiscreteBayesianNetwork'>

Instanciamos el muestreador BayesianNetworkSampler de partículas.

sampler = BayesianModelSampling(student_model)

Calculamos n_samples (Hoeffding)

\[ M \geq \frac{\ln(2/\delta)}{2\varepsilon^2} \quad \]

# Calcular el número de muestras para una precisión $\varepsilon$
# y un nivel de confianza $1 - \delta$

delta = 0.01
eps = 0.01

n_samples = np.log(2 / delta) / (2 * eps**2)
n_samples

np.float64(26491.58683274018)

# Tamaño de muestra Hoeffding
n_samples = np.ceil(n_samples).astype(int)
n_samples

np.int64(26492)

Generamos partículas mediante forward_sample

# Generamos $xi[m]$ muestras independientes del modelo
samples = sampler.forward_sample(size=n_samples)
samples.head()

	S	B	E	L	N
0	0	1	0	1	0
1	1	1	0	0	1
2	0	1	0	0	1
3	0	0	0	0	1
4	1	1	1	0	1

Monte Carlo: calcular frecuencias/expectativas.

# Estimación de $P(Y=y)$ para el caso de $Y = C$
p_C_approx = samples.groupby('C').size() / n_samples
p_C_approx

C
0    0.765401
1    0.234599
dtype: float64

1.6. Comparación con Inferencia exacta#

inference = VariableElimination(student_model)
p_C_exact = inference.query(variables=['C'])
print(p_C_exact)

+------+----------+
| C    |   phi(C) |
+======+==========+
| C(0) |   0.7650 |
+------+----------+
| C(1) |   0.2350 |
+------+----------+

Ejemplo con evidencia (Hoeffding)

\[ P(C = i^1 \mid I = i^0) \]

p_C_giv_I0 = (
    samples[samples['I'] == 0]
    .groupby('C')
    .size()
    / len(samples[samples['I'] == 0])
)
p_C_giv_I0

---------------------------------------------------------------------------
KeyError                                  Traceback (most recent call last)
File ~/checkouts/readthedocs.org/user_builds/modelos-graficos-probabilisticos/envs/v2026/lib/python3.13/site-packages/pandas/core/indexes/base.py:3641, in Index.get_loc(self, key)
   3640 try:
-> 3641     return self._engine.get_loc(casted_key)
   3642 except KeyError as err:

File pandas/_libs/index.pyx:168, in pandas._libs.index.IndexEngine.get_loc()
--> 168 'Could not get source, probably due dynamically evaluated source code.'

File pandas/_libs/index.pyx:197, in pandas._libs.index.IndexEngine.get_loc()
--> 197 'Could not get source, probably due dynamically evaluated source code.'

File pandas/_libs/hashtable_class_helper.pxi:7668, in pandas._libs.hashtable.PyObjectHashTable.get_item()
-> 7668 'Could not get source, probably due dynamically evaluated source code.'

File pandas/_libs/hashtable_class_helper.pxi:7676, in pandas._libs.hashtable.PyObjectHashTable.get_item()
-> 7676 'Could not get source, probably due dynamically evaluated source code.'

KeyError: 'I'

The above exception was the direct cause of the following exception:

KeyError                                  Traceback (most recent call last)
Cell In[9], line 2
      1 p_C_giv_I0 = (
----> 2     samples[samples['I'] == 0]
      3     .groupby('C')
      4     .size()
      5     / len(samples[samples['I'] == 0])

File ~/checkouts/readthedocs.org/user_builds/modelos-graficos-probabilisticos/envs/v2026/lib/python3.13/site-packages/pandas/core/frame.py:4378, in DataFrame.__getitem__(self, key)
   4374 
   4375         if is_single_key:
   4376             if self.columns.nlevels > 1:
   4377                 return self._getitem_multilevel(key)
-> 4378             indexer = self.columns.get_loc(key)
   4379             if is_integer(indexer):
   4380                 indexer = [indexer]
   4381         else:

File ~/checkouts/readthedocs.org/user_builds/modelos-graficos-probabilisticos/envs/v2026/lib/python3.13/site-packages/pandas/core/indexes/base.py:3648, in Index.get_loc(self, key)
   3643     if isinstance(casted_key, slice) or (
   3644         isinstance(casted_key, abc.Iterable)
   3645         and any(isinstance(x, slice) for x in casted_key)
   3646     ):
   3647         raise InvalidIndexError(key) from err
-> 3648     raise KeyError(key) from err
   3649 except TypeError:
   3650     # If we have a listlike key, _check_indexing_error will raise
   3651     #  InvalidIndexError. Otherwise we fall through and re-raise
   3652     #  the TypeError.
   3653     self._check_indexing_error(key)

KeyError: 'I'

# para contrastar con el valor exacto
p_C_given_I0 = inference.query(variables=['C'], evidence={'I': 0})
print(p_C_given_I0)

+------+----------+
| C    |   phi(C) |
+======+==========+
| C(0) |   0.4600 |
+------+----------+
| C(1) |   0.3400 |
+------+----------+
| C(2) |   0.2000 |
+------+----------+

1.7. Análisis con cota de Chernoff#

Ahora, ¿qué pasa si necesitamos estar seguros de cuántas observaciones hacen falta para que la estimación de una clase poco frecuente sea precisa?

Si $P(y)$ es muy pequeño (evento raro o clase minoritaria), entonces se necesitan muchos más ejemplos para que el error relativo sea pequeño.

De hecho, si $P(y)$ es tan pequeño que las $M$ muestras casi nunca ocurre, la estimación empírica puede ser $0$, lo que está infinitamente lejos del valor verdader en términos relativos.

La desigualdad de Chernoff refina el análisis anterior para errores relativos, obteniendo:

\[ P_D\left( |\hat{P}_\mathcal{D}(y) - P(y)| \ge \epsilon P(y) \right) \le 2 e^{-M P(y) \epsilon^2 / 3} \]

Despejando $M$:

\[ M \ge \frac{3 \ln(2 / \delta)}{P(y) \epsilon^2} \]

Esto implica que si el evento $Y=y$ es poco probable (por ejemplo, $P(y) \ll 1$), se requerirá un número mucho mayor de muestras para observarlo lo suficiente y estimarlo con precisión.

# Calcular el número de muestras para una precisión $\varepsilon$
# y un nivel de confianza $1 - \delta$
# además, para una variable Y específica

delta = 0.01
eps_relativo = 0.01
y_state = 2  # queremos estimar P(C = c2)

# Probabilidad "verdadera" de C
infer = VariableElimination(student_model)
dist_C = infer.query(variables=['C'])
p_y = dist_C.get_value(C=2)

p_y

np.float64(0.36200000000000004)

# Tamaño de muestra Chernoff
n_samples = np.ceil( 3 * np.log(2/delta) / (p_y * eps_relativo**2) ).astype(int)
n_samples

np.int64(439088)

# Generamos $xi[m]$ muestras independientes del modelo
samples = sampler.forward_sample(size=n_samples)
samples.head()

	D	C	I	E	R
0	0	2	0	0	1
1	0	2	1	1	1
2	0	1	0	0	1
3	0	1	0	0	1
4	1	0	0	0	0

# Estimación de P(C)
p_C_approx = samples.groupby('C').size() / n_samples

print(p_C_approx)

C
0    0.349857
1    0.288097
2    0.362046
dtype: float64

print(dist_C)

+------+----------+
| C    |   phi(C) |
+======+==========+
| C(0) |   0.3496 |
+------+----------+
| C(1) |   0.2884 |
+------+----------+
| C(2) |   0.3620 |
+------+----------+

Usar Chernoff para $C=c^2$ nos da una cota conservadora (segura) que también garantiza buena precisión para los demás valores de $C$.

Es decir, si tenemos suficientes muestras para estimar el caso más dificil, también estamos cubiertos para los demás casos.

Ejemplo con evidencia (Chernoff)

# Params
delta = 0.01         
eps_relativo = 0.01       
y_state = 2

evidencia_tup = [('I', 0)]

infer = VariableElimination(student_model)
dist_C_given_I0 = infer.query(variables=['C'],
                             evidence={'I': 0})

p_y = dist_C_given_I0.get_value(C=y_state)

n_samples = np.ceil( 3 * np.log(2/delta) / (p_y * eps_relativo**2) ).astype(int)
n_samples

np.int64(794748)

Generamos partículas mediante rejection_sample

samples = sampler.rejection_sample(evidence=evidencia_tup,
                                   size=n_samples)
samples

	D	C	I	E	R
0	0	0	0	0	0
1	0	0	0	0	0
2	1	2	0	0	1
3	0	0	0	0	0
4	1	1	0	0	0
...	...	...	...	...	...
794743	1	2	0	0	1
794744	0	1	0	0	0
794745	0	2	0	0	1
794746	1	0	0	0	0
794747	0	2	0	0	1

794748 rows × 5 columns

samples.I.value_counts()

I
0    794748
Name: count, dtype: int64

p_C_given_E_approx = samples.groupby('C').size() / samples.shape[0]
p_C_given_E_approx

C
0    0.459727
1    0.340429
2    0.199844
dtype: float64

print(dist_C_given_I0)

+------+----------+
| C    |   phi(C) |
+======+==========+
| C(0) |   0.4600 |
+------+----------+
| C(1) |   0.3400 |
+------+----------+
| C(2) |   0.2000 |
+------+----------+

Partícula	D (Dificultad)	I (Inteligencia)	C (Calificación)	E (Prueba)	R (Carta)
\(\xi[1]\)	Fácil (\(d^0\))	Inteligente (\(i^1\))	Alta (\(c^2\))	Bueno (\(e^1\))	Buena (\(r^1\))
\(\xi[2]\)	Difícil (\(d^1\))	No mucho (\(i^0\))	Media (\(c^1\))	Malo (\(e^0\))	Mala (\(r^0\))
\(\xi[3]\)	Fácil (\(d^0\))	No mucho (\(i^0\))	Baja (\(c^0\))	Malo (\(e^0\))	Mala (\(r^0\))
…	…	…	…	…	…
\(\xi[M]\)	Difícil (\(d^1\))	Inteligente (\(i^1\))	Alta (\(c^2\))	Bueno (\(e^1\))	Buena (\(r^1\))

Muestra (\(m\))	\(R[m]\)	\(\mathbf{I}\{R[m] = r^1\}\)
1	\(r^0\)	0
2	\(r^1\)	1
3	\(r^1\)	1
4	\(r^0\)	0
5	\(r^1\)	1

Muestra (\(m\))	\(C[m]\)	\(R[m]\)	\(\mathbf{I}\{R[m] = r^1\}\)
1	\(c^2\)	\(r^1\)	1
2	\(c^0\)	\(r^0\)	—
3	\(c^2\)	\(r^1\)	1
4	\(c^1\)	\(r^0\)	—
5	\(c^2\)	\(r^1\)	1
6	\(c^0\)	\(r^0\)	—
7	\(c^2\)	\(r^1\)	1
8	\(c^2\)	\(r^0\)	0

Paso	Descripción
(1) Orden topológico	Se ordenan las variables \(X_1, \dots, X_n\) de modo que los padres aparezcan antes que sus hijos. Esto garantiza que cuando llegamos a una variable \(X_i\), ya conocemos los valores de sus padres.
(2–3) Obtener valores de los padres	Para cada variable \(X_i\), se obtienen los valores actuales \(u_i\) de sus padres \(\text{Pa}(X_i)\), que ya fueron muestreados.
(4) Muestreo condicional	Se genera un valor \(x_i\) de acuerdo con la distribución condicional \(P(X_i \mid u_i)\). En la práctica, se genera un número aleatorio \(r \sim \text{Uniform}(0,1)\) y se selecciona el valor de \(X_i\) correspondiente al intervalo acumulado que contiene \(r\).
(5) Retorno de la partícula	El algoritmo devuelve la asignación completa \((x_1, \dots, x_n)\), que constituye una partícula del modelo.