14 - 贝叶斯方法¶

import numpy as np
from scipy import stats
import matplotlib.pyplot as plt

# 先验: Beta(α, β)
alpha_prior, beta_prior = 2, 2

# 观测数据: 10次抛掷，7次正面
n_trials = 10
n_heads = 7

# 后验: Beta(α + k, β + n - k)
alpha_post = alpha_prior + n_heads
beta_post = beta_prior + n_trials - n_heads

# 可视化
x = np.linspace(0, 1, 100)
plt.plot(x, stats.beta.pdf(x, alpha_prior, beta_prior), label='Prior')
plt.plot(x, stats.beta.pdf(x, alpha_post, beta_post), label='Posterior')
plt.legend()
plt.show()

import numpy as np

class BayesianLinearRegression:
    def __init__(self, alpha=1.0, beta=1.0):  # __init__构造方法，创建对象时自动调用
        """
        alpha: 先验精度 (1/先验方差)
        beta: 噪声精度 (1/噪声方差)
        """
        self.alpha = alpha
        self.beta = beta
        self.mean = None
        self.cov = None

    def fit(self, X, y):
        """计算后验分布"""
        n_features = X.shape[1]

        # 后验协方差
        self.cov = np.linalg.inv(
            self.alpha * np.eye(n_features) +
            self.beta * X.T @ X
        )

        # 后验均值
        self.mean = self.beta * self.cov @ X.T @ y

        return self

    def predict(self, X, return_std=True):
        """预测分布"""
        y_mean = X @ self.mean

        if return_std:
            # 预测方差
            y_var = 1/self.beta + np.sum(X @ self.cov * X, axis=1)
            y_std = np.sqrt(y_var)
            return y_mean, y_std

        return y_mean

# 使用示例
from sklearn.datasets import make_regression
from sklearn.preprocessing import StandardScaler

X, y = make_regression(n_samples=100, n_features=2, noise=10)
X = StandardScaler().fit_transform(X)

model = BayesianLinearRegression(alpha=1.0, beta=0.1)
model.fit(X, y)

# 预测带不确定性
y_pred, y_std = model.predict(X)

import numpy as np

def rbf_kernel(x1, x2, length_scale=1.0, sigma_f=1.0):
    """RBF/高斯核"""
    sqdist = (np.sum(x1**2, 1).reshape(-1, 1) +  # reshape重塑数组形状
              np.sum(x2**2, 1) - 2 * np.dot(x1, x2.T))  # np.dot计算点积/矩阵乘法
    return sigma_f**2 * np.exp(-0.5 * sqdist / length_scale**2)

def linear_kernel(x1, x2, sigma_b=0.0, sigma_v=1.0):
    """线性核"""
    return sigma_b**2 + sigma_v**2 * (x1 @ x2.T)

def periodic_kernel(x1, x2, period=1.0, length_scale=1.0):
    """周期核"""
    dist = np.abs(x1[:, None] - x2[None, :])
    return np.exp(-2 * (np.sin(np.pi * dist / period) / length_scale)**2)

class GaussianProcess:
    def __init__(self, kernel=rbf_kernel, noise=1e-5):
        self.kernel = kernel
        self.noise = noise
        self.X_train = None
        self.y_train = None
        self.K = None

    def fit(self, X, y):
        self.X_train = X
        self.y_train = y

        # 计算训练集核矩阵
        self.K = self.kernel(X, X) + self.noise * np.eye(len(X))
        return self

    def predict(self, X_test, return_std=True):
        """预测后验分布"""
        # 计算测试集与训练集的核
        K_s = self.kernel(self.X_train, X_test)
        K_ss = self.kernel(X_test, X_test)

        # 计算后验均值和协方差
        K_inv = np.linalg.inv(self.K)
        mu_s = K_s.T @ K_inv @ self.y_train

        if return_std:
            cov_s = K_ss - K_s.T @ K_inv @ K_s
            std_s = np.sqrt(np.diag(cov_s))
            return mu_s, std_s

        return mu_s

# 使用示例
import matplotlib.pyplot as plt

# 训练数据
X_train = np.array([[1], [2], [3], [4], [5]])
y_train = np.sin(X_train).ravel()

# 测试点
X_test = np.linspace(0, 6, 100).reshape(-1, 1)

# 拟合GP
gp = GaussianProcess()
gp.fit(X_train, y_train)
mu, std = gp.predict(X_test)

# 可视化
plt.plot(X_test, mu, 'b-', label='Prediction')
plt.fill_between(X_test.ravel(), mu - 2*std, mu + 2*std, alpha=0.2)
plt.scatter(X_train, y_train, c='r', label='Training data')
plt.legend()
plt.show()

from skopt import gp_minimize
from skopt.space import Real, Integer
from sklearn.ensemble import RandomForestClassifier
from sklearn.model_selection import cross_val_score

# 定义搜索空间
space = [
    Integer(10, 200, name='n_estimators'),
    Integer(2, 20, name='max_depth'),
    Real(0.01, 1.0, name='min_samples_split')
]

# 目标函数
def objective(params):
    n_estimators, max_depth, min_samples_split = params
    model = RandomForestClassifier(
        n_estimators=n_estimators,
        max_depth=max_depth,
        min_samples_split=min_samples_split,
        random_state=42
    )
    score = cross_val_score(model, X, y, cv=5).mean()
    return -score  # 最小化负准确率

# 贝叶斯优化
result = gp_minimize(
    objective,
    space,
    n_calls=50,
    random_state=42
)

print(f"最佳参数: {result.x}")
print(f"最佳得分: {-result.fun:.4f}")