📊 概率统计基础¶

import numpy as np

# 经典例子：疾病检测
# 疾病发病率1%，检测准确率99%，假阳性率2%
# 检测阳性时，真的有病的概率？

P_disease = 0.01
P_positive_given_disease = 0.99    # 灵敏度
P_positive_given_healthy = 0.02    # 假阳性率

# 贝叶斯公式
P_positive = P_positive_given_disease * P_disease + \
             P_positive_given_healthy * (1 - P_disease)
P_disease_given_positive = P_positive_given_disease * P_disease / P_positive

print(f"检测阳性时，实际有病概率: {P_disease_given_positive:.4f}")
# ≈ 0.3333 — 只有33%！这就是贝叶斯的反直觉之处

import matplotlib.pyplot as plt
from scipy import stats

fig, axes = plt.subplots(1, 3, figsize=(15, 4))
x = np.linspace(-5, 5, 200)

# 高斯分布族
for mu, sigma in [(0, 1), (0, 2), (2, 0.5)]:
    axes[0].plot(x, stats.norm.pdf(x, mu, sigma),
                 label=f'μ={mu}, σ={sigma}')
axes[0].legend(); axes[0].set_title("高斯分布")

# Beta分布（贝叶斯先验）
x_beta = np.linspace(0, 1, 200)
for a, b in [(1, 1), (5, 2), (2, 5), (0.5, 0.5)]:
    axes[1].plot(x_beta, stats.beta.pdf(x_beta, a, b),
                 label=f'α={a}, β={b}')
axes[1].legend(); axes[1].set_title("Beta分布(贝叶斯先验)")

# 泊松分布
x_pois = np.arange(0, 20)
for lam in [1, 4, 10]:
    axes[2].bar(x_pois + lam*0.1, stats.poisson.pmf(x_pois, lam),
                width=0.3, alpha=0.7, label=f'λ={lam}')
axes[2].legend(); axes[2].set_title("泊松分布")
plt.tight_layout(); plt.show()

# 手写高斯MLE
def gaussian_mle(data):
    mu_hat = np.mean(data)
    sigma2_hat = np.mean((data - mu_hat) ** 2)
    return mu_hat, sigma2_hat

data = np.random.normal(loc=5, scale=2, size=1000)
mu, sigma2 = gaussian_mle(data)
print(f"MLE估计: μ={mu:.4f}, σ²={sigma2:.4f}")
print(f"真实值:  μ=5.0000, σ²=4.0000")

class GaussianMixtureEM:
    """手写高斯混合模型EM算法"""

    def __init__(self, n_components=3, max_iter=100, tol=1e-6):  # __init__构造方法，创建对象时自动调用
        self.K = n_components
        self.max_iter = max_iter
        self.tol = tol

    def fit(self, X):
        n, d = X.shape

        # 初始化参数
        self.weights = np.ones(self.K) / self.K  # π_k
        self.means = X[np.random.choice(n, self.K, replace=False)]
        self.covs = [np.eye(d) for _ in range(self.K)]  # 列表推导式，简洁创建列表

        log_likelihood_old = -np.inf

        for iteration in range(self.max_iter):
            # === E步：计算责任度 γ(z_nk) ===
            responsibilities = np.zeros((n, self.K))
            for k in range(self.K):
                responsibilities[:, k] = self.weights[k] * \
                    self._gaussian_pdf(X, self.means[k], self.covs[k])
            responsibilities /= responsibilities.sum(axis=1, keepdims=True)

            # === M步：更新参数 ===
            N_k = responsibilities.sum(axis=0)  # 每个分量的有效样本数

            for k in range(self.K):
                # 更新均值
                self.means[k] = (responsibilities[:, k:k+1].T @ X) / N_k[k]
                # 更新协方差
                diff = X - self.means[k]
                self.covs[k] = (responsibilities[:, k:k+1] * diff).T @ diff / N_k[k]
                self.covs[k] += 1e-6 * np.eye(d)  # 数值稳定
            # 更新权重
            self.weights = N_k / n

            # 检查收敛
            log_likelihood = np.sum(np.log(
                sum(self.weights[k] * self._gaussian_pdf(X, self.means[k], self.covs[k])
                    for k in range(self.K))))
            if abs(log_likelihood - log_likelihood_old) < self.tol:
                print(f"EM在第{iteration}轮收敛")
                break
            log_likelihood_old = log_likelihood

    def _gaussian_pdf(self, X, mean, cov):
        d = X.shape[1]
        diff = X - mean
        inv_cov = np.linalg.inv(cov)  # np.linalg线性代数运算
        exponent = -0.5 * np.sum(diff @ inv_cov * diff, axis=1)
        norm = 1 / np.sqrt((2*np.pi)**d * np.linalg.det(cov))
        return norm * np.exp(exponent)

