Split View: 확률적 머신러닝 & 베이지안 방법론: 베이즈 추론부터 가우시안 프로세스, 불확실성 정량화까지

import numpy as np
import scipy.stats as stats

# Key probability distributions
x = np.linspace(-4, 4, 200)

# Gaussian distribution
gaussian = stats.norm(loc=0, scale=1)
print(f"Gaussian mean: {gaussian.mean():.2f}, variance: {gaussian.var():.2f}")

# Bernoulli distribution
p = 0.7
bern_samples = stats.bernoulli(p).rvs(1000)
print(f"Bernoulli success rate: {bern_samples.mean():.3f} (expected: {p})")

# Dirichlet distribution sampling
alpha = np.array([2.0, 3.0, 5.0])
dirichlet_samples = stats.dirichlet(alpha).rvs(1000)
print(f"Dirichlet sample mean: {dirichlet_samples.mean(axis=0)}")
print(f"Theoretical mean: {alpha / alpha.sum()}")

import pymc as pm
import numpy as np
import arviz as az

# Bayesian linear regression with PyMC
np.random.seed(42)
X = np.random.randn(100)
true_alpha, true_beta, true_sigma = 1.5, 2.3, 0.5
y = true_alpha + true_beta * X + np.random.randn(100) * true_sigma

with pm.Model() as linear_model:
    # Prior distributions
    alpha = pm.Normal('alpha', mu=0, sigma=10)
    beta = pm.Normal('beta', mu=0, sigma=10)
    sigma = pm.HalfNormal('sigma', sigma=1)

    # Likelihood function
    mu = alpha + beta * X
    likelihood = pm.Normal('y', mu=mu, sigma=sigma, observed=y)

    # Sample posterior with NUTS sampler
    trace = pm.sample(2000, tune=1000, return_inferencedata=True, random_seed=42)

# Summarize results
summary = az.summary(trace, var_names=['alpha', 'beta', 'sigma'])
print(summary)
# Should converge near alpha ~ 1.5, beta ~ 2.3, sigma ~ 0.5

import numpy as np

def metropolis_hastings(log_posterior, initial, n_samples=10000, proposal_std=0.5):
    """Metropolis-Hastings MCMC implementation"""
    samples = np.zeros((n_samples, len(initial)))
    current = np.array(initial, dtype=float)
    accepted = 0

    for i in range(n_samples):
        # Sample candidate from proposal distribution
        proposal = current + np.random.randn(len(current)) * proposal_std

        # Compute acceptance probability (in log space)
        log_ratio = log_posterior(proposal) - log_posterior(current)
        accept_prob = min(1.0, np.exp(log_ratio))

        # Accept/reject decision
        if np.random.rand() < accept_prob:
            current = proposal
            accepted += 1

        samples[i] = current

    print(f"Acceptance rate: {accepted / n_samples:.3f}")
    return samples

# Sample from a Gaussian mixture distribution
def log_posterior_mixture(theta):
    log_p1 = -0.5 * ((theta[0] - 2)**2 + (theta[1] - 2)**2)
    log_p2 = -0.5 * ((theta[0] + 2)**2 + (theta[1] + 2)**2)
    return np.log(0.5 * np.exp(log_p1) + 0.5 * np.exp(log_p2) + 1e-10)

samples = metropolis_hastings(log_posterior_mixture, [0.0, 0.0], n_samples=50000)

# Analyze after removing burn-in
burn_in = 5000
posterior_samples = samples[burn_in:]
print(f"Posterior mean: {posterior_samples.mean(axis=0)}")
print(f"Posterior std: {posterior_samples.std(axis=0)}")

import torch
import torch.nn as nn

class VAE(nn.Module):
    def __init__(self, input_dim=784, latent_dim=20, hidden_dim=400):
        super().__init__()
        # Encoder
        self.encoder = nn.Sequential(
            nn.Linear(input_dim, hidden_dim),
            nn.ReLU()
        )
        self.fc_mu = nn.Linear(hidden_dim, latent_dim)
        self.fc_logvar = nn.Linear(hidden_dim, latent_dim)

        # Decoder
        self.decoder = nn.Sequential(
            nn.Linear(latent_dim, hidden_dim),
            nn.ReLU(),
            nn.Linear(hidden_dim, input_dim),
            nn.Sigmoid()
        )

    def encode(self, x):
        h = self.encoder(x)
        return self.fc_mu(h), self.fc_logvar(h)

    def reparameterize(self, mu, logvar):
        """Reparameterization Trick: z = mu + eps * std"""
        std = torch.exp(0.5 * logvar)
        eps = torch.randn_like(std)
        return mu + eps * std

    def forward(self, x):
        mu, logvar = self.encode(x)
        z = self.reparameterize(mu, logvar)
        recon = self.decoder(z)
        return recon, mu, logvar

