贝叶斯方法（Bayesian Methods）详解

1

后验信念 = 新证据（似然）更新 × 原始信念（先验）

1
2
3

给定 p(Y|θ), p(θ)
  计算 p(θ|Y) = 贝叶斯定理
  后验有闭式，可直接求积分

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13

需要数值方法：
  1. MCMC（Markov Chain Monte Carlo）
     - Gibbs抽样
     - Metropolis-Hastings
     - 哈密尔顿蒙特卡洛（HMC）

  2. 变分推断（Variational Inference）
     - 用简单分布近似后验
     - 比MCMC快，但近似

  3. 拉普拉斯近似（Laplace Approximation）
     - 用高斯近似后验
     - 快速但简陋

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17

初始化：θ₀ = 某个初值

for iter = 1 to n_iter:

  # 1. 提议新值
  θ_new ~ q(θ_new | θ_old)  # q是提议分布，如N(θ_old, σ²I)

  # 2. 计算接受概率
  α = min(1, [p(Y|θ_new)p(θ_new)] / [p(Y|θ_old)p(θ_old)] × q(θ_old|θ_new)/q(θ_new|θ_old))

  # 3. 随机接受或拒绝
  if rand() < α:
    θ_old = θ_new  # 接受
  # else: θ_old = θ_old  # 拒绝，保持原值

  # 4. 存储样本
  samples[iter] = θ_old

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30

初始化参数 θ₀

for iter = 1 to n_mcmc:

  # 1. 提议新参数（例如随机游走）
  θ_new = θ_old + ε, 其中 ε ~ N(0, Σ_proposal)

  # 2. 计算后验比（接受概率）

  ## 2.1 计算CCP (通过求解贝尔曼方程)
  P_old = solve_bellman(θ_old)
  P_new = solve_bellman(θ_new)

  ## 2.2 计算似然
  lik_old = compute_likelihood(Y, P_old)
  lik_new = compute_likelihood(Y, P_new)

  ## 2.3 计算先验比
  prior_old = prior_density(θ_old)
  prior_new = prior_density(θ_new)

  ## 2.4 计算接受概率
  α = min(1, (lik_new × prior_new) / (lik_old × prior_old))

  # 3. 接受/拒绝
  if rand() < α:
    θ_old = θ_new

  # 4. 存储
  samples[iter] = θ_old

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21

# 阶段1：非参数估计经验CCP
P̂(x) = 从数据非参数估计（核方法、局部多项式等）

# 阶段2：贝叶斯反演参数
初始化 θ₀

for iter = 1 to n_mcmc:

  θ_new = θ_old + ε

  # 使用固定的经验CCP P̂，而非求解DP
  # （P̂不随θ变化）

  lik_new = compute_likelihood(Y, P̂)  # P̂固定，快速

  α = min(1, (lik_new × prior_new) / (lik_old × prior_old))

  if rand() < α:
    θ_old = θ_new

  samples[iter] = θ_old

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18