# 测试
np.random.seed(42)
X = np.vstack([
    np.random.randn(100, 2) + [0, 0],
    np.random.randn(100, 2) + [5, 5],
    np.random.randn(100, 2) + [0, 5],
])
gmm = GaussianMixtureEM(n_components=3)
gmm.fit(X)
print(f"估计的均值:\n{gmm.means}")
print(f"估计的权重: {gmm.weights}")

def entropy(probs):
    """计算离散分布的熵"""
    probs = np.array(probs)  # np.array创建NumPy数组
    probs = probs[probs > 0]  # 避免log(0)
    return -np.sum(probs * np.log(probs))

# 公平硬币 vs 不公平硬币
print(f"公平硬币熵: {entropy([0.5, 0.5]):.4f} nats")   # 0.6931
print(f"偏硬币熵:   {entropy([0.9, 0.1]):.4f} nats")   # 0.325
print(f"确定事件熵: {entropy([1.0, 0.0]):.4f} nats")   # 0.0

def cross_entropy(p, q):
    """交叉熵H(p, q)"""
    q = np.clip(q, 1e-10, 1)  # 数值稳定
    return -np.sum(p * np.log(q))

def kl_divergence(p, q):
    """KL散度D_KL(p||q)"""
    p, q = np.array(p), np.array(q)
    mask = p > 0
    return np.sum(p[mask] * np.log(p[mask] / q[mask]))

p = [0.7, 0.2, 0.1]  # 真实分布
q = [0.6, 0.3, 0.1]  # 模型分布

print(f"H(p):       {entropy(p):.4f}")
print(f"H(p, q):    {cross_entropy(p, q):.4f}")
print(f"D_KL(p||q): {kl_divergence(p, q):.4f}")
print(f"验证: H(p) + D_KL = {entropy(p) + kl_divergence(p, q):.4f}")  # 应等于H(p,q)

import torch

def vae_loss(x_recon, x, mu, log_var):
    """VAE损失 = 重构损失 + KL散度"""
    # 重构损失 = E_q[log p(x|z)]
    recon_loss = torch.nn.functional.binary_cross_entropy(
        x_recon, x, reduction='sum')

    # KL散度 = D_KL(q(z|x) || p(z))
    # 解析解推导见上：
    # D_KL = -0.5 * sum(1 + log_var - mu^2 - exp(log_var))
    kl_loss = -0.5 * torch.sum(1 + log_var - mu.pow(2) - log_var.exp())

    return recon_loss + kl_loss

import numpy as np
from scipy.stats import multivariate_normal

class GaussianMixtureEM:
    """高斯混合模型的EM实现"""
    def __init__(self, n_components=3, max_iter=100, tol=1e-6):
        self.K = n_components
        self.max_iter = max_iter
        self.tol = tol

    def fit(self, X):
        n, d = X.shape

        # 初始化
        np.random.seed(42)
        self.pi = np.ones(self.K) / self.K  # 混合权重
        indices = np.random.choice(n, self.K, replace=False)
        self.mu = X[indices].copy()  # 均值
        self.sigma = [np.eye(d) for _ in range(self.K)]  # 协方差（不能用[]*K，会共享引用）

        log_likelihoods = []

        for iteration in range(self.max_iter):
            # === E步：计算责任度 γ(z_nk) ===
            gamma = np.zeros((n, self.K))
            for k in range(self.K):
                gamma[:, k] = self.pi[k] * multivariate_normal.pdf(
                    X, mean=self.mu[k], cov=self.sigma[k]
                )
            gamma /= gamma.sum(axis=1, keepdims=True)  # 归一化

            # === M步：更新参数 ===
            N_k = gamma.sum(axis=0)  # 每个分量的有效样本数

            for k in range(self.K):
                # 更新均值
                self.mu[k] = (gamma[:, k:k+1].T @ X / N_k[k]).flatten()
                # 更新协方差
                diff = X - self.mu[k]
                self.sigma[k] = (diff.T @ (diff * gamma[:, k:k+1])) / N_k[k]
                self.sigma[k] += 1e-6 * np.eye(d)  # 正则化防奇异
            # 更新混合权重
            self.pi = N_k / n

            # 计算对数似然
            ll = 0
            for k in range(self.K):
                ll += self.pi[k] * multivariate_normal.pdf(
                    X, mean=self.mu[k], cov=self.sigma[k]
                )
            log_likelihood = np.log(ll + 1e-300).sum()  # 链式调用，连续执行多个方法
            log_likelihoods.append(log_likelihood)

            if iteration > 0 and abs(log_likelihoods[-1] - log_likelihoods[-2]) < self.tol:  # [-1]负索引取最后一个元素
                print(f"EM在第{iteration}轮收敛")
                break

        return gamma, log_likelihoods