def vae_loss(recon_x, x, mu, logvar):
    # Reconstruction loss (BCE)
    recon_loss = nn.functional.binary_cross_entropy(recon_x, x, reduction='sum')
    # KL divergence: -0.5 * sum(1 + log(sigma^2) - mu^2 - sigma^2)
    kl_div = -0.5 * torch.sum(1 + logvar - mu.pow(2) - logvar.exp())
    return recon_loss + kl_div

import torch
import gpytorch

class ExactGPModel(gpytorch.models.ExactGP):
    def __init__(self, train_x, train_y, likelihood):
        super().__init__(train_x, train_y, likelihood)
        self.mean_module = gpytorch.means.ConstantMean()
        self.covar_module = gpytorch.kernels.ScaleKernel(
            gpytorch.kernels.RBFKernel()
        )

    def forward(self, x):
        mean_x = self.mean_module(x)
        covar_x = self.covar_module(x)
        return gpytorch.distributions.MultivariateNormal(mean_x, covar_x)

# Prepare data
torch.manual_seed(42)
train_x = torch.linspace(0, 1, 20)
train_y = torch.sin(train_x * 2 * 3.14159) + torch.randn(20) * 0.1

likelihood = gpytorch.likelihoods.GaussianLikelihood()
model = ExactGPModel(train_x, train_y, likelihood)

# Optimize hyperparameters (maximize MLL)
model.train()
likelihood.train()
optimizer = torch.optim.Adam(model.parameters(), lr=0.1)
mll = gpytorch.mlls.ExactMarginalLogLikelihood(likelihood, model)

for i in range(200):
    optimizer.zero_grad()
    output = model(train_x)
    loss = -mll(output, train_y)
    loss.backward()
    optimizer.step()

# Predict
model.eval()
likelihood.eval()
test_x = torch.linspace(0, 1, 100)
with torch.no_grad(), gpytorch.settings.fast_pred_var():
    predictions = likelihood(model(test_x))
    mean = predictions.mean
    lower, upper = predictions.confidence_region()

print(f"Learned lengthscale: {model.covar_module.base_kernel.lengthscale.item():.4f}")
print(f"Learned noise: {likelihood.noise.item():.4f}")

import torch
import torch.nn as nn
import numpy as np

class BayesianMLP(nn.Module):
    def __init__(self, input_dim, hidden_dim, output_dim, dropout_p=0.1):
        super().__init__()
        self.dropout_p = dropout_p
        self.net = nn.Sequential(
            nn.Linear(input_dim, hidden_dim),
            nn.ReLU(),
            nn.Dropout(p=dropout_p),
            nn.Linear(hidden_dim, hidden_dim),
            nn.ReLU(),
            nn.Dropout(p=dropout_p),
            nn.Linear(hidden_dim, output_dim)
        )

    def forward(self, x):
        return self.net(x)

def mc_dropout_predict(model, x, n_samples=100):
    """Estimate prediction uncertainty using MC Dropout"""
    model.train()  # Keep dropout active
    predictions = []
    with torch.no_grad():
        for _ in range(n_samples):
            pred = model(x)
            predictions.append(pred)

    predictions = torch.stack(predictions)  # (n_samples, batch, output)
    mean = predictions.mean(dim=0)
    # Uncertainty = variance of predictions
    uncertainty = predictions.var(dim=0)
    return mean, uncertainty

# Usage example
model = BayesianMLP(input_dim=10, hidden_dim=64, output_dim=1, dropout_p=0.1)
x_test = torch.randn(32, 10)
mean_pred, uncertainty = mc_dropout_predict(model, x_test, n_samples=200)
print(f"Prediction mean shape: {mean_pred.shape}, uncertainty: {uncertainty.shape}")
print(f"Mean uncertainty: {uncertainty.mean().item():.4f}")