for iter = 1 to n_mcmc:

  θ_new = θ_old + ε

  # 模拟虚拟数据并计算矩
  m_new = simulate_and_compute_moments(θ_new, N_sim=5000)
  m_old = simulate_and_compute_moments(θ_old, N_sim=5000)

  # 计算伪似然
  pseudo_lik_new = exp(-0.5 * (m_new - m̂)' W (m_new - m̂))
  pseudo_lik_old = exp(-0.5 * (m_old - m̂)' W (m_old - m̂))

  α = min(1, (pseudo_lik_new × prior_new) / (pseudo_lik_old × prior_old))

  if rand() < α:
    θ_old = θ_new

  samples[iter] = θ_old

1
2
3

选3-5个先验（从保守到信息性强）
用各先验做贝叶斯估计
比较后验

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11

初始化 Σ₀

for iter = 1 to n:

  # M-H步骤，用当前的Σ
  ...

  # 每K步，更新提议方差
  if mod(iter, K) == 0:
    Σ = empirical_cov(samples[1:iter])
    Σ = 0.234/0.234 × 2.38²/dim × Σ  # 自适应缩放

1
2
3
4
5
6
7

θ
│     ╱╲      ╱╲╱╲
│    ╱  ╲____╱      ╲____
│___╱
│
└─────────────────── iter
  burn-in    后验样本

1
2
3
4
5
6
7

n_mcmc = 10000
burn_in = 2000
thin = 10

# 使用第2000-10000的每10个样本
post_samples = samples[2001:10000:10]
n_post = len(post_samples) = 800

1

theta_mean = mean(post_samples)

1

theta_std = std(post_samples)

1
2

ci_lower = quantile(post_samples, 0.025)
ci_upper = quantile(post_samples, 0.975)

1
2
3
4
5
6

import matplotlib.pyplot as plt

fig, axes = plt.subplots(1, k)
for j in range(k):
  axes[j].hist(post_samples[:, j], bins=50, density=True)
  axes[j].set_title(f'θ_{j}')

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
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

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

class BayesianDynamicModel:

  def __init__(self, data, prior):
    self.data = data
    self.prior = prior
    self.N = len(data)

  # ===== 模型 =====
  def compute_CCP(self, theta):
    """
    从参数计算CCP（通过求解DP或其他方式）
    简化版本：假设已有某个映射
    """
    theta_rc, theta_cost, beta = theta
    # ... 求解贝尔曼方程 ...
    return P  # CCP

  def log_likelihood(self, theta):
    """计算对数似然"""
    P = self.compute_CCP(theta)

    # 对数似然 = Σ d*log(P) + (1-d)*log(1-P)
    log_lik = 0
    for i in range(self.N):
      d = self.data[i]['action']
      log_lik += d * np.log(P[i] + 1e-10) + (1-d) * np.log(1 - P[i] + 1e-10)

    return log_lik

  def log_prior(self, theta):
    """计算对数先验"""
    theta_rc, theta_cost, beta = theta

    # 先验：RC ~ Gamma, cost ~ Normal, beta ~ Beta
    lp = self.prior['RC'].logpdf(theta_rc) + \
         self.prior['cost'].logpdf(theta_cost) + \
         self.prior['beta'].logpdf(beta)

    return lp

  def log_posterior(self, theta):
    """计算对数后验 ∝ 似然 + 先验"""
    return self.log_likelihood(theta) + self.log_prior(theta)

  # ===== MCMC: Metropolis-Hastings =====
  def mcmc(self, theta_init, n_iter=10000, proposal_scale=0.1):
    """
    Metropolis-Hastings算法
    """
    dim = len(theta_init)
    samples = np.zeros((n_iter, dim))
    accepts = 0

    theta_current = theta_init
    log_post_current = self.log_posterior(theta_current)

    for iter in range(n_iter):

      # 1. 随机游走提议
      epsilon = np.random.normal(0, proposal_scale, dim)
      theta_new = theta_current + epsilon

      # 2. 边界检验（某些参数有约束）
      if not self.in_bounds(theta_new):
        samples[iter] = theta_current
        continue

      # 3. 计算接受概率
      log_post_new = self.log_posterior(theta_new)
      log_alpha = log_post_new - log_current

      # 4. 接受/拒绝
      if np.log(np.random.uniform()) < log_alpha:
        theta_current = theta_new
        log_post_current = log_post_new
        accepts += 1

      # 5. 存储样本
      samples[iter] = theta_current

      # 进度报告
      if (iter + 1) % 1000 == 0:
        accept_rate = accepts / (iter + 1)
        print(f"Iter {iter+1}: accept_rate = {accept_rate:.2%}")

    self.accept_rate = accepts / n_iter
    return samples

  def in_bounds(self, theta):
    """检查参数是否在有效范围"""
    theta_rc, theta_cost, beta = theta
    return theta_rc > 0 and 0 < beta < 1

  # ===== 诊断 =====
  def plot_trace(self, samples, burn_in=2000):
    """绘制trace图"""
    dim = samples.shape[1]
    fig, axes = plt.subplots(dim, 1, figsize=(10, 3*dim))

    for j in range(dim):
      axes[j].plot(samples[:, j], alpha=0.7)
      axes[j].axvline(burn_in, color='r', linestyle='--', label='burn-in')
      axes[j].set_ylabel(f'θ_{j}')

    plt.tight_layout()
    plt.show()

  def plot_posterior(self, samples, burn_in=2000):
    """绘制后验分布"""
    post_samples = samples[burn_in:, :]
    dim = post_samples.shape[1]

    fig, axes = plt.subplots(dim, 1, figsize=(10, 3*dim))

    for j in range(dim):
      axes[j].hist(post_samples[:, j], bins=50, density=True)
      axes[j].set_ylabel(f'p(θ_{j} | data)')

    plt.tight_layout()
    plt.show()

  def posterior_summary(self, samples, burn_in=2000):
    """计算后验统计"""
    post_samples = samples[burn_in:, :]

    mean = np.mean(post_samples, axis=0)
    std = np.std(post_samples, axis=0)
    ci_lower = np.quantile(post_samples, 0.025, axis=0)
    ci_upper = np.quantile(post_samples, 0.975, axis=0)

    print("="*60)
    print("POSTERIOR SUMMARY")
    print("="*60)
    print(f"{'Parameter':<15} {'Mean':<12} {'Std':<12} {'95% CI':<25}")
    print("-"*60)

    param_names = ['theta_RC', 'theta_cost', 'beta']
    for i, name in enumerate(param_names):
      ci_str = f"[{ci_lower[i]:.4f}, {ci_upper[i]:.4f}]"
      print(f"{name:<15} {mean[i]:<12.4f} {std[i]:<12.4f} {ci_str:<25}")

    return {
      'mean': mean,
      'std': std,
      'ci_lower': ci_lower,
      'ci_upper': ci_upper
    }

# ===== 使用示例 =====

# 定义先验
prior = {
  'RC': gamma(a=4, scale=5),      # 均值20
  'cost': norm(loc=0, scale=0.5),  # 均值0，标准差0.5
  'beta': beta(a=10, b=2)          # 均值0.833
}

# 创建模型
model = BayesianDynamicModel(data, prior)

# 初值
theta_init = np.array([15, -0.1, 0.8])

# MCMC
samples = model.mcmc(theta_init, n_iter=10000, proposal_scale=0.05)

# 诊断
model.plot_trace(samples, burn_in=2000)
model.plot_posterior(samples, burn_in=2000)
summary = model.posterior_summary(samples, burn_in=2000)

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55

def adaptive_mcmc(self, theta_init, n_iter=10000,
                  adapt_interval=500, target_accept=0.234):
  """
  自适应Metropolis-Hastings
  根据接受率自动调整提议的尺度
  """
  dim = len(theta_init)
  samples = np.zeros((n_iter, dim))

  theta_current = theta_init
  log_post_current = self.log_posterior(theta_current)

  # 初始提议方差
  proposal_var = np.ones(dim) * 0.1

  accepts = 0

  for iter in range(n_iter):

    # 提议：独立维度
    epsilon = np.random.normal(0, np.sqrt(proposal_var))
    theta_new = theta_current + epsilon

    if not self.in_bounds(theta_new):
      samples[iter] = theta_current
      continue

    # 接受/拒绝
    log_post_new = self.log_posterior(theta_new)
    log_alpha = log_post_new - log_post_current

    if np.log(np.random.uniform()) < log_alpha:
      theta_current = theta_new
      log_post_current = log_post_new
      accepts += 1

    samples[iter] = theta_current

    # 每adapt_interval步，调整提议方差
    if (iter + 1) % adapt_interval == 0:
      accept_rate = accepts / (iter + 1)

      # 自适应缩放：如果接受率太高，扩大方差；太低，缩小方差
      scale_factor = 1.0
      if accept_rate > target_accept:
        scale_factor = 1.1  # 扩大
      else:
        scale_factor = 0.9  # 缩小

      proposal_var *= scale_factor**2

      print(f"Iter {iter+1}: accept_rate = {accept_rate:.2%}, " +
            f"proposal_scale = {np.sqrt(proposal_var[0]):.4f}")

  return samples

1

用CCP方法估计 θ̂_freq

1
2
3

先验：p(θ) = Normal(θ̂_freq, Σ_freq)

然后用MCMC得到完整的后验

1
2
3

for s = 1 to n_post:
  θ^(s) ~ 后验
  Y^rep,s ~ p(Y | θ^(s))

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29

def posterior_predictive_check(self, samples, test_stat_func, burn_in=2000):
  """
  后验预测检验
  """
  post_samples = samples[burn_in:, :]

  # 计算真实数据的统计量
  T_real = test_stat_func(self.data)

  # 计算后验预测的统计量
  T_rep = []

  for theta in post_samples:
    Y_rep = self.simulate_data(theta)
    T_rep.append(test_stat_func(Y_rep))

  T_rep = np.array(T_rep)

  # 绘图
  plt.hist(T_rep, bins=50, density=True, alpha=0.7, label='Posterior Predictive')
  plt.axvline(T_real, color='r', linewidth=2, label='Observed')
  plt.legend()
  plt.show()

  # p值（后验预测与观测一致的概率）
  p_value = np.mean(T_rep >= T_real)
  print(f"Posterior Predictive p-value: {p_value:.3f}")

  return T_rep, p_value

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14

def sensitivity_analysis(self, scenarios):
  """
  运行多个情景的贝叶斯估计
  scenarios: [{'name': '...', 'prior': ..., 'data': ...}, ...]
  """
  results = {}

  for scenario in scenarios:
    model = BayesianDynamicModel(scenario['data'], scenario['prior'])
    samples = model.mcmc(theta_init, n_iter=5000)
    summary = model.posterior_summary(samples, burn_in=2000)
    results[scenario['name']] = summary

  return results

1
2

RC ~ Gamma(4, 0.2)        # 均值20，标准差10
β ~ Beta(10, 2)           # 均值0.833，标准差0.085

1
2
3
4
5
6
7

Accept rate: 25.3%
Effective sample size: 3847 / 5000

Parameter        Mean      Std      95% CI
─────────────────────────────────────────
RC              19.2      3.1    [13.5, 25.8]
beta            0.829     0.018   [0.794, 0.864]

1
2
3

1. 非参数估计经验CCP P̂(x₁, x₂)
2. 定义伪似然：p(Y|θ) ∝ Π P̂(d_t | x_t)^{d_t} [1-P̂(d_t|x_t)]^{1-d_t}
3. 运行MCMC估计参数后验

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14

START

有强的先验信息吗？
├─ YES, 有: 贝叶斯方法
│  ├─ 计算资源充足 → 贝叶斯NFXP
│  └─ 计算资源有限 → 贝叶斯CCP/SMM
│
└─ NO, 没有: 频率派方法
   ├─ 样本很大 → NFXP（最有效）
   ├─ 样本中等 → CCP或GMM
   ├─ 样本小或模型复杂 → SMM或间接推断
   └─ 维度高(≥3) → GMM/SMM/间接推断

END

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32

┌─────────────────────────────────────┐
│   贝叶斯方法的快速检查清单           │
├─────────────────────────────────────┤
│                                     │
│ 【前期准备】                        │
│ ☐ 先验已指定？                      │
│ ☐ 先验的合理性已检验？              │
│ ☐ 似然函数正确实现？                │
│                                     │
│ 【MCMC诊断】                        │
│ ☐ 接受率在20-50%？                  │
│ ☐ Rhat < 1.1？                      │
│ ☐ ESS充分？                         │
│ ☐ Trace图稳定？                     │
│                                     │
│ 【后验推断】                        │
│ ☐ 后验统计计算无误？                │
│ ☐ 可信区间合理？                    │
│ ☐ 后验预测检验完成？                │
│                                     │
│ 【敏感性与报告】                    │
│ ☐ 多先验敏感性检验？                │
│ ☐ 模型比较（若需）？                │
│ ☐ 完整报告先验与诊断？              │
│                                     │
└─────────────────────────────────────┘

【贝叶斯方法的三行总结】

1. 定义：参数是随机变量，目标是后验分布
2. 方法：Bayes定理 + MCMC计算
3. 优势：完整分布，自然模型比较，稳健小样本