# 生成测试数据
np.random.seed(0)
X = np.vstack([
    np.random.randn(100, 2) + [2, 2],
    np.random.randn(100, 2) + [-2, -2],
    np.random.randn(100, 2) + [2, -2]
])

gmm = GaussianMixtureEM(n_components=3)
gamma, lls = gmm.fit(X)
print(f"混合权重: {gmm.pi.round(3)}")
print(f"最终对数似然: {lls[-1]:.2f}")

import torch
import torch.nn as nn

class VAE(nn.Module):  # 继承nn.Module定义神经网络层
    def __init__(self, input_dim=784, latent_dim=20):
        super().__init__()  # super()调用父类方法
        # 编码器 q(z|x)
        self.encoder = nn.Sequential(
            nn.Linear(input_dim, 256), nn.ReLU(),
            nn.Linear(256, 128), nn.ReLU()
        )
        self.fc_mu = nn.Linear(128, latent_dim)
        self.fc_logvar = nn.Linear(128, latent_dim)

        # 解码器 p(x|z)
        self.decoder = nn.Sequential(
            nn.Linear(latent_dim, 128), nn.ReLU(),
            nn.Linear(128, 256), nn.ReLU(),
            nn.Linear(256, input_dim), nn.Sigmoid()
        )

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

    def forward(self, x):
        h = self.encoder(x)
        mu, logvar = self.fc_mu(h), self.fc_logvar(h)
        z = self.reparameterize(mu, logvar)
        return self.decoder(z), mu, logvar

def vae_loss(recon_x, x, mu, logvar):
    """ELBO = 重构损失 + KL散度"""
    recon = nn.functional.binary_cross_entropy(recon_x, x, reduction='sum')
    # KL(N(mu, sigma) || N(0, I)) 的解析解
    kl = -0.5 * torch.sum(1 + logvar - mu.pow(2) - logvar.exp())
    return recon + kl

vae = VAE()
x = torch.randn(32, 784).sigmoid()  # 模拟输入
recon, mu, logvar = vae(x)
loss = vae_loss(recon, x, mu, logvar)
print(f"ELBO loss: {loss.item():.2f}")  # .item()将单元素张量转为Python数值

import numpy as np

def metropolis_hastings(log_target, initial, n_samples=10000, proposal_std=1.0):
    """
    Metropolis-Hastings MCMC采样

    Args:
        log_target: 目标分布的对数密度（不需要归一化常数）
        initial: 初始点
        n_samples: 采样数量
        proposal_std: 提议分布标准差
    """
    samples = [initial]
    current = initial
    accepted = 0

    for _ in range(n_samples):
        # 提议新样本（对称随机游走）
        proposal = current + np.random.randn(*current.shape) * proposal_std

        # 计算接受概率（对数空间避免溢出）
        log_alpha = log_target(proposal) - log_target(current)

        # 接受/拒绝
        if np.log(np.random.rand()) < log_alpha:
            current = proposal
            accepted += 1

        samples.append(current.copy())

    acceptance_rate = accepted / n_samples
    return np.array(samples), acceptance_rate

# 示例：采样二维高斯后验
def log_posterior(theta):
    """模拟一个多模态后验分布"""
    mu1, mu2 = np.array([2, 2]), np.array([-2, -2])
    return np.logaddexp(
        -0.5 * np.sum((theta - mu1)**2),
        -0.5 * np.sum((theta - mu2)**2)
    )

samples, acc_rate = metropolis_hastings(
    log_posterior, np.array([0.0, 0.0]), n_samples=50000, proposal_std=1.5
)
print(f"接受率: {acc_rate:.2%} (经验范围: 20-50%)")
print(f"后验均值: {samples[5000:].mean(axis=0).round(3)}")  # burn-in 5000
print(f"后验标准差: {samples[5000:].std(axis=0).round(3)}")