import torch
import pyro
import pyro.distributions as dist
from pyro.nn import PyroModule, PyroSample

class BayesianLinear(PyroModule):
    def __init__(self, in_features, out_features):
        super().__init__()
        self.in_features = in_features
        self.out_features = out_features

        # Define prior distributions on weights and biases
        self.weight = PyroSample(
            dist.Normal(0., 1.).expand([out_features, in_features]).to_event(2)
        )
        self.bias = PyroSample(
            dist.Normal(0., 10.).expand([out_features]).to_event(1)
        )

    def forward(self, x):
        return x @ self.weight.T + self.bias

class BayesianNN(PyroModule):
    def __init__(self, in_dim, hidden_dim, out_dim):
        super().__init__()
        self.layer1 = BayesianLinear(in_dim, hidden_dim)
        self.layer2 = BayesianLinear(hidden_dim, out_dim)

    def forward(self, x, y=None):
        x = torch.relu(self.layer1(x))
        mu = self.layer2(x).squeeze(-1)
        sigma = pyro.sample("sigma", dist.Uniform(0., 1.))
        with pyro.plate("data", x.shape[0]):
            obs = pyro.sample("obs", dist.Normal(mu, sigma), obs=y)
        return mu

import jax
import jax.numpy as jnp
import numpyro
import numpyro.distributions as dist
from numpyro.infer import MCMC, NUTS

def logistic_regression_model(X, y=None):
    """NumPyro Bayesian logistic regression model"""
    n_features = X.shape[1]

    # Prior distributions
    alpha = numpyro.sample("alpha", dist.Normal(0., 10.))
    beta = numpyro.sample("beta", dist.Normal(jnp.zeros(n_features), jnp.ones(n_features)))

    # Compute logits
    logits = alpha + X @ beta

    # Likelihood
    with numpyro.plate("obs", X.shape[0]):
        numpyro.sample("y", dist.Bernoulli(logits=logits), obs=y)

# Run inference with NUTS sampler
def run_mcmc(X_train, y_train, num_samples=2000, num_warmup=1000):
    kernel = NUTS(logistic_regression_model)
    mcmc = MCMC(kernel, num_warmup=num_warmup, num_samples=num_samples)
    rng_key = jax.random.PRNGKey(42)
    mcmc.run(rng_key, X_train, y_train)
    return mcmc.get_samples()

# Example usage (generate data)
from sklearn.datasets import make_classification
from sklearn.preprocessing import StandardScaler

X, y = make_classification(n_samples=200, n_features=4, random_state=42)
X = StandardScaler().fit_transform(X)
X_jax = jnp.array(X)
y_jax = jnp.array(y)

samples = run_mcmc(X_jax, y_jax)
print(f"Inferred alpha: {samples['alpha'].mean():.3f} +/- {samples['alpha'].std():.3f}")
print(f"Inferred beta: {samples['beta'].mean(axis=0)}")

import pymc as pm
import numpy as np

# Hierarchical Bayesian model
# Model data from multiple groups simultaneously
n_groups = 5
n_per_group = 20
np.random.seed(42)

group_means = np.random.randn(n_groups) * 2
y_obs = np.concatenate([
    np.random.randn(n_per_group) + group_means[g]
    for g in range(n_groups)
])
group_idx = np.repeat(np.arange(n_groups), n_per_group)

with pm.Model() as hierarchical_model:
    # Hyperprior distributions
    mu_global = pm.Normal('mu_global', mu=0, sigma=10)
    sigma_global = pm.HalfNormal('sigma_global', sigma=5)

    # Group-level means (Partial Pooling)
    mu_group = pm.Normal('mu_group', mu=mu_global, sigma=sigma_global, shape=n_groups)

    # Observation likelihood
    sigma_obs = pm.HalfNormal('sigma_obs', sigma=2)
    y = pm.Normal('y', mu=mu_group[group_idx], sigma=sigma_obs, observed=y_obs)

    trace = pm.sample(2000, tune=1000, return_inferencedata=True, random_seed=42)

import arviz as az
print(az.summary(trace, var_names=['mu_group', 'mu_global', 'sigma_global']))