Arquivo de Python - Página 9 de 31

Inferência Aproximada (MCMC, Gibbs Sampling)

16/07/202618/06/2026 Por antonino

O que é inferência aproximada (MCMC e Gibbs)?

Inferência aproximada por MCMC é uma técnica para amostrar de distribuições complexas sem fórmula fechada. MCMC significa Markov Chain Monte Carlo, que gera uma cadeia de amostras dependentes. A cadeia é construída para convergir para a distribuição alvo (posteriori). Gibbs sampling é um caso especial do MCMC onde cada variável é amostrada condicionalmente. A cada passo, atualiza-se uma variável por vez, usando as demais fixas. Isso é feito quando a distribuição condicional completa é conhecida. Diferentemente da inferência exata, MCMC é escalável para redes com muitas variáveis. Ele não sofre com explosão combinatória, mas produz amostras correlacionadas. Portanto, ele é amplamente usado em problemas de alta dimensão.

Características fundamentais do MCMC

O MCMC possui três características principais que o definem. Primeiro, ele é assintoticamente exato: converge para a verdadeira posteriori com muitas amostras. Segundo, as amostras iniciais (burn-in) são descartadas para eliminar o efeito do início. Terceiro, a autocorrelação entre amostras reduz a eficiência, exigindo thinning. A convergência é monitorada por diagnósticos como o fator de escala de Gelman-Rubin. O MCMC não requer conhecimento da constante de normalização da posteriori. Isso é uma grande vantagem sobre métodos analíticos.

Vantagens e aplicações típicas

A principal vantagem é a capacidade de lidar com modelos hierárquicos e não-lineares. Ele é usado em aprendizado de máquina, econometria e genética populacional. Também é aplicado em análise de imagens e processamento de sinais. Contudo, a convergência pode ser lenta e difícil de diagnosticar.

O algoritmo de Metropolis-Hastings é o MCMC mais geral. Ele propõe um novo estado a partir de uma distribuição de proposta. Aceita-se o novo estado com probabilidade dada pela razão de verossimilhanças. Gibbs sampling é mais eficiente quando as condicionais são fáceis de amostrar. Por exemplo, em modelos com prioris conjugados, as condicionais são distribuições conhecidas. A cadeia de Gibbs amostra cada variável sequencialmente, mantendo as outras fixas. Após um número suficiente de iterações, as amostras representam a posteriori. A média das amostras estima a esperança a posteriori. Intervalos de credibilidade são obtidos dos quantis empíricos das amostras. O MCMC pode ser implementado em bibliotecas como PyMC3 e Stan. Essas ferramentas automatizam a escolha do amostrador e diagnósticos. O usuário só precisa definir o modelo e os dados. A inferência aproximada é a espinha dorsal da estatística bayesiana moderna. Assim, o MCMC é uma ferramenta indispensável para problemas reais.

Um exemplo clássico é estimar a média e variância de uma distribuição normal com dados observados. A posteriori conjunta não tem forma simples, mas as condicionais são normais e gamma inversa. Gibbs sampling amostra média dado variância, e variância dado média, alternadamente. Após algumas milhares de iterações, as amostras convergem para a posteriori verdadeira.

Enunciado do exemplo clássico

Implemente Gibbs sampling para estimar os parâmetros (μ, σ²) de uma normal com dados sintéticos. Gerar 50 dados de N(μ=5, σ²=4). Use prioris: μ ~ N(0, 1000) e σ² ~ Inverse-Gamma(3, 2). As condicionais completas são: μ | dados, σ² ~ N( média ponderada, variância ) e σ² | dados, μ ~ Inverse-Gamma( α + n/2, β + 0.5*Σ(xᵢ-μ)² ). Execute 5000 iterações, descarte os primeiros 1000 (burn-in). Plote as cadeias e os histogramas das amostras para μ e σ².

import numpy as np  
import matplotlib.pyplot as plt  
from scipy.stats import norm, invgamma  

# Gerar dados sintéticos  
np.random.seed(42)  
n = 50  
mu_real = 5.0  
sigma2_real = 4.0  
dados = np.random.normal(mu_real, np.sqrt(sigma2_real), n)  

# Prioris:  
# μ ~ N(0, 1000)  
mu_prior_mean = 0.0  
mu_prior_var = 1000.0  
# σ² ~ Inv-Gamma(alpha=3, beta=2)  
alpha_prior = 3.0  
beta_prior = 2.0  

# Gibbs sampling  
n_iter = 5000  
burn_in = 1000  
amostras_mu = np.zeros(n_iter)  
amostras_sigma2 = np.zeros(n_iter)  

# Inicialização  
mu_atual = 0.0  
sigma2_atual = 1.0  

for i in range(n_iter):  
    # Amostrar μ | dados, σ²  
    # Variância condicional  
    var_mu = 1.0 / (1.0/mu_prior_var + n/sigma2_atual)  
    # Média condicional  
    mean_mu = var_mu * (mu_prior_mean/mu_prior_var + np.sum(dados)/sigma2_atual)  
    mu_atual = np.random.normal(mean_mu, np.sqrt(var_mu))  

    # Amostrar σ² | dados, μ  
    alpha_post = alpha_prior + n/2.0  
    beta_post = beta_prior + 0.5 * np.sum((dados - mu_atual)**2)  
    sigma2_atual = invgamma.rvs(alpha_post, scale=beta_post)  

    amostras_mu[i] = mu_atual  
    amostras_sigma2[i] = sigma2_atual  

# Descartar burn-in  
mu_amostras = amostras_mu[burn_in:]  
sigma2_amostras = amostras_sigma2[burn_in:]  

# Estatísticas  
media_mu = np.mean(mu_amostras)  
ic_mu = np.percentile(mu_amostras, [2.5, 97.5])  
media_sigma2 = np.mean(sigma2_amostras)  
ic_sigma2 = np.percentile(sigma2_amostras, [2.5, 97.5])  

print(f"μ real = {mu_real}, estimado = {media_mu:.3f} (IC 95%: [{ic_mu[0]:.3f}, {ic_mu[1]:.3f}])")  
print(f"σ² real = {sigma2_real}, estimado = {media_sigma2:.3f} (IC 95%: [{ic_sigma2[0]:.3f}, {ic_sigma2[1]:.3f}])")  

# Gráficos  
plt.figure(figsize=(12, 10))  

# Cadeias de μ  
plt.subplot(2, 2, 1)  
plt.plot(mu_amostras, 'b-', linewidth=0.5, alpha=0.7)  
plt.axhline(mu_real, color='red', linestyle='--', label=f'μ real = {mu_real}')  
plt.axhline(media_mu, color='green', linestyle='-', label=f'Média amostral = {media_mu:.3f}')  
plt.title('Cadeia de Gibbs - μ')  
plt.xlabel('Iteração (após burn-in)')  
plt.ylabel('μ')  
plt.legend()  
plt.grid(True)  

# Histograma de μ  
plt.subplot(2, 2, 2)  
plt.hist(mu_amostras, bins=40, density=True, color='blue', alpha=0.6, edgecolor='black')  
x = np.linspace(mu_real-2, mu_real+2, 200)  
# Aproximar posteriori por normal (teórica)  
plt.plot(x, norm.pdf(x, media_mu, np.std(mu_amostras)), 'r-', linewidth=2, label='Aprox. normal')  
plt.axvline(mu_real, color='red', linestyle='--', label=f'μ real')  
plt.title('Distribuição posteriori de μ')  
plt.xlabel('μ')  
plt.ylabel('Densidade')  
plt.legend()  
plt.grid(True)  

# Cadeias de σ²  
plt.subplot(2, 2, 3)  
plt.plot(sigma2_amostras, 'g-', linewidth=0.5, alpha=0.7)  
plt.axhline(sigma2_real, color='red', linestyle='--', label=f'σ² real = {sigma2_real}')  
plt.axhline(media_sigma2, color='darkgreen', linestyle='-', label=f'Média amostral = {media_sigma2:.3f}')  
plt.title('Cadeia de Gibbs - σ²')  
plt.xlabel('Iteração (após burn-in)')  
plt.ylabel('σ²')  
plt.legend()  
plt.grid(True)  

# Histograma de σ²  
plt.subplot(2, 2, 4)  
plt.hist(sigma2_amostras, bins=40, density=True, color='green', alpha=0.6, edgecolor='black')  
x2 = np.linspace(0, 10, 200)  
# Aproximação por Inv-Gamma (teórica)  
plt.plot(x2, invgamma.pdf(x2, alpha_prior + n/2, scale=beta_prior + 0.5*np.var(dados)*n),  
         'r-', linewidth=2, label='Aprox. Inv-Gamma')  
plt.axvline(sigma2_real, color='red', linestyle='--', label=f'σ² real')  
plt.title('Distribuição posteriori de σ²')  
plt.xlabel('σ²')  
plt.ylabel('Densidade')  
plt.legend()  
plt.grid(True)  

plt.tight_layout()  
plt.show()

100

101

102

103

104

105

106

107

108

109

110

111

112

import numpy as np

import matplotlib.pyplot as plt

from scipy.stats import norm, invgamma

# Gerar dados sintéticos

np.random.seed(42)

n = 50

mu_real = 5.0

sigma2_real = 4.0

dados = np.random.normal(mu_real, np.sqrt(sigma2_real), n)

# Prioris:

# μ ~ N(0, 1000)

mu_prior_mean = 0.0

mu_prior_var = 1000.0

# σ² ~ Inv-Gamma(alpha=3, beta=2)

alpha_prior = 3.0

beta_prior = 2.0

# Gibbs sampling

n_iter = 5000

burn_in = 1000

amostras_mu = np.zeros(n_iter)

amostras_sigma2 = np.zeros(n_iter)

# Inicialização

mu_atual = 0.0

sigma2_atual = 1.0

for i in range(n_iter):

# Amostrar μ | dados, σ²

# Variância condicional

var_mu = 1.0 / (1.0/mu_prior_var + n/sigma2_atual)

# Média condicional

mean_mu = var_mu * (mu_prior_mean/mu_prior_var + np.sum(dados)/sigma2_atual)

mu_atual = np.random.normal(mean_mu, np.sqrt(var_mu))

# Amostrar σ² | dados, μ

alpha_post = alpha_prior + n/2.0

beta_post = beta_prior + 0.5 * np.sum((dados - mu_atual)**2)

sigma2_atual = invgamma.rvs(alpha_post, scale=beta_post)

amostras_mu[i] = mu_atual

amostras_sigma2[i] = sigma2_atual

# Descartar burn-in

mu_amostras = amostras_mu[burn_in:]

sigma2_amostras = amostras_sigma2[burn_in:]

# Estatísticas

media_mu = np.mean(mu_amostras)

ic_mu = np.percentile(mu_amostras, [2.5, 97.5])

media_sigma2 = np.mean(sigma2_amostras)

ic_sigma2 = np.percentile(sigma2_amostras, [2.5, 97.5])

print(f"μ real = {mu_real}, estimado = {media_mu:.3f} (IC 95%: [{ic_mu[0]:.3f}, {ic_mu[1]:.3f}])")

print(f"σ² real = {sigma2_real}, estimado = {media_sigma2:.3f} (IC 95%: [{ic_sigma2[0]:.3f}, {ic_sigma2[1]:.3f}])")

# Gráficos

plt.figure(figsize=(12, 10))

# Cadeias de μ

plt.subplot(2, 2, 1)

plt.plot(mu_amostras, 'b-', linewidth=0.5, alpha=0.7)

plt.axhline(mu_real, color='red', linestyle='--', label=f'μ real = {mu_real}')

plt.axhline(media_mu, color='green', linestyle='-', label=f'Média amostral = {media_mu:.3f}')

plt.title('Cadeia de Gibbs - μ')

plt.xlabel('Iteração (após burn-in)')

plt.ylabel('μ')

plt.legend()

plt.grid(True)

# Histograma de μ

plt.subplot(2, 2, 2)

plt.hist(mu_amostras, bins=40, density=True, color='blue', alpha=0.6, edgecolor='black')

x = np.linspace(mu_real-2, mu_real+2, 200)

# Aproximar posteriori por normal (teórica)

plt.plot(x, norm.pdf(x, media_mu, np.std(mu_amostras)), 'r-', linewidth=2, label='Aprox. normal')

plt.axvline(mu_real, color='red', linestyle='--', label=f'μ real')

plt.title('Distribuição posteriori de μ')

plt.xlabel('μ')

plt.ylabel('Densidade')

plt.legend()

plt.grid(True)

# Cadeias de σ²

plt.subplot(2, 2, 3)

plt.plot(sigma2_amostras, 'g-', linewidth=0.5, alpha=0.7)

plt.axhline(sigma2_real, color='red', linestyle='--', label=f'σ² real = {sigma2_real}')

plt.axhline(media_sigma2, color='darkgreen', linestyle='-', label=f'Média amostral = {media_sigma2:.3f}')

plt.title('Cadeia de Gibbs - σ²')

plt.xlabel('Iteração (após burn-in)')

plt.ylabel('σ²')

plt.legend()

plt.grid(True)

# Histograma de σ²

plt.subplot(2, 2, 4)

plt.hist(sigma2_amostras, bins=40, density=True, color='green', alpha=0.6, edgecolor='black')

x2 = np.linspace(0, 10, 200)

# Aproximação por Inv-Gamma (teórica)

plt.plot(x2, invgamma.pdf(x2, alpha_prior + n/2, scale=beta_prior + 0.5*np.var(dados)*n),

'r-', linewidth=2, label='Aprox. Inv-Gamma')

plt.axvline(sigma2_real, color='red', linestyle='--', label=f'σ² real')

plt.title('Distribuição posteriori de σ²')

plt.xlabel('σ²')

plt.ylabel('Densidade')

plt.legend()

plt.grid(True)

plt.tight_layout()

plt.show()

Este código implementa Gibbs sampling com prioris conjugados. As cadeias mostram a convergência rápida para os valores reais. Os histogramas exibem as distribuições posteriores aproximadas. Os intervalos de credibilidade contêm os valores verdadeiros. Para iniciantes, este exemplo demonstra o poder do MCMC prático. A inferência aproximada é, portanto, essencial para modelos complexos.

Inferência Exata (Eliminação de Variáveis)

16/07/202618/06/2026 Por antonino

O que é inferência exata por eliminação de variáveis?

Inferência exata por eliminação de variáveis é um algoritmo para calcular probabilidades marginais em redes bayesianas. Ele elimina variáveis não observadas (e não consultadas) uma por uma, somando-as. O processo transforma a distribuição conjunta em uma expressão com somatórios aninhados. Cada eliminação gera um novo fator (função) que captura a influência da variável removida. Esses fatores são combinados com os fatores restantes antes de cada eliminação. A ordem de eliminação afeta drasticamente a complexidade computacional. Uma ordem ruim pode criar fatores com muitas variáveis (explosão exponencial). Porém, com uma ordem ótima, o algoritmo é polinomial para redes de árvore. A eliminação de variáveis é a base de algoritmos mais avançados, como árvores de junção.

Características fundamentais

A eliminação de variáveis possui três características principais que a definem. Primeiro, ela é exata: não há aproximações, apenas manipulações algébricas. Segundo, ela usa a propriedade de fatoração da distribuição conjunta. Terceiro, a complexidade é exponencial no tamanho do maior fator gerado (treewidth). O algoritmo é completo para redes com treewidth pequeno. Além disso, ele pode ser usado para calcular probabilidades condicionais e evidências. A ordem de eliminação pode ser determinada por heurísticas (ex.: menor grau).

Vantagens e limitações

A principal vantagem é a garantia de resultado exato para redes moderadas. Ela é usada em diagnósticos, sistemas especialistas e validação de modelos. Também é útil para comparar com métodos aproximados (MCMC). Contudo, para redes com treewidth alto (ex.: redes densas), ela é inviável. Nesses casos, usa-se inferência aproximada por amostragem ou variational.

O algoritmo funciona fatorando a conjunta em fatores locais (CPTs). Cada fator é uma tabela sobre um subconjunto de variáveis. Para eliminar uma variável X, multiplicam-se todos os fatores que contêm X. Em seguida, soma-se X dessa tabela, criando um novo fator sem X. Esse novo fator é adicionado ao conjunto de fatores restantes. O processo repete até que apenas a(s) variável(eis) de interesse permaneçam. Finalmente, normaliza-se o fator resultante para obter a distribuição condicional. Se houver evidência, os fatores são condicionados (fixando os valores). A ordem de eliminação pode ser escolhida para minimizar o custo. Heurísticas como “elimine a variável com menor número de fatores” são comuns. A eliminação de variáveis é equivalente a reordenar os somatórios na conjunta. Ela é ensinada como primeiro passo para entender inferência em grafos. Para redes pequenas, é implementada facilmente com dicionários. Assim, a eliminação de variáveis é um algoritmo fundamental e didático.

Um exemplo clássico é uma rede com três variáveis: A → B → C. Queremos P(C | A) sem evidências. Eliminamos B somando sobre seus valores. A conjunta é P(A)P(B|A)P(C|B). Somando B, obtemos P(A)P(C|A). Isso é feito em um passo, com fator resultante envolvendo A e C.

Enunciado do exemplo clássico

Implemente a eliminação de variáveis para a rede: Variáveis: A, B, C, D (todas binárias). Relações: A → B, A → C, B → D, C → D. CPTs fornecidas manualmente (valores arbitrários). Calcule P(D | A=1) e P(C | D=0, B=1) usando eliminação. Plote a estrutura da rede e os fatores gerados durante a eliminação.

import numpy as np  
import matplotlib.pyplot as plt  
import networkx as nx  
from itertools import product  

# Definição da rede (todas binárias: 0/1)  
# CPTs como dicionários: P(var | pais)  
# P(A)  
P_A = {0: 0.6, 1: 0.4}  

# P(B | A)  
P_B_dado_A = {  
    (0,0): 0.7, (0,1): 0.3,  # A=0, B=0/1  
    (1,0): 0.2, (1,1): 0.8   # A=1, B=0/1  
}  

# P(C | A)  
P_C_dado_A = {  
    (0,0): 0.6, (0,1): 0.4,  
    (1,0): 0.1, (1,1): 0.9  
}  

# P(D | B, C)  
P_D_dado_B_C = {  
    (0,0,0): 0.9, (0,0,1): 0.1,  # B=0,C=0, D=0/1  
    (0,1,0): 0.4, (0,1,1): 0.6,  # B=0,C=1  
    (1,0,0): 0.3, (1,0,1): 0.7,  # B=1,C=0  
    (1,1,0): 0.2, (1,1,1): 0.8   # B=1,C=1  
}  

# Representar fatores como dicionários {tupla_de_variaveis: valor}  
# Fator é uma função de (var1, var2, ...) -> probabilidade  

def fator_prior(var, dist):  
    f = {}  
    for v in dist:  
        f[(v,)] = dist[v]  
    return f  

def fator_condicional(var, pais, dist):  
    f = {}  
    for chave, prob in dist.items():  
        # chave é (pai1, pai2, ..., var_valor) no formato definido  
        # Precisamos extrair os valores na ordem correta  
        # Vamos assumir que a ordem é (pais..., var)  
        f[tuple(list(chave[:-1]) + [chave[-1]])] = prob  
    return f  

# Construir fatores iniciais  
fatores = []  
# Fator A  
fatores.append({(a,): P_A[a] for a in (0,1)})  
# Fator B|A  
f_B = {}  
for a in (0,1):  
    for b in (0,1):  
        f_B[(a,b)] = P_B_dado_A[(a,b)]  
fatores.append(f_B)  
# Fator C|A  
f_C = {}  
for a in (0,1):  
    for c in (0,1):  
        f_C[(a,c)] = P_C_dado_A[(a,c)]  
fatores.append(f_C)  
# Fator D|B,C  
f_D = {}  
for b in (0,1):  
    for c in (0,1):  
        for d in (0,1):  
            f_D[(b,c,d)] = P_D_dado_B_C[(b,c,d)]  
fatores.append(f_D)  

# Função para multiplicar fatores que compartilham variáveis  
def multiplicar_fatores(f1, f2):  
    # f1 e f2: dict com chaves = tupla de variáveis  
    # Retorna novo fator com união das variáveis  
    vars1 = set(f1.keys()).pop() if len(f1)==1 else set(next(iter(f1.keys())))  
    # Na verdade, melhor: pegar a primeira chave para ver as variáveis  
    vars1 = set(next(iter(f1.keys())))  
    vars2 = set(next(iter(f2.keys())))  
    vars_union = tuple(sorted(vars1 | vars2))  
    novo = {}  
    # Gerar todas as combinações para as variáveis da união  
    # Como todas são binárias, usamos product  
    for vals in product([0,1], repeat=len(vars_union)):  
        chave = tuple(vals)  
        # Extrair valores para f1 e f2  
        sub1 = tuple(vals[vars_union.index(v)] for v in sorted(vars1))  
        sub2 = tuple(vals[vars_union.index(v)] for v in sorted(vars2))  
        # Multiplicar (se existir, senao zero)  
        val1 = f1.get(sub1, 0.0)  
        val2 = f2.get(sub2, 0.0)  
        novo[chave] = val1 * val2  
    return novo  

# Função para eliminar uma variável (somar)  
def eliminar_variavel(fatores, var):  
    # Encontrar todos os fatores que contêm a variável  
    fatores_contem = []  
    outros = []  
    for f in fatores:  
        vars_f = set(next(iter(f.keys())))  
        if var in vars_f:  
            fatores_contem.append(f)  
        else:  
            outros.append(f)  
    if not fatores_contem:  
        return outros  # variável não está presente  
    # Multiplicar todos esses fatores  
    f_prod = fatores_contem[0]  
    for f in fatores_contem[1:]:  
        f_prod = multiplicar_fatores(f_prod, f)  
    # Somar sobre a variável var  
    novo_fator = {}  
    for chave, val in f_prod.items():  
        # chave é uma tupla com var e outras variáveis  
        # criar nova chave sem a variável var  
        idx = list(chave).index(var)  
        nova_chave = tuple(chave[:idx] + chave[idx+1:])  
        novo_fator[nova_chave] = novo_fator.get(nova_chave, 0.0) + val  
    outros.append(novo_fator)  
    return outros  

# Consulta 1: P(D | A=1)  
# Evidência: A=1. Primeiro, condicionamos fixando A=1 em todos os fatores  
fatores_evid = []  
for f in fatores:  
    novo_f = {}  
    for chave, val in f.items():  
        # Se A está na chave, verificar se é 1, senão descartar  
        if 'A' in str(chave):  # forma simples de verificar  
            # encontrar posição de A  
            # como usamos nomes, vamos converter para lista de strings  
            # melhor: assumir ordem fixa: (A,B,C,D)  
            # Para simplificar, usamos índices: A é 0, B é 1, C é 2, D é 3  
            pass  
    # Vamos reimplementar de forma mais robusta usando índices  
# Devido à complexidade de manipular nomes, vamos usar índices (0=A,1=B,2=C,3=D)  
# Refazer fatores com índices  
fatores_idx = []  
# Fator A (índice 0)  
fatores_idx.append({(0,): P_A[0], (1,): P_A[1]})  # chave (valor de A,)  
# Fator B|A (A,B)  
f_B_idx = {}  
for a in (0,1):  
    for b in (0,1):  
        f_B_idx[(a,b)] = P_B_dado_A[(a,b)]  
fatores_idx.append(f_B_idx)  
# Fator C|A (A,C)  
f_C_idx = {}  
for a in (0,1):  
    for c in (0,1):  
        f_C_idx[(a,c)] = P_C_dado_A[(a,c)]  
fatores_idx.append(f_C_idx)  
# Fator D|B,C (B,C,D)  
f_D_idx = {}  
for b in (0,1):  
    for c in (0,1):  
        for d in (0,1):  
            f_D_idx[(b,c,d)] = P_D_dado_B_C[(b,c,d)]  
fatores_idx.append(f_D_idx)  

# Função de multiplicação para índices (reusar a anterior, mas com tuplas)  
def mult_fatores(f1, f2):  
    vars1 = set(next(iter(f1.keys())))  
    vars2 = set(next(iter(f2.keys())))  
    vars_union = tuple(sorted(vars1 | vars2))  
    novo = {}  
    for vals in product([0,1], repeat=len(vars_union)):  
        chave = tuple(vals)  
        sub1 = tuple(vals[vars_union.index(v)] for v in sorted(vars1))  
        sub2 = tuple(vals[vars_union.index(v)] for v in sorted(vars2))  
        val1 = f1.get(sub1, 0.0)  
        val2 = f2.get(sub2, 0.0)  
        novo[chave] = val1 * val2  
    return novo  

def eliminar_idx(fatores, var):  
    contem = []  
    outros = []  
    for f in fatores:  
        vars_f = set(next(iter(f.keys())))  
        if var in vars_f:  
            contem.append(f)  
        else:  
            outros.append(f)  
    if not contem:  
        return outros  
    f_prod = contem[0]  
    for f in contem[1:]:  
        f_prod = mult_fatores(f_prod, f)  
    novo = {}  
    for chave, val in f_prod.items():  
        idx = list(chave).index(var)  
        nova_chave = tuple(chave[:idx] + chave[idx+1:])  
        novo[nova_chave] = novo.get(nova_chave, 0.0) + val  
    outros.append(novo)  
    return outros  

# Consulta 1: P(D | A=1)  
# Condicionar: fixar A=1  
fatores_cond = []  
for f in fatores_idx:  
    novo = {}  
    for chave, val in f.items():  
        if 0 in chave:  # A está presente  
            idxA = list(chave).index(0)  
            if chave[idxA] == 1:  
                nova_chave = tuple(chave[:idxA] + chave[idxA+1:])  
                novo[nova_chave] = novo.get(nova_chave, 0.0) + val  
        else:  
            novo[chave] = val  
    fatores_cond.append(novo)  

# Agora eliminar B e C (variáveis não consultadas)  
fatores_temp = fatores_cond  
fatores_temp = eliminar_idx(fatores_temp, 1)  # eliminar B  
fatores_temp = eliminar_idx(fatores_temp, 2)  # eliminar C  
# Agora sobra fator sobre D (índice 3)  
fator_D = fatores_temp[0]  
# Normalizar para obter P(D | A=1)  
soma = sum(fator_D.values())  
if soma > 0:  
    fator_D = {k: v/soma for k,v in fator_D.items()}  
P_D_dado_A1 = fator_D  
print("P(D | A=1):")  
for d, prob in sorted(P_D_dado_A1.items()):  
    print(f"  P(D={d[0]}) = {prob:.4f}")  

# Consulta 2: P(C | D=0, B=1)  
# Evidência: D=0, B=1. Fixar esses valores.  
fatores_cond2 = []  
for f in fatores_idx:  
    novo = {}  
    for chave, val in f.items():  
        # Verificar se D (3) e B (1) estão na chave  
        ok = True  
        chave_lista = list(chave)  
        if 3 in chave:  
            idxD = chave.index(3)  
            if chave[idxD] != 0:  
                ok = False  
        if 1 in chave:  
            idxB = chave.index(1)  
            if chave[idxB] != 1:  
                ok = False  
        if ok:  
            # Remover B e D da chave  
            nova_chave = []  
            for i, v in enumerate(chave):  
                if i != idxB and i != idxD:  
                    nova_chave.append(v)  
            nova_chave = tuple(nova_chave)  
            novo[nova_chave] = novo.get(nova_chave, 0.0) + val  
    fatores_cond2.append(novo)  

# Eliminar A (índice 0)  
fatores_temp2 = eliminar_idx(fatores_cond2, 0)  
# Agora sobra fator sobre C (índice 2)  
fator_C = fatores_temp2[0]  
soma2 = sum(fator_C.values())  
if soma2 > 0:  
    fator_C = {k: v/soma2 for k,v in fator_C.items()}  
print("\nP(C | D=0, B=1):")  
for c, prob in sorted(fator_C.items()):  
    print(f"  P(C={c[0]}) = {prob:.4f}")  

# Plotar estrutura da rede  
G = nx.DiGraph()  
G.add_edges_from([('A','B'), ('A','C'), ('B','D'), ('C','D')])  
pos = {'A': (0, 1), 'B': (0, 0), 'C': (1, 1), 'D': (1, 0)}  

plt.figure(figsize=(8, 4))  
nx.draw(G, pos, with_labels=True, node_color='lightblue',  
        node_size=2500, font_size=12, font_weight='bold',  
        arrowsize=20, edge_color='gray', width=2)  
plt.title('Rede Bayesiana para Eliminação de Variáveis')  
plt.show()

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

133

134

135

136

137

138

139

140

141

142

143

144

145

146

147

148

149

150

151

152

153

154

155

156

157

158

159

160

161

162

163

164

165

166

167

168

169

170

171

172

173

174

175

176

177

178

179

180

181

182

183

184

185

186

187

188

189

190

191

192

193

194

195

196

197

198

199

200

201

202

203

204

205

206

207

208

209

210

211

212

213

214

215

216

217

218

219

220

221

222

223

224

225

226

227

228

229

230

231

232

233

234

235

236

237

238

239

240

241

242

243

244

245

246

247

248

249

250

251

252

253

254

255

256

257

258

259

260

261

262

263

264

265

266

267

268

269

270

271

272

273

274

275

276

277

278

import numpy as np

import matplotlib.pyplot as plt

import networkx as nx

from itertools import product

# Definição da rede (todas binárias: 0/1)

# CPTs como dicionários: P(var | pais)

# P(A)

P_A = {0: 0.6, 1: 0.4}

# P(B | A)

P_B_dado_A = {

(0,0): 0.7, (0,1): 0.3, # A=0, B=0/1

(1,0): 0.2, (1,1): 0.8 # A=1, B=0/1

}

# P(C | A)

P_C_dado_A = {

(0,0): 0.6, (0,1): 0.4,

(1,0): 0.1, (1,1): 0.9

}

# P(D | B, C)

P_D_dado_B_C = {

(0,0,0): 0.9, (0,0,1): 0.1, # B=0,C=0, D=0/1

(0,1,0): 0.4, (0,1,1): 0.6, # B=0,C=1

(1,0,0): 0.3, (1,0,1): 0.7, # B=1,C=0

(1,1,0): 0.2, (1,1,1): 0.8 # B=1,C=1

}

# Representar fatores como dicionários {tupla_de_variaveis: valor}

# Fator é uma função de (var1, var2, ...) -> probabilidade

def fator_prior(var, dist):

f = {}

for v in dist:

f[(v,)] = dist[v]

return f

def fator_condicional(var, pais, dist):

f = {}

for chave, prob in dist.items():

# chave é (pai1, pai2, ..., var_valor) no formato definido

# Precisamos extrair os valores na ordem correta

# Vamos assumir que a ordem é (pais..., var)

f[tuple(list(chave[:-1]) + [chave[-1]])] = prob

return f

# Construir fatores iniciais

fatores = []

# Fator A

fatores.append({(a,): P_A[a] for a in (0,1)})

# Fator B|A

f_B = {}

for a in (0,1):

for b in (0,1):

f_B[(a,b)] = P_B_dado_A[(a,b)]

fatores.append(f_B)

# Fator C|A

f_C = {}

for a in (0,1):

for c in (0,1):

f_C[(a,c)] = P_C_dado_A[(a,c)]

fatores.append(f_C)

# Fator D|B,C

f_D = {}

for b in (0,1):

for c in (0,1):

for d in (0,1):

f_D[(b,c,d)] = P_D_dado_B_C[(b,c,d)]

fatores.append(f_D)

# Função para multiplicar fatores que compartilham variáveis

def multiplicar_fatores(f1, f2):

# f1 e f2: dict com chaves = tupla de variáveis

# Retorna novo fator com união das variáveis

vars1 = set(f1.keys()).pop() if len(f1)==1 else set(next(iter(f1.keys())))

# Na verdade, melhor: pegar a primeira chave para ver as variáveis

vars1 = set(next(iter(f1.keys())))

vars2 = set(next(iter(f2.keys())))

vars_union = tuple(sorted(vars1 | vars2))

novo = {}

# Gerar todas as combinações para as variáveis da união

# Como todas são binárias, usamos product

for vals in product([0,1], repeat=len(vars_union)):

chave = tuple(vals)

# Extrair valores para f1 e f2

sub1 = tuple(vals[vars_union.index(v)] for v in sorted(vars1))

sub2 = tuple(vals[vars_union.index(v)] for v in sorted(vars2))

# Multiplicar (se existir, senao zero)

val1 = f1.get(sub1, 0.0)

val2 = f2.get(sub2, 0.0)

novo[chave] = val1 * val2

return novo

# Função para eliminar uma variável (somar)

def eliminar_variavel(fatores, var):

# Encontrar todos os fatores que contêm a variável

fatores_contem = []

outros = []

for f in fatores:

vars_f = set(next(iter(f.keys())))

if var in vars_f:

fatores_contem.append(f)

else:

outros.append(f)

if not fatores_contem:

return outros # variável não está presente

# Multiplicar todos esses fatores

f_prod = fatores_contem[0]

for f in fatores_contem[1:]:

f_prod = multiplicar_fatores(f_prod, f)

# Somar sobre a variável var

novo_fator = {}

for chave, val in f_prod.items():

# chave é uma tupla com var e outras variáveis

# criar nova chave sem a variável var

idx = list(chave).index(var)

nova_chave = tuple(chave[:idx] + chave[idx+1:])

novo_fator[nova_chave] = novo_fator.get(nova_chave, 0.0) + val

outros.append(novo_fator)

return outros

# Consulta 1: P(D | A=1)

# Evidência: A=1. Primeiro, condicionamos fixando A=1 em todos os fatores

fatores_evid = []

for f in fatores:

novo_f = {}

for chave, val in f.items():

# Se A está na chave, verificar se é 1, senão descartar

if 'A' in str(chave): # forma simples de verificar

# encontrar posição de A

# como usamos nomes, vamos converter para lista de strings

# melhor: assumir ordem fixa: (A,B,C,D)

# Para simplificar, usamos índices: A é 0, B é 1, C é 2, D é 3

pass

# Vamos reimplementar de forma mais robusta usando índices

# Devido à complexidade de manipular nomes, vamos usar índices (0=A,1=B,2=C,3=D)

# Refazer fatores com índices

fatores_idx = []

# Fator A (índice 0)

fatores_idx.append({(0,): P_A[0], (1,): P_A[1]}) # chave (valor de A,)

# Fator B|A (A,B)

f_B_idx = {}

for a in (0,1):

for b in (0,1):

f_B_idx[(a,b)] = P_B_dado_A[(a,b)]

fatores_idx.append(f_B_idx)

# Fator C|A (A,C)

f_C_idx = {}

for a in (0,1):

for c in (0,1):

f_C_idx[(a,c)] = P_C_dado_A[(a,c)]

fatores_idx.append(f_C_idx)

# Fator D|B,C (B,C,D)

f_D_idx = {}

for b in (0,1):

for c in (0,1):

for d in (0,1):

f_D_idx[(b,c,d)] = P_D_dado_B_C[(b,c,d)]

fatores_idx.append(f_D_idx)

# Função de multiplicação para índices (reusar a anterior, mas com tuplas)

def mult_fatores(f1, f2):

vars1 = set(next(iter(f1.keys())))

vars2 = set(next(iter(f2.keys())))

vars_union = tuple(sorted(vars1 | vars2))

novo = {}

for vals in product([0,1], repeat=len(vars_union)):

chave = tuple(vals)

sub1 = tuple(vals[vars_union.index(v)] for v in sorted(vars1))

sub2 = tuple(vals[vars_union.index(v)] for v in sorted(vars2))

val1 = f1.get(sub1, 0.0)

val2 = f2.get(sub2, 0.0)

novo[chave] = val1 * val2

return novo

def eliminar_idx(fatores, var):

contem = []

outros = []

for f in fatores:

vars_f = set(next(iter(f.keys())))

if var in vars_f:

contem.append(f)

else:

outros.append(f)

if not contem:

return outros

f_prod = contem[0]

for f in contem[1:]:

f_prod = mult_fatores(f_prod, f)

novo = {}

for chave, val in f_prod.items():

idx = list(chave).index(var)

nova_chave = tuple(chave[:idx] + chave[idx+1:])

novo[nova_chave] = novo.get(nova_chave, 0.0) + val

outros.append(novo)

return outros

# Consulta 1: P(D | A=1)

# Condicionar: fixar A=1

fatores_cond = []

for f in fatores_idx:

novo = {}

for chave, val in f.items():

if 0 in chave: # A está presente

idxA = list(chave).index(0)

if chave[idxA] == 1:

nova_chave = tuple(chave[:idxA] + chave[idxA+1:])

novo[nova_chave] = novo.get(nova_chave, 0.0) + val

else:

novo[chave] = val

fatores_cond.append(novo)

# Agora eliminar B e C (variáveis não consultadas)

fatores_temp = fatores_cond

fatores_temp = eliminar_idx(fatores_temp, 1) # eliminar B

fatores_temp = eliminar_idx(fatores_temp, 2) # eliminar C

# Agora sobra fator sobre D (índice 3)

fator_D = fatores_temp[0]

# Normalizar para obter P(D | A=1)

soma = sum(fator_D.values())

if soma > 0:

fator_D = {k: v/soma for k,v in fator_D.items()}

P_D_dado_A1 = fator_D

print("P(D | A=1):")

for d, prob in sorted(P_D_dado_A1.items()):

print(f" P(D={d[0]}) = {prob:.4f}")

# Consulta 2: P(C | D=0, B=1)

# Evidência: D=0, B=1. Fixar esses valores.

fatores_cond2 = []

for f in fatores_idx:

novo = {}

for chave, val in f.items():

# Verificar se D (3) e B (1) estão na chave

ok = True

chave_lista = list(chave)

if 3 in chave:

idxD = chave.index(3)

if chave[idxD] != 0:

ok = False

if 1 in chave:

idxB = chave.index(1)

if chave[idxB] != 1:

ok = False

if ok:

# Remover B e D da chave

nova_chave = []

for i, v in enumerate(chave):

if i != idxB and i != idxD:

nova_chave.append(v)

nova_chave = tuple(nova_chave)

novo[nova_chave] = novo.get(nova_chave, 0.0) + val

fatores_cond2.append(novo)

# Eliminar A (índice 0)

fatores_temp2 = eliminar_idx(fatores_cond2, 0)

# Agora sobra fator sobre C (índice 2)

fator_C = fatores_temp2[0]

soma2 = sum(fator_C.values())

if soma2 > 0:

fator_C = {k: v/soma2 for k,v in fator_C.items()}

print("\nP(C | D=0, B=1):")

for c, prob in sorted(fator_C.items()):

print(f" P(C={c[0]}) = {prob:.4f}")

# Plotar estrutura da rede

G = nx.DiGraph()

G.add_edges_from([('A','B'), ('A','C'), ('B','D'), ('C','D')])

pos = {'A': (0, 1), 'B': (0, 0), 'C': (1, 1), 'D': (1, 0)}

plt.figure(figsize=(8, 4))

nx.draw(G, pos, with_labels=True, node_color='lightblue',

node_size=2500, font_size=12, font_weight='bold',

arrowsize=20, edge_color='gray', width=2)

plt.title('Rede Bayesiana para Eliminação de Variáveis')

plt.show()

Este código implementa eliminação de variáveis com fatores e evidências. A estrutura da rede é exibida graficamente para referência. Os resultados mostram as probabilidades condicionais calculadas exatamente. A eliminação de B e C para a primeira consulta reduz a complexidade. Para iniciantes, este exemplo ilustra o mecanismo interno da inferência exata. A eliminação de variáveis é, portanto, um algoritmo fundamental e transparente.