diff --git a/Figure_2.png b/Figure_2.png new file mode 100644 index 0000000..b4a7380 Binary files /dev/null and b/Figure_2.png differ diff --git a/README.md b/README.md new file mode 100644 index 0000000..ef0d906 --- /dev/null +++ b/README.md @@ -0,0 +1,99 @@ +# MOTIVATION + +论文中的辨识是认为系统的阶数已知的辨识，而且 MCMC 的抽样的变量是特征多项式矩阵，我想使用特征值来抽样，并同时考虑系统阶数未知的情况。 + +采用RJMMC的方法 + +## 细致平稳条件 + +$$ +\alpha \left( x,y \right) =\min \left\{ 1,R \right\} +\\ +\pi \left( x \right) p\left( x,y \right) =\pi \left( y \right) p\left( y,x \right) +$$ + +## 跨维度游动 + +### k <-> k+1 + +1. Birth move + $$ + \log R=\log \frac{P\left( x_{k+1} \right) P\left( y|x_{k+1} \right)}{P\left( x_k \right) P\left( y|x_k \right)}+\log \frac{P\left( dead \right)}{P\left( birth \right)}+\log \frac{\small{\frac{1}{k+1}}}{q\left( u \right)}+\log \left| J_1 \right| + $$ +2. Death move + $$ + \log R=\log \frac{P\left( x_k \right) P\left( y|x_k \right)}{P\left( x_{k+1} \right) P\left( y|x_{k+1} \right)}+\log \frac{P\left( birth \right)}{P\left( dead \right)}+\log \frac{\small{q\left( u \right)}}{\frac{1}{k+1}}-\log \left| J_1 \right| + $$ + +### k <-> k+2 + +1. Birth move + $$ + \log R=\log \frac{P\left( x_{k+2} \right) P\left( y|x_{k+2} \right)}{P\left( x_k \right) P\left( y|x_k \right)}+\log \frac{P\left( dead \right)}{P\left( birth \right)}+\log \frac{\small{\frac{1}{m+1}}}{q\left( \theta \right) q\left( r \right)}+\log \left| J_2 \right| + $$ +2. Death move + $$ + \log R=\log \frac{P\left( x_k \right) P\left( y|x_k \right)}{P\left( x_{k+2} \right) P\left( y|x_{k+2} \right)}+\log \frac{P\left( birth \right)}{P\left( dead \right)}+\log \frac{\small{q\left( \theta \right) q\left( r \right)}}{\frac{1}{m+1}}-\log \left| J_2 \right| + $$ + +$其中 m$ 是当前系统中复共轭极点对的数量 +$J_1$ 和 $J_2$ 是雅可比矩阵行列式 +$\left| J_1 \right|=\left| \prod_{\boldsymbol{i}=1}^{\boldsymbol{k}}{\left( \lambda _i-u_i \right)} \right|$ +$\left| J_2 \right|=\left| \prod_{\boldsymbol{i}=1}^{\boldsymbol{k}}{\left( \lambda _{k+1}-\lambda _i \right) \left( \lambda _{k+2}-\lambda _i \right)} \right|\left| \lambda _{k+1}-\lambda _{k+1} \right|$ + +## 同维度游动 + +### 根的类型转换 + +先确定映射关系，令 + +$$ +a=\frac{r_1+r_2}{2}\text{、}b=\frac{r_1-r_2}{2} +\\ +r_1=a+b\text{、}r_2=a-b +$$ + + +再确定雅可比矩阵行列式 +$$ +J_{C\rightarrow R}=\left| \det \left( \begin{matrix} + \frac{\partial a}{\partial r_1}& \frac{\partial a}{\partial r_2}\\ + \frac{\partial b}{\partial r_1}& \frac{\partial b}{\partial r_2}\\ +\end{matrix} \right) \right|=\frac{1}{2} +\\ +J_{R\rightarrow C}={J_{C\rightarrow R}}^{-1}=2 +$$ + +1. 实根 <-> 复共轭根对 + 在$n_r$个实根中任选两个实根$r_1$和$r_2$，转换为复共轭根对$a\pm jb$，其中$a=\frac{r_1+r_2}{2}$，$b=\frac{r_1-r_2}{2}$。 + +$$ +q\left( x|x' \right) =p_m\times \frac{1}{C\left( n_r,2 \right)} +$$ + +此时的接受率写作:$Merge\text{：}\alpha \left( x',x \right) =\min \left\{ 1,\frac{\pi \left( x \right)}{\pi \left( x' \right)}\cdot \frac{1}{q\left( x|x' \right)}\cdot \frac{1}{2} \right\} $ + +2. 复共轭根对 <-> 实根 + 在$n_c$个复共轭根对中任选一个复共轭根对$a\pm jb$，转换为两个实根$r_1=a+b$和$r_2=a-b$。 + +$$ +q\left( x'|x \right) =p_c\times \frac{1}{n_c} +$$ + +此时的接受率写作:$Split\text{：}\alpha \left( x,x' \right) =\min \left\{ 1,\frac{\pi \left( x' \right)}{\pi \left( x \right)}\cdot \frac{1}{q\left( x'|x \right)}\cdot 2 \right\} $ + +### 同类型根的调整 + +1. 实根调整 + 选择一个实根$r$，通过添加噪声$\epsilon \sim N\left( 0,\sigma ^2 \right)$调整该实根的位置为$r' = r + \epsilon$。 + +$$ +\alpha \left( x,x' \right) =\min \left\{ 1,\frac{\pi \left( x' \right)}{\pi \left( x \right)} \right\} +$$ + +2. 复共轭根对调整 + 选择一个复共轭根对$a\pm jb$，通过添加噪声$\epsilon _a \sim N\left( 0,\sigma _a^2 \right)$和$\epsilon _b \sim N\left( 0,\sigma _b^2 \right)$调整该复共轭根对的位置为$a' = a + \epsilon _a$和$b' = b + \epsilon _b$。 + +$$ +\alpha \left( x,x' \right) =\min \left\{ 1,\frac{\pi \left( x' \right)}{\pi \left( x \right)} \right\} +$$ diff --git a/RJMCMC_ModelSelector.py b/RJMCMC_ModelSelector.py new file mode 100644 index 0000000..03ccd45 --- /dev/null +++ b/RJMCMC_ModelSelector.py @@ -0,0 +1,648 @@ +import numpy as np +from scipy import stats + + +class RJMCMC_Sampler: + def __init__(self, y, u, k_min=1, k_max=10, initial_k=None, + initial_a=None, initial_b=None): + """ + 初始化采样器 + + 参数: + y (array): 观测到的输出序列 + u (array): 已知的输入序列 + k_min (int): 探索的最小模型阶数 + k_max (int): 探索的最大模型阶数 + initial_k (int, optional): 初始模型阶数 + initial_a (dict, optional): 初始多项式系数 + initial_b (array, optional): 初始多项式系数 + """ + + # 定义观测数据 + self.y = y + self.u = u + self.N = len(y) + + # 定义模型阶数 + self.k_min = k_min + self.k_max = k_max + gamma_sample = np.random.gamma(shape=2, scale=1) + self.current_k = initial_k if initial_k is not None else min(max(int(gamma_sample), k_min), k_max) + + # 初始化多项式系数 + if initial_a is not None: + self.current_a = initial_a + else: + self.current_a = self._generate_stable_a(self.current_k) + + # 定义标准型的b系数 + if initial_b is not None: + self.current_b = initial_b + else: + self.current_b = self._generate_observable_b(self.current_k) + + # 初始化噪声先验 (假设为 Gamma 先验) + self.lambda_a_w = 1e-3 # 过程噪声精度的先验 + self.lambda_b_w = 1e-3 + self.lambda_a_z = 1e-3 # 观测噪声精度的先验 + self.lambda_b_z = 1e-3 + + # 初始化噪声方差 (从先验采样) + self.current_sigma2_w = 1.0 / np.random.gamma(self.lambda_a_w, 1.0/self.lambda_b_w) + self.current_sigma2_z = 1.0 / np.random.gamma(self.lambda_a_z, 1.0/self.lambda_b_z) + + + def _update_current_eigenvalues(self): + """一个辅助函数，根据 current_a 更新 current_eigenvalues。""" + if self.current_k > 0: + poly_coeffs = np.concatenate(([1], -self.current_a[::-1])) + self.current_eigenvalues = np.roots(poly_coeffs) + else: + self.current_eigenvalues = np.array([]) + + def _cal_current_eigenvalues(self, a_coeffs_temp): + """一个辅助函数，根据给定的 a_coeffs 计算对应的特征值。""" + k = len(a_coeffs_temp) + if k > 0: + poly_coeffs = np.concatenate(([1], -a_coeffs_temp[::-1])) + return np.roots(poly_coeffs) + else: + return np.array([]) + + def _cal_current_a_coeffs(self, eigenvalues_temp): + """一个辅助函数，根据给定的特征值计算对应的 a_coeffs。""" + k = len(eigenvalues_temp) + if k > 0: + # np.poly 返回 [1, a_{k-1}, ..., a_0]，去掉首项“1”，反转剩下的 + poly_coeffs = np.poly(eigenvalues_temp) + return np.real(poly_coeffs[1:][::-1]) + else: + return np.array([]) + + + def _log_prior_eigenvalues(self, eigenvalues): + """ + 计算 k 个特征值 (Lambda) 的对数先验概率 log p(Lambda)。 + + 该先验基于论文 4.1.1 节 讨论的先验： + 1. 稳定性：如果任何 |lambda| >= 1，先验概率为 0 (log(0) = -inf)。 + 2. 实数根 (p_real)：在 (-1, 1) 上均匀分布, p(lambda_r) = 1/2。 + log p(lambda_r) = -log(2)。 + 3. 复数根 (p_complex)：在单位圆盘上均匀分布, p(lambda_c) = 1/pi。 + log p(lambda_c) = -log(pi)。 + + 注意：复共轭对 (lambda, lambda_conj) 只计算一次 (例如上半平面)。 + 我们通过只计算 imag(lambda) >= 0 的复数根来实现。 + """ + + log_prior = 0.0 + + # 跟踪已处理的复数根 (避免重复计算共轭对) + processed_complex = set() + + for eig in eigenvalues: + if np.abs(eig) >= 1.0: + # 违反稳定性约束 + return -np.inf + + if np.isclose(eig.imag, 0): + # 是实数 + # 这是 "Uniform real-eigenvalue prior" + log_prior += -np.log(2.0) + + else: + # 是复数 + # 检查是否已处理过 + if eig in processed_complex or np.conjugate(eig) in processed_complex: + continue + + # 我们只计算 imag >= 0 的根 (代表一对) + if eig.imag >= 0: + # 这是 "Polar-coordinate prior" + log_prior += -np.log(np.pi) + processed_complex.add(eig) + + return log_prior + + def _log_prior_eigenvalues_to_a(self,eigenvalues): + """ + 计算从特征值 (eigenvalues) 到多项式系数 a = [a_0, ..., a_{k-1}] 的对数先验概率。 + + 该方法基于 Proposition 4.2 (Change of variables) + log p(a) = log p(Lambda) - log(|det J|) + + 其中 log p(Lambda) 是特征值的先验 + 而 log(|det J|) 是维塔变换的对数雅可比行列式。 + """ + # 计算 log p(Lambda) + log_prior_lambda = self._log_prior_eigenvalues(eigenvalues) + + # 如果特征值不稳定 (log p(Lambda) = -inf)，则 p(a) 也不可能 + if log_prior_lambda == -np.inf: + return -np.inf + + # 计算对数雅可比行列式 log(|det J|) + log_det_jacobian = 0.0 + k = len(eigenvalues) + if k < 2: + return log_prior_lambda + + for i in range(k): + + for j in range(i + 1, k): + diff = eigenvalues[i] - eigenvalues[j] + abs_diff = np.abs(diff) + if abs_diff == 0: + return -np.inf # 重根导致雅可比行列式为零 + log_det_jacobian += np.log(abs_diff) + + if log_det_jacobian == -np.inf: + return -np.inf + + # 应用变量替换公式 + log_prior_a = log_prior_lambda - log_det_jacobian + + return log_prior_a + + def _log_likelihood(self, k, a_coeffs, b_coeffs, sigma2_w, sigma2_z): + """ + 计算给定参数下的对数似然 log p(y | k, a, b, sigma_w, sigma_z)。 + + 使用标准卡尔曼滤波器 (Kalman Filter) (如论文 Appendix B ) + 和预测误差分解 (Prediction Error Decomposition) (如论文 Eq. (D.37) )。 + """ + + # 如果 k=0，模型无法运行 + if k == 0: + # 返回一个非常小的似然 + return -np.inf + + # 获取当前参数下的状态空间矩阵（以参数为输入，避免依赖实例状态） + A, B, C, D = self.get_controller_canonical_form(k, a_coeffs, b_coeffs) + + # 过程噪声 Sigma (k x k) + Sigma = np.zeros((k, k)) + Sigma[k-1, k-1] = sigma2_w + + # 测量噪声 Gamma (1 x 1) + Gamma = np.array([[sigma2_z]]) + + # 初始化卡尔曼滤波器 + + # 状态 x_t|t-1 + x_pred = np.zeros((k, 1)) + + # 状态协方差 P_t|t-1 + P_pred = np.eye(k) * 1e6 + + total_log_likelihood = 0.0 + + # 运行卡尔曼滤波器 N 步 + for t in range(self.N): + y_t = self.y[t] + u_t = self.u[t] + # --- 预测步 (Prediction) --- + # (在 t=0 时，x_pred 和 P_pred 是我们初始化的 x_0|-1, P_0|-1) + + # --- 计算似然 --- + # 预测误差 nu_t + y_pred = C @ x_pred + D @ u_t + nu_t = y_t - y_pred + + # 预测误差协方差 S_t + S_t = C @ P_pred @ C.T + Gamma + S_t_inv = 1.0 / S_t + + # 累加对数似然 + total_log_likelihood += -0.5 * (np.log(2 * np.pi) + np.log(S_t) + (nu_t * S_t_inv * nu_t)) + + # --- 更新步 (Update) --- + # 卡尔曼增益 K_t (k x 1) + K_t = P_pred @ C.T * S_t_inv + + # 更新状态 x_t|t + x_update = x_pred + K_t @ nu_t + + # 更新协方差 P_t|t + I_KC = np.eye(k) - K_t @ C + P_update = I_KC @ P_pred @ I_KC.T + K_t @ Gamma @ K_t.T + + # --- 为下一次循环准备预测 (t+1) --- + x_pred = A @ x_update + B @ u_t + P_pred = A @ P_update @ A.T + Sigma + + return total_log_likelihood # 返回标量值 + + def _propose_birth(self): + """ + 提议一个 "诞生" 转移 (k -> k+1 或 k -> k+2)。 + 返回包含提议状态和计算接受率所需的对数比率的字典。 + """ + k = self.current_k + current_eigs = self.current_eigenvalues + + # 决定是诞生一个实数根还是一对复共轭根 + can_add_complex = (k + 2) <= self.k_max + add_real = True + if can_add_complex: + # 以0.5的概率选择诞生一对复共轭根 + if np.random.rand() < 0.5: + add_real = False + + if add_real: + k_new = k + 1 + + # 从提议分布 q(u) 中采样辅助变量 u = (u_lambda, u_b) + u_lambda = np.random.uniform(-1.0, 1.0) # 新特征值 + u_b = np.random.normal(0, 1) # 新 b 系数 + + # 计算对数提议密度 log(q(u)) + log_q_forward = -np.log(2.0) + stats.norm.logpdf(u_b, 0, 1) + + # 计算对数雅可比行列式 log|J| + if k == 0: + log_det_jacobian = 0.0 + else: + log_det_jacobian = np.sum(np.log(np.abs(current_eigs - u_lambda))) + + # 构造新状态 + new_eigs = np.append(current_eigs, u_lambda) + new_a = np.real(np.poly(new_eigs)[1:][::-1]) + new_b = np.append(self.current_b, u_b) + + # 接受率中的项为 |J| / q(u)，在对数空间中为 log|J| - log(q(u)) + log_ratio = log_det_jacobian - log_q_forward + + return {"k_new": k_new, "a_new": new_a, "b_new": new_b, "log_ratio": log_ratio, "type": "birth_real", "new_eigs": new_eigs} + + else: + # --- 诞生一对复共轭根 (k -> k+2) --- + k_new = k + 2 + + # 采样辅助变量 u = (rho, theta, u_b1, u_b2) + rho = np.sqrt(np.random.uniform(0, 1.0)) + theta = np.random.uniform(0, np.pi) + u_lambda = rho * (np.cos(theta) + 1j * np.sin(theta)) + u_b1, u_b2 = np.random.normal(0, 1, 2) + + # 计算对数提议密度 log(q(u)) + log_q_forward = -np.log(np.pi) + stats.norm.logpdf(u_b1, 0, 1) + stats.norm.logpdf(u_b2, 0, 1) + + # 计算对数雅可比行列式 log|J| + # |J| = |product(|u_lambda - lambda_i|^2) * (2*Im(u_lambda))| + if k == 0: + log_jacobian = np.log(np.abs(2 * u_lambda.imag)) + else: + log_jacobian = np.sum(np.log(np.abs(u_lambda - current_eigs)**2)) + \ + np.log(np.abs(2 * u_lambda.imag)) + + # 构造新状态 + new_eigs = np.append(current_eigs, [u_lambda, np.conjugate(u_lambda)]) + new_a = np.real(np.poly(new_eigs)[1:][::-1]) + new_b = np.append(self.current_b, [u_b1, u_b2]) + + log_ratio = log_jacobian - log_q_forward + + return {"k_new": k_new, "a_new": new_a, "b_new": new_b, "log_ratio": log_ratio, "type": "birth_real", "new_eigs": new_eigs} + + def _propose_death(self): + """ + 提议一个 "消亡" 转移 (k -> k-1 或 k -> k-2)。 + 返回包含提议状态和计算接受率所需的对数比率的字典。 + """ + k = self.current_k + current_eigs = self.current_eigenvalues + + # 识别出可移除的实数根和复共轭对 + real_eigs_indices = np.where(np.isreal(current_eigs)) + complex_eigs_indices = np.where((np.iscomplex(current_eigs)) & (current_eigs.imag > 0)) #type: ignore + + can_remove_real = len(real_eigs_indices) > 0 + can_remove_complex = len(complex_eigs_indices) > 0 + + if not can_remove_real and not can_remove_complex: + return None # 无法执行消亡 + + # 决定是移除实数根还是复共轭对 + remove_real = True + if can_remove_real and can_remove_complex: + if np.random.rand() < 0.5: + remove_real = False + elif can_remove_complex: + remove_real = False + + if remove_real: + # --- 消亡一个实数根 (k -> k-1) --- + k_new = k - 1 + + # 1. 随机选择一个实数根移除 + idx_to_remove = np.random.choice(real_eigs_indices) + lambda_removed = current_eigs[idx_to_remove] + + # 移除的根和b系数构成了逆向（诞生）提议的辅助变量 u + u_lambda = lambda_removed + u_b = self.current_b[-1] + + # 2. 计算逆向提议的对数密度 log(q(u)) + log_q_reverse = -np.log(2.0) + stats.norm.logpdf(u_b, 0, 1) + + # 3. 计算对应诞生过程的对数雅可比行列式 log|J| + remaining_eigs = np.delete(current_eigs, idx_to_remove) + if k_new == 0: + log_jacobian_birth = 0.0 + else: + log_jacobian_birth = np.sum(np.log(np.abs(u_lambda - remaining_eigs))) + + # 构造新状态 + new_a = np.real(np.poly(remaining_eigs)[1:][::-1]) + new_b = self.current_b[:-1] + + # 接受率中的项为 q_reverse(u) / |J_birth|，在对数空间中为 log(q_reverse) - log|J_birth| + log_ratio = log_q_reverse - log_jacobian_birth + + return {"k_new": k_new, "a_new": new_a, "b_new": new_b, "log_ratio": log_ratio, "type": "death_real", "new_eigs": remaining_eigs} + + else: + # --- 消亡一对复共轭根 (k -> k-2) --- + k_new = k - 2 + + # 随机选择一对共轭根移除 + complex_idx_to_remove = np.random.choice(complex_eigs_indices) + lambda_removed = current_eigs[complex_idx_to_remove] + + # 找到其共轭对 + conjugate_idx_to_remove = np.where(current_eigs == np.conjugate(lambda_removed)) + + # 逆向提议的辅助变量 u + u_lambda = lambda_removed + u_b1, u_b2 = self.current_b[-2:] + + # 计算逆向提议的对数密度 log(q(u)) + log_q_reverse = -np.log(np.pi) + stats.norm.logpdf(u_b1, 0, 1) + stats.norm.logpdf(u_b2, 0, 1) + + # 计算对应诞生过程的对数雅可比行列式 + remaining_eigs = np.delete(current_eigs, [complex_idx_to_remove, conjugate_idx_to_remove]) + if k_new == 0: + log_jacobian_birth = np.log(np.abs(2 * u_lambda.imag)) #type: ignore + else: + log_jacobian_birth = np.sum(np.log(np.abs(u_lambda - remaining_eigs)**2)) + \ + np.log(np.abs(2 * u_lambda.imag)) #type: ignore + + # 构造新状态 + new_a = np.real(np.poly(remaining_eigs)[1:][::-1]) + new_b = self.current_b[:-2] + + log_ratio = log_q_reverse - log_jacobian_birth + + return {"k_new": k_new, "a_new": new_a, "b_new": new_b, "log_ratio": log_ratio, "type": "death_complex", "new_eigs": remaining_eigs} + + def _log_prior_b(self, b_coeffs): + """ + 计算 C_c 系数 (b_0, ..., b_{k-1}) 的对数先验概率。 + 根据论文 4.2 节，我们使用独立标准正态先验。 + b_i ~ N(0, 1) + """ + # + return stats.norm.logpdf(b_coeffs, 0, 1).sum() + + def _generate_stable_a(self, k): + """ + 生成 k 个稳定的 a = [a_0, a_1, ..., a_{k-1}] 系数, 并保存对应的特征值。 + + 该方法基于论文 4.1.1 节 [cite: 314-320] 和 6.2 节 的讨论， + 生成 k 个稳定的特征值 (根)，然后使用维塔定理 (Proposition 4.1) [cite: 290] + 计算多项式系数。 + + 为确保系数a为实数，特征值必须是实数或复共轭对 。 + """ + + eigenvalues = [] + i = 0 + + while i < k: + # 如果 k-i > 1，随机决定采样一个实数根还是一个复共轭对 + # 如果 k-i == 1，则只能采样一个实数根 + if i == k - 1 or np.random.rand() < 0.5: + # --- 1. 采样一个实数根 --- + # 使用 "Uniform real-eigenvalue prior" (在 (-1, 1) 上均匀分布) + real_eig = np.random.uniform(-1.0, 1.0) + eigenvalues.append(real_eig) + i += 1 + else: + # --- 2. 采样一个复共轭对 --- + # 使用 "Polar-coordinate prior" (极坐标先验) + # 为确保在单位圆盘内面积均匀，rho = sqrt(U(0,1)) + rho = np.sqrt(np.random.uniform(0, 1.0)) + + # theta 在 (0, pi) 均匀分布 (上半平面) + theta = np.random.uniform(0, np.pi) + + # 计算复特征值 + complex_eig = rho * (np.cos(theta) + 1j * np.sin(theta)) + + # 添加该值及其共轭 + eigenvalues.append(complex_eig) + eigenvalues.append(np.conjugate(complex_eig)) + i += 2 + + # --- 3. (维塔定理) 从根计算多项式系数 --- + # np.poly(roots) 返回 [1, a_{k-1}, a_{k-2}, ..., a_0] + coefficients_descending = np.poly(eigenvalues) + + # 返回 [a_0, a_1, ..., a_{k-1}] + a_coefficients = coefficients_descending[1:][::-1] + + self.current_eigenvalues = np.array(eigenvalues) # 顺便保存特征值 + return np.real(a_coefficients) + + + def _generate_observable_b(self, k): + """ + 生成 b = [b_0, ..., b_{k-1}] 系数 (共 k 个)。 + """ + b_array = np.random.normal(0, 1, size=k) + return b_array + + + def get_controller_canonical_form(self, k=None, a_coeffs=None, b_coeffs=None): + """ + 返回给定参数的规范型 (Controller Canonical Form) 矩阵 A_c, B_c, C_c, D_c。 + + 如果提供了 k、a_coeffs、b_coeffs，则使用这些输入（使函数成为纯函数）； + 否则回退到实例的当前状态（向后兼容）。 + """ + # 回退到实例状态以保持向后兼容 + if k is None: + k = self.current_k + if a_coeffs is None: + a_coeffs = self.current_a + if b_coeffs is None: + b_coeffs = self.current_b + + # k=0 是无效情况，但 k=1 时 np.diag(np.ones(0), 1) 会创建一个空 1x1 矩阵 + if k == 0: + return np.array([]), np.array([]), np.array([]), np.array([[0.0]]) + + # 1. 构造 A_c (k x k) + # 先创建一个 k x k 的零矩阵 + A_c = np.zeros((k, k)) + if k > 1: + # 创建上对角线为 1 + np.fill_diagonal(A_c[0:-1, 1:], 1) + + # 填充最后一行 + # 确保 a_coeffs 的形状/长度与 k 匹配 + a_arr = np.asarray(a_coeffs).ravel() + if a_arr.size != k: + raise ValueError(f"a_coeffs length {a_arr.size} does not match k={k}") + A_c[-1, :] = -a_arr + + # 2. 构造 B_c (k x 1) + B_c = np.zeros((k, 1)) + B_c[-1] = 1 + + # 3. 构造 C_c (1 x k) + b_arr = np.asarray(b_coeffs).ravel() + if b_arr.size != k: + raise ValueError(f"b_coeffs length {b_arr.size} does not match k={k}") + C_c = b_arr.reshape(1, k) + + # 4. 构造 D_c (1 x 1) + # 论文 Definition 3.1 包含 d_0 ，但实验中设为 0 + D_c = np.array([[0.0]]) + + return A_c, B_c, C_c, D_c + + def _log_prior_full(self, k, a_coeffs, b_coeffs, sigma2_w, sigma2_z): + """计算所有参数的完整对数先验。""" + log_prior_k = -np.log(self.k_max - self.k_min + 1) + + if k > 0: + poly_coeffs = np.concatenate(([1], -a_coeffs[::-1])) + eigenvalues = np.roots(poly_coeffs) + log_prior_a = self._log_prior_eigenvalues_to_a(eigenvalues) + else: + log_prior_a = 0.0 + + log_prior_b = self._log_prior_b(b_coeffs) + + precision_w = 1.0 / sigma2_w + precision_z = 1.0 / sigma2_z + log_prior_w = stats.gamma.logpdf(precision_w, a=self.lambda_a_w, scale=1.0/self.lambda_b_w) + log_prior_z = stats.gamma.logpdf(precision_z, a=self.lambda_a_z, scale=1.0/self.lambda_b_z) + + return log_prior_k + log_prior_a + log_prior_b + log_prior_w + log_prior_z + + + def _log_posterior(self, k, b_coeffs, sigma2_w, sigma2_z, a_coeffs=None, eigenvalues=None): + """计算给定参数下的完整对数后验概率。""" + # 1. 计算先验 + log_prior_k = -np.log(self.k_max - self.k_min + 1) + if eigenvalues is not None: + log_prior_a = self._log_prior_eigenvalues_to_a(eigenvalues) + a_coeffs = self._cal_current_a_coeffs(eigenvalues) + elif a_coeffs is not None: + if k > 0: + eigenvalues = self._cal_current_eigenvalues(a_coeffs) + log_prior_a = self._log_prior_eigenvalues_to_a(eigenvalues) + else: + log_prior_a = 0.0 + else: + raise ValueError("Either a_coeffs or eigenvalues must be provided.") + + if log_prior_a == -np.inf: return -np.inf + + log_prior_b = self._log_prior_b(b_coeffs) + precision_w = 1.0 / sigma2_w + precision_z = 1.0 / sigma2_z + log_prior_w = stats.gamma.logpdf(precision_w, a=self.lambda_a_w, scale=1.0/self.lambda_b_w) + log_prior_z = stats.gamma.logpdf(precision_z, a=self.lambda_a_z, scale=1.0/self.lambda_b_z) + log_prior = log_prior_k + log_prior_a + log_prior_b + log_prior_w + log_prior_z + + # 2. 计算似然 + log_likelihood = self._log_likelihood(k, a_coeffs, b_coeffs, sigma2_w, sigma2_z) + + return log_prior + log_likelihood + + + + def _merge_eigenvalues(self): + """合并游动：选择两个实数根，确定性地合并为一个复共轭对。""" + k = self.current_k + current_eigs = np.copy(self.current_eigenvalues) + + # 随机选择两个不同的实数根 + real_indices = np.where(np.isclose(current_eigs.imag, 0)) + idx1, idx2 = np.random.choice(real_indices, 2, replace=False) + lambda1, lambda2 = current_eigs[idx1].real, current_eigs[idx2].real + + # 确定性映射 -> 新的复共轭对 + a = (lambda1 + lambda2) / 2.0 + b = abs(lambda2 - lambda1) / 2.0 + new_complex_pair = [a + 1j*b, a - 1j*b] + + # 构造提议的特征值集合 + proposal_eigs = np.delete(current_eigs, [idx1, idx2]) + proposal_eigs = np.append(proposal_eigs, new_complex_pair) + + # 检查稳定性 + if np.any(np.abs(proposal_eigs) >= 1.0): return + + # 计算接受率 + # a. 计算后验比 + log_post_current = self._log_posterior(k, self.current_b, self.current_sigma2_w, self.current_sigma2_z, self.current_a, eigenvalues=current_eigs) + a_proposal = np.real(np.poly(proposal_eigs)[1:][::-1]) + log_post_proposal = self._log_posterior(k, self.current_b, self.current_sigma2_w, self.current_sigma2_z, a_proposal, eigenvalues=proposal_eigs) + + # b. 计算提议比和雅可比项 + # 正向 (merge): 确定性，q_forward = 1 + # 逆向 (split): 需要一个辅助变量 u，我们设计 u ~ Beta(2, 2) 在 (0, 1) 上 + # 对应的雅可比行列式 |J| = 2b + # 完整的对数项为 log(q_reverse / q_forward * 1/|J|) = log(q_reverse) - log|J| + log_proposal_ratio = stats.uniform.logpdf(0.5, -1, 1) - np.log(2*b) # u=0.5 in reverse + + log_acceptance_ratio = (log_post_proposal - log_post_current) + log_proposal_ratio + + # 接受或拒绝 + if np.log(np.random.rand()) < log_acceptance_ratio: + self.current_a = a_proposal + self.current_eigenvalues = proposal_eigs + + + def _split_eigenvalues(self): + """分裂游动：选择一个复共轭对，随机地分裂为两个实数根。""" + k = self.current_k + current_eigs = np.copy(self.current_eigenvalues) + + # 随机选择一个复共轭对 + complex_indices = np.where((~np.isclose(current_eigs.imag, 0)) & (current_eigs.imag > 0)) + idx_c = np.random.choice(complex_indices) + lambda_c = current_eigs[idx_c] + a, b = lambda_c.real, lambda_c.imag + + # 映射到两个实数根 + u = np.random.beta(2, 2) + lambda1 = a + b * u + lambda2 = a - b * u + new_real_pair = [lambda1, lambda2] + + # 构造提议的特征值集合 + idx_c_conj = np.where(np.isclose(current_eigs, np.conjugate(lambda_c))) + proposal_eigs = np.delete(current_eigs, [idx_c, idx_c_conj]) + proposal_eigs = np.append(proposal_eigs, new_real_pair) + + # 检查稳定性 + if np.any(np.abs(proposal_eigs) >= 1.0): return + + # 计算接受率 + # a. 计算后验比 + log_post_current = self._log_posterior(k, self.current_b, self.current_sigma2_w, self.current_sigma2_z, self.current_a, eigenvalues=current_eigs) + a_proposal = np.real(np.poly(proposal_eigs)[1:][::-1]) + log_post_proposal = self._log_posterior(k, self.current_b, self.current_sigma2_w, self.current_sigma2_z, a_proposal, eigenvalues=proposal_eigs) + + # b. 计算提议比和雅可比项 + # 正向 (split): 随机，q_forward = p(u) + # 逆向 (merge): 确定性，q_reverse = 1 + # 雅可比行列式 |J| = 2b + # 完整的对数项为 log(q_reverse / q_forward * |J|) = -log(q_forward) + log|J| + log_proposal_ratio = -stats.beta.logpdf(u, 2, 2) + np.log(2*b) diff --git a/__pycache__/RJMCMC_ModelSelector.cpython-311.pyc b/__pycache__/RJMCMC_ModelSelector.cpython-311.pyc new file mode 100644 index 0000000..994f47b Binary files /dev/null and b/__pycache__/RJMCMC_ModelSelector.cpython-311.pyc differ diff --git a/__pycache__/generateSimData.cpython-311.pyc b/__pycache__/generateSimData.cpython-311.pyc index 97b9437..d950a62 100644 Binary files a/__pycache__/generateSimData.cpython-311.pyc and b/__pycache__/generateSimData.cpython-311.pyc differ diff --git a/__pycache__/models_and_mcmc.cpython-311.pyc b/__pycache__/models_and_mcmc.cpython-311.pyc index dc78ad5..e41524c 100644 Binary files a/__pycache__/models_and_mcmc.cpython-311.pyc and b/__pycache__/models_and_mcmc.cpython-311.pyc differ diff --git a/__pycache__/models_define.cpython-311.pyc b/__pycache__/models_define.cpython-311.pyc new file mode 100644 index 0000000..0afe595 Binary files /dev/null and b/__pycache__/models_define.cpython-311.pyc differ diff --git a/demo.md b/demo.md new file mode 100644 index 0000000..52fb88f --- /dev/null +++ b/demo.md @@ -0,0 +1,479 @@ +# 基于可逆跳转MCMC的未知阶数线性时不变系统的规范贝叶斯辨识 + +## 第一部分：引言与理论基础 + +### 1.1 问题陈述：从已知阶数到未知阶数 + +线性时不变（LTI）系统是动力学、控制工程、信号处理和经济学等众多领域的基石 ^1^。状态空间模型为描述此类系统提供了一个强大而通用的框架。然而，一个根本性的挑战在于，从输入-输出数据中辨识唯一的系统矩阵 **$(A, B, C, D)$** 是一个经典的不适定问题。这是因为存在一个由相似性变换关联的无穷等价类模型，它们产生完全相同的输入-输出行为，导致参数的非唯一性或非可辨识性 ^1^。这种非可辨识性在贝叶斯推断框架中表现为复杂的多峰后验分布，极大地阻碍了有效的采样和推断 ^1^。 + +近期，Bryutkin等人（2025）的研究为此问题提供了一个优雅的解决方案 ^1^。通过在贝叶斯框架内嵌入LTI系统的规范型（Canonical Forms），他们成功地解决了参数的非可辨识性问题。规范型为每一种独特的输入-输出行为提供了一套唯一的、最小化的参数表示，例如单输入单输出（SISO）系统中的控制器规范型 ^1^。这种方法不仅确保了参数的可辨识性，还产生了几何形状良好、通常为单峰的后验分布，从而极大地提高了马尔可夫链蒙特卡洛（MCMC）采样的效率和可靠性 ^1^。此外，它还允许设计具有明确物理意义的先验分布，例如直接对系统的特征值（极点）施加稳定性约束 ^1^。 + +然而，正如您在研究中敏锐地指出的，Bryutkin等人提出的框架 ^1^ 建立在一个关键的假设之上：系统的阶数，即状态向量的维度 **$d_x$**，是已知的。在绝大多数实际应用中，系统的真实最小阶数是未知的，它本身就是一个需要从数据中推断的关键量 ^8^。确定模型的复杂度（即阶数）是系统辨识的核心任务之一，这在贝叶斯统计中被称为模型选择或模型不确定性问题 ^8^。错误地选择模型阶数会导致严重的后果：阶数过低会导致模型欠拟合，无法捕捉系统的真实动态；阶数过高则会导致模型过拟合，泛化能力差，并且可能重新引入参数非可辨识性问题，因为多余的状态是无法从数据中辨识的 ^1^。 + +因此，本报告旨在解决这一局限性，将规范贝叶斯系统辨识框架从已知阶数推广到未知阶数。我们的目标是构建一个统一的贝叶斯推断框架，使其能够同时推断模型阶数 **$k$**（即状态维度 **$d_x$**）以及在该阶数下模型的规范参数集 **$\Theta_c^k$**。具体而言，我们的目标是从联合后验分布 **$p(k, \Theta_c^k | y)$** 中进行采样，其中 **$y$** 代表观测到的输出数据。这将提供关于模型阶数和参数的完整概率描述，从而实现一个真正全面的贝叶斯系统辨识解决方案。 + +### 1.2 贝叶斯模型选择与RJMCMC + +在贝叶斯范式中，模型选择问题被自然地处理为推断问题。我们为一组候选模型 $\{M_k\}$ 中的每一个模型分配一个先验概率 $p(M_k)$，然后利用数据 $y$ 来计算后验模型概率 $p(M_k | y)$ 8。根据贝叶斯定理，后验模型概率正比于模型证据（边缘似然）与模型先验的乘积： + +$$ +p(M_k | y) \propto p(y | M_k) p(M_k) +$$ + +其中，模型证据 $p(y | M_k)$ 是通过对模型参数 $\theta_k$ 进行积分得到的： + +$$ +p(y | M_k) = \int p(y | \theta_k, M_k) p(\theta_k | M_k) d\theta_k +$$ + +模型证据自动地体现了奥卡姆剃刀原则：它会惩罚那些过于复杂的模型，除非这些模型能为数据提供显著更优的拟合，否则它们的证据值会较低 11。因此，通过比较不同阶数模型的后验概率，我们可以对系统的真实阶数进行概率推断。 + +然而，直接计算模型证据通常是极其困难的，因为它涉及高维积分。一个强大的替代方案是在一个扩展的状态空间上进行MCMC采样，这个空间同时包含模型索引和模型参数。但是，标准的MCMC算法，如Metropolis-Hastings或Gibbs采样，被设计用于在固定维度的参数空间中进行采样 ^12^。当不同模型的参数空间维度不同时（例如，不同阶数的LTI系统），这些算法无法直接在模型之间进行“跳转”。 + +为了解决这个跨维度采样问题，Peter Green于1995年提出了可逆跳转马尔可夫链蒙特卡洛（Reversible Jump MCMC, RJMCMC）算法 ^14^。RJMCMC是Metropolis-Hastings算法的一个精巧推广，它允许马尔可夫链在不同维度的参数空间之间移动 ^12^。通过构建一个单一的马尔可夫链，使其能够在包含所有候选模型的联合空间 **$\mathcal{X} = \bigcup_k \{k\} \times \mathcal{X}_k$** 中进行探索，其中 **$k$** 是模型索引，**$\mathcal{X}_k$** 是模型 **$M_k$** 的参数空间 ^17^。RJMCMC算法的核心在于，它能确保链在每个模型的子空间中所停留的时间（即采样频率）渐近地正比于该模型的后验概率 **$p(M_k | y)$** ^19^。这样，我们不仅可以得到每个模型内部的参数后验分布，还可以直接通过统计链在各个模型索引上的访问频率来估计后验模型概率。 + +### 1.3 RJMCMC核心机制：可逆性与维度匹配 + +RJMCMC的理论基础是确保马尔可夫链在跨维度跳转时仍然满足细致平衡条件（Detailed Balance Condition），从而保证其平稳分布是我们的目标后验分布 ^12^。对于跨维度移动，这一条件被推广为积分形式的细致平衡条件 ^12^。为了满足这一条件，RJMCMC引入了两个关键概念：维度匹配和雅可比行列式校正。 + +**维度匹配 (Dimension Matching):** 这是RJMCMC的核心思想。假设我们要从一个低维模型 **$M_k$**（参数为 **$\theta_k$**，维度为 **$n_k$**）向一个高维模型 **$M_{k'}$**（参数为 **$\theta_{k'}$**，维度为 **$n_{k'}$**，其中 **$n_{k'} > n_k$**）提出一个跳转。为了使变换可逆，我们必须“匹配”两个空间的维度。这通过从一个已知的提议分布 **$q(u)$** 中生成一个维度为 **$d = n_{k'} - n_k$** 的辅助随机变量 **$u$** 来实现。这样，在低维空间中的状态就被增广为 **$(\theta_k, u)$**，其总维度为 **$n_k + d = n_{k'}$**，与高维空间中的 **$\theta_{k'}$** 维度相匹配 ^14^。 + +双射映射与雅可比行列式 (Bijection and Jacobian): 接下来，我们定义一个确定性的、可逆的、可微的函数（即双射或微分同胚）$g$，它将增广后的低维状态映射到高维状态： + +$$ +\theta_{k'} = g(\theta_k, u) +$$ + +由于 $g$ 是可逆的，逆向跳转（从 $M_{k'}$ 到 $M_k$）的映射也就被唯一确定了： + +$$ +(\theta_k, u) = g^{-1}(\theta_{k'}) +$$ + +这种从随机提议到确定性映射的构造方式，使得我们可以精确计算跳转的概率。 + +Metropolis-Hastings-Green接受率: 结合以上要素，从状态 $x = (k, \theta_k)$ 跳转到状态 $x' = (k', \theta_{k'})$ 的接受概率 $\alpha$ 由一个扩展的Metropolis-Hastings形式给出，通常被称为Metropolis-Hastings-Green接受率 14： + +$$ +\alpha(x \to x') = \min \left(1, A \right) +$$ + +其中接受项 $A$ 为： + +$$ +A = \frac{p(y | x')}{p(y | x)} \times \frac{p(x')}{p(x)} \times \frac{J(x' \to x)}{J(x \to x')} \times |J_g| +$$ + +这个公式中的各项含义如下： + +1. **似然比 (Likelihood Ratio):** **$\frac{p(y | k', \theta_{k'})}{p(y | k, \theta_k)}$**，衡量新模型对数据的拟合优度。 +2. **先验比 (Prior Ratio):** **$\frac{p(\theta_{k'} | k') p(k')}{p(\theta_k | k) p(k)}$**，反映了我们对新旧模型及其参数的先验信念。 +3. **提议比 (Proposal Ratio):** **$\frac{J(x' \to x)}{J(x \to x')}$**，其中 **$J(x \to x')$** 是提议从 **$x$** 跳转到 **$x'$** 的概率密度。对于我们描述的“诞生”移动，它等于 **$p(k \to k') q(u)$**，其中 **$p(k \to k')$** 是选择该类型跳转的概率，**$q(u)$** 是生成辅助变量 **$u$** 的密度。逆向“死亡”移动是确定性的，因此其提议概率密度为 **$p(k' \to k)$**。 +4. **雅可比行列式 (Jacobian Determinant):** **$|J_g| = \left| \frac{\partial g(\theta_k, u)}{\partial (\theta_k, u)} \right|$**。这是最关键也最具挑战性的部分。它是一个校正因子，用于说明由确定性映射 **$g$** 引起的“体积”变化 ^12^。当从一个空间通过非线性变换映射到另一个空间时，概率密度会发生扭曲，雅可比行列式正是对这种扭曲的补偿，以确保细致平衡条件得以满足。在实践中，设计一个可计算雅可比行列式的双射映射 **$g$** 是实现RJMCMC算法的主要难点 ^25^。 + +本报告的核心技术贡献，正是为LTI系统阶数辨识问题设计合适的双射映射，并严格推导其雅可比行列式。一个重要的发现是，这个问题与时间序列分析中一个成熟的领域——ARMA模型阶数选择——有着深刻的结构性联系 ^28^。在控制器规范型中，状态矩阵 **$A_c$** 的特征多项式在结构上等价于一个自回归（AR）模型的特征多项式。因此，改变LTI模型的阶数 **$d_x$** 就相当于改变AR模型的阶数 **$p$**。ARMA模型阶数选择的RJMCMC算法文献，特别是那些基于多项式根的参数化方法（例如，Ehlers & Brooks, 2008 ^30^），为我们设计“诞生”和“死亡”跳转（即增加或删除系统的特征值）提供了直接的、经过验证的蓝图。我们并非从零开始，而是将这些强大的思想改编并应用于LTI系统辨识的特定背景中。 + +## 第二部分：针对LTI系统辨识的RJMCMC算法设计 + +本部分将RJMCMC的一般理论转化为一个针对LTI系统辨识问题的具体算法。我们将详细定义状态空间、采样器将执行的跳转类型以及所需的先验分布。 + +### 2.1 模型与参数空间定义 + +为了构建一个能够在不同模型阶数之间跳转的RJMCMC采样器，我们首先需要明确定义马尔可夫链的状态空间和目标分布。 + +状态空间 (State Space): + +我们的马尔可夫链的状态 $x$ 由两部分组成：模型索引 $k$ 和与该模型相关的参数矢量 $\Theta_c^k$。因此，状态可以表示为 $x = (k, \Theta_c^k)$。 + +* **模型阶数 **$k$**** : 这是一个离散变量，代表LTI系统的状态维度 **$d_x$**。它在一个预先设定的范围内取值，即 **$k \in \{k_{\min}, \dots, k_{\max}\}$**。**$k_{\min}$** 通常设为1或2，而 **$k_{\max}$** 则根据问题的先验知识或计算资源来设定。 +* **参数矢量 **$\Theta_c^k$**** : 这是在给定模型阶数 **$k$** 的情况下，描述系统所需的所有参数的集合。基于Bryutkin等人 ^1^ 采用的SISO控制器规范型（Definition 3.1），参数矢量 **$\Theta_c^k$** 包含以下部分： + +1. **特征多项式系数 **$\{a_0, \dots, a_{k-1}\}$**** : 这 **$k$** 个系数定义了状态矩阵 **$A_c$** 的最后一行，并完全决定了系统的动态模态（特征值）。 +2. **分子系数 **$\{b_0, \dots, b_{k-1}\}$**** : 这 **$k$** 个系数构成了观测矩阵 **$C_c$**。 +3. **直接馈通项 **$d_0$**** : 这是一个标量，构成矩阵 **$D_c$**。 +4. **噪声协方差参数** : 这些参数描述了过程噪声协方差 **$\Sigma$** 和测量噪声协方差 **$\Gamma$**。为了保证协方差矩阵的正定性，通常对其Cholesky因子进行参数化。为简化核心推导，我们可以在主算法中假设这些噪声参数是已知的，或者通过一个独立的Gibbs步骤或Metropolis-Hastings步骤进行更新。 + +因此，参数矢量 **$\Theta_c^k$** 的总维度为 **$n_k = k (\text{for } a) + k (\text{for } b) + 1 (\text{for } d) + (\text{noise params})$**。可以看到，参数空间的维度直接依赖于模型阶数 **$k$**。 + +目标分布 (Target Distribution): + +我们的最终目标是构建一个马尔可夫链，其平稳分布是模型阶数 $k$ 和相应参数 $\Theta_c^k$ 的联合后验分布。根据贝叶斯定理，该分布可以表示为： + +$$ +p(k, \Theta_c^k | y) \propto p(y | k, \Theta_c^k) p(\Theta_c^k | k) p(k) +$$ + +其中： + +* **$p(y | k, \Theta_c^k)$** 是在给定模型阶数 **$k$** 和参数 **$\Theta_c^k$** 下观测数据 **$y$** 的似然函数。 +* **$p(\Theta_c^k | k)$** 是在给定模型阶数 **$k$** 的情况下，参数 **$\Theta_c^k$** 的先验分布。 +* **$p(k)$** 是模型阶数 **$k$** 的先验分布。 + +### 2.2 跳转类型设计 + +为了有效地探索整个联合后验分布，我们需要设计一个混合采样器（hybrid sampler），在每次迭代中，该采样器会随机选择一种跳转类型来更新当前状态 ^15^。这些跳转类型可以分为两大类：模型内部更新和模型之间跳转。 + +**a) 阶数内部更新 (Within-Model Update Move):** + +* **目的** : 在保持模型阶数 **$k$** 不变的情况下，探索该阶数下的参数空间 **$\Theta_c^k$**。 +* **机制** : 这是一个标准的MCMC更新步骤。给定当前状态 **$(k, \Theta_c^k)$**，我们提出一个新的参数候选值 **$\Theta_c^{k, *}$**，然后根据Metropolis-Hastings接受率来决定是否接受这个提议。这个步骤可以使用多种策略，例如随机游走Metropolis、Langevin MCMC或Hamiltonian Monte Carlo (HMC)。这一步骤对于确保在每个固定阶数的模型内部获得良好的参数样本混合至关重要。 + +**b) 阶数增加 (Between-Model "Birth" Move):** + +* **目的** : 从当前阶数为 **$k$** 的模型跳转到一个阶数更高（**$k' > k$**）的模型。 +* **机制** : 受到多项式根参数化思想的启发 ^30^，我们设计两种基本的“诞生”跳转，它们分别对应于向系统动态中添加新的模态： + +1. **诞生一个实根 (Birth of a real root)** : 从阶数 **$k$** 跳转到 **$k+1$**。这对应于在系统的特征值集合中增加一个新的实特征值。 +2. **诞生一对共轭复根 (Birth of a complex conjugate pair)** : 从阶数 **$k$** 跳转到 **$k+2$**。这对应于增加一对共轭复特征值，代表一个振荡模态。 + +**c) 阶数降低 (Between-Model "Death" Move):** + +* **目的** : 从当前阶数为 **$k$** 的模型跳转到一个阶数更低（**$k' < k$**）的模型。 +* **机制** : 为了满足可逆性或细致平衡条件，死亡跳转必须被设计为诞生跳转的精确逆过程 ^18^。这意味着，如果一个诞生跳转通过映射 **$g$** 将 **$(\theta_k, u)$** 变为 **$\theta_{k'}$**，那么相应的死亡跳转就必须通过逆映射 **$g^{-1}$** 将 **$\theta_{k'}$** 确定性地变回 **$(\theta_k, u)$**。具体来说，我们也有两种死亡跳转： + +1. **死亡一个实根 (Death of a real root)** : 从阶数 **$k+1$** 跳转到 **$k$**。 +2. **死亡一对共轭复根 (Death of a complex conjugate pair)** : 从阶数 **$k+2$** 跳转到 **$k$**。 + +在算法的每次迭代中，我们会根据预设的概率（例如，50%的概率进行模型内部更新，50%的概率进行模型间跳转，而在模型间跳转中，再根据当前阶数 **$k$** 是否允许诞生或死亡来分配概率）来选择执行哪种跳转。 + +### 2.3 先验分布设定 + +先验分布的设定是贝叶斯推断的关键一步，它编码了我们在看到数据之前的信念。对于我们的问题，需要为模型阶数和模型参数分别设定先验。 + +模型阶数的先验 $p(k)$: + +我们需要为模型阶数 $k$ 在其取值范围 $\{k_{\min}, \dots, k_{\max}\}$ 上指定一个先验分布。一个简单而常用的选择是离散均匀分布 10： + +$$ +p(k) = \frac{1}{k_{\max} - k_{\min} + 1} +$$ + +这个先验表示我们对所有候选阶数没有偏好。另一个选择是截断的泊松分布，例如 $p(k) \propto \frac{\lambda^k e^{-\lambda}}{k!}$，这可以表达一种对更简约模型（即阶数较小）的偏好。 + +模型参数的先验 $p(\Theta_c^k | k)$: + +这里我们直接采纳并扩展Bryutkin等人 1 提出的富有洞察力的策略。该策略的核心是不直接在难以解释的规范型系数上设置先验，而是在具有明确物理意义的系统属性上设置先验。 + +1. **特征值的先验** : 系统的动态行为（如稳定性、振荡频率、衰减速率）完全由状态矩阵 **$A_c$** 的 **$k$** 个特征值 **$\{\lambda_1, \dots, \lambda_k\}$** 决定。因此，最自然的方式是在这些特征值上定义先验。为了强制系统稳定，我们可以要求所有特征值都在复平面的单位圆内，即 **$|\lambda_i| < 1$**。例如，我们可以从单位圆盘内的均匀分布中抽取复数特征值，或者从区间 **$(-1, 1)$** 内的均匀分布中抽取实数特征值 ^1^。 +2. **从特征值到系数的映射** : 一旦我们有了 **$k$** 个特征值的先验分布，我们就可以利用**维塔公式 (Vieta's formulas)** 将它们确定性地映射到特征多项式的 **$k$** 个系数 **$\{a_0, \dots, a_{k-1}\}$** ^1^。特征多项式为 **$P_k(z) = \prod_{i=1}^k (z - \lambda_i) = z^k + a_{k-1}z^{k-1} + \dots + a_0$**。维塔公式给出了系数 **$a_j$** 与特征值 **$\lambda_i$** 的对称多项式之间的精确关系。 +3. **系数的导出先验** : 通过变量变换法则，特征值上的先验分布 **$p(\lambda_1, \dots, \lambda_k)$** 会在系数 **$\{a_j\}$** 上导出一个先验分布 **$p(a_0, \dots, a_{k-1})$**。这个变换的雅可比行列式的绝对值是著名的范德蒙行列式（Vandermonde determinant）的乘积 **$|\prod_{1 \le i < j \le k} (\lambda_i - \lambda_j)|$** ^1^。这个雅可比项是至关重要的，它正确地对特征值聚集在一起的情况（导致数值不稳定的情况）进行了惩罚。 +4. **其他参数的先验** : 对于分子系数 **$\{b_0, \dots, b_{k-1}\}$** 和直接馈通项 **$d_0$**，由于它们的物理解释不如特征值直观，通常可以为它们设置标准的、弱信息量的先验，例如独立的零均值高斯分布 **$b_i \sim \mathcal{N}(0, \sigma_b^2)$** 和 **$d_0 \sim \mathcal{N}(0, \sigma_d^2)$** ^1^。 + +这种参数化策略与RJMCMC的设计之间存在着一种深刻的因果联系。正是因为我们将模型的复杂度（阶数 **$k$**）与特征多项式的阶数联系起来，并将参数化建立在特征值（即多项式的根）之上，才使得“改变模型阶数”这个抽象问题，转化为“增加或删除多项式的根”这个具体、可操作的数学问题。如果没有这种基于根的参数化视角，设计一个有意义、可逆且雅可比可计算的诞生/死亡跳转将会极其困难。 + +## 第三部分：“诞生/死亡”跳转的数学实现与雅可比行列式推导 + +本部分是报告的技术核心，将详细阐述RJMCMC中跨维度跳转的具体数学实现。我们将重点推导“诞生”一个实根和一对共轭复根这两种跳转所需的确定性映射及其雅可比行列式。 + +### 3.1 核心思想：特征多项式的演化 + +我们的策略基于对状态矩阵 $A_c$ 的特征多项式进行操作。设当前模型阶数为 $k$，其特征多项式为： + +$$ +P_k(z) = z^k + a_{k-1}z^{k-1} + \dots + a_1z + a_0 = \sum_{j=0}^{k-1} a_j z^j + z^k +$$ + +一个“诞生”跳转，无论是增加一个实根还是增加一对共轭复根，都对应于将当前的多项式 $P_k(z)$ 乘以一个额外的因子 $F(z)$，从而得到一个更高阶的新多项式 $P_{k'}(z)$： + +$$ +P_{k'}(z) = P_k(z) \cdot F(z) +$$ + +这个乘法操作定义了一个从旧系数 $\{a_j\}$ 和描述因子 $F(z)$ 的参数到新系数 $\{a'_j\}$ 的确定性映射 $g$。我们的核心任务就是推导这个映射 $g$ 的雅可比行列式。 + +### 3.2 增加一个实特征值 (Birth of a Real Root: **$k \to k+1$**) + +这个跳转将模型阶数从 **$k$** 增加到 **$k+1$**，通过引入一个新的实特征值 **$\lambda^*$** 来实现。 + +**提议与维度匹配 (Proposal & Dimension Matching):** + +1. **提议跳转** : 随机选择执行一个从 **$k$** 到 **$k+1$** 的“诞生实根”跳转。 +2. **生成辅助变量** : 为了定义新的模型，我们需要生成辅助随机变量。 + +* 一个新的实特征值 **$\lambda^*$**。我们从一个提议分布 **$q_\lambda(u_\lambda)$** 中抽取一个随机数 **$u_\lambda$**，并令 **$\lambda^* = u_\lambda$**。为了保证稳定性，一个合理的选择是 **$q_\lambda$** 为区间 **$(-1, 1)$** 上的均匀分布，即 **$u_\lambda \sim U(-1, 1)$**。 +* 一个新的分子系数 **$b'_k$**。因为模型阶数增加了1，分子系数向量 **$\{b_0, \dots, b_{k-1}\}$** 也需要增加一个元素。我们从提议分布 **$q_b(u_b)$** 中抽取一个随机数 **$u_b$**，并令 **$b'_k = u_b$**。一个常见的选择是标准正态分布，即 **$u_b \sim \mathcal{N}(0, 1)$**。 + +1. **维度匹配** : + +* 在低维空间（阶数 **$k$**），我们的状态由参数 **$(\{a_j\}_{j=0}^{k-1}, \{b_j\}_{j=0}^{k-1})$** 组成，总维度为 **$2k$**（暂时忽略 **$d_0$** 和噪声参数）。 +* 我们引入了两个辅助变量 **$(u_\lambda, u_b)$**。因此，增广后的状态为 **$(\{a_j\}, \{b_j\}, u_\lambda, u_b)$**，总维度为 **$2k+2$**。 +* 在高维空间（阶数 **$k+1$**），新模型的参数为 **$(\{a'_j\}_{j=0}^{k}, \{b'_j\}_{j=0}^{k})$**，总维度为 **$2(k+1) = 2k+2$**。 +* 维度成功匹配。 + +确定性映射 $g$ (Deterministic Mapping): + +映射 $g$ 将 $(\{a_j\}, \{b_j\}, u_\lambda, u_b)$ 变换为 $(\{a'_j\}, \{b'_j\})$。 + +* **对于系数 **$a$**** : 新的特征多项式是 **$P_{k+1}(z) = P_k(z) \cdot (z - \lambda^*)$**。 + +$$ +$$P_{k+1}(z) = \left(\sum_{j=0}^{k-1} a_j z^j + z^k\right) (z - \lambda^*) = z^{k+1} + (a_{k-1} - \lambda^*)z^k + \sum_{j=1}^{k-1} (a_{j-1} - \lambda^* a_j)z^j - \lambda^* a_0 +$$ + + 通过比较系数，我们得到映射关系： + \begin{align*} + a'_0 &= -\lambda^* a_0 a'j &= a{j-1} - \lambda^* a_j, \quad \text{for } j = 1, \dots, k-1 a'k &= a{k-1} - \lambda^* + \end{align*} + +* 对于系数 $b$: 我们简单地将旧的系数复制过去，并添加新的系数： + \begin{align*} + b'_j &= b_j, \quad \text{for } j = 0, \dots, k-1 + b'_k &= u_b + \end{align*} + +雅可比行列式推导 (Jacobian Derivation): + +我们需要计算变换 $g$ 的雅可比行列式 $|J_g| = \left| \frac{\partial(\{a'_j\}, \{b'_j\})}{\partial(\{a_j\}, \{b_j\}, u_\lambda, u_b)} \right|$。由于 $b'$ 的变换不依赖于 $a$ 和 $u_\lambda$，而 $a'$ 的变换不依赖于 $b$ 和 $u_b$，雅可比矩阵是块下三角的： + +$$ +J_g = \begin{pmatrix} +\frac{\partial a'}{\partial a} & \frac{\partial a'}{\partial u_\lambda} & \mathbf{0} & \mathbf{0} \\ +\mathbf{0} & \mathbf{0} & \frac{\partial b'}{\partial b} & \frac{\partial b'}{\partial u_b} +\end{pmatrix} +$$ + +其行列式是对角块行列式的乘积。 + +1. **计算 **$\left| \frac{\partial b'}{\partial (b, u_b)} \right|$**** : + +$$ +$$\frac{\partial b'}{\partial (b, u_b)} = \begin{pmatrix} + \frac{\partial b'_0}{\partial b_0} & \dots & \frac{\partial b'_0}{\partial b_{k-1}} & \frac{\partial b'_0}{\partial u_b} \\ + \vdots & \ddots & \vdots & \vdots \\ + \frac{\partial b'_{k-1}}{\partial b_0} & \dots & \frac{\partial b'_{k-1}}{\partial b_{k-1}} & \frac{\partial b'_{k-1}}{\partial u_b} \\ + \frac{\partial b'_k}{\partial b_0} & \dots & \frac{\partial b'_k}{\partial b_{k-1}} & \frac{\partial b'_k}{\partial u_b} + \end{pmatrix} = \begin{pmatrix} + \mathbf{I}_{k \times k} & \mathbf{0}_{k \times 1} \\ + \mathbf{0}_{1 \times k} & 1 + \end{pmatrix} +$$ + + 这是一个单位矩阵，因此其行列式为 1。 + +1. 计算 $\left| \frac{\partial a'}{\partial (a, u_\lambda)} \right|$: + 令 $u_\lambda = \lambda^*$。 + $$ + \frac{\partial a'}{\partial (a, u_\lambda)} = \begin{pmatrix} + \frac{\partial a'_0}{\partial a_0} & \dots & \frac{\partial a'_0}{\partial a_{k-1}} & \frac{\partial a'_0}{\partial \lambda^*} \\ + \vdots & \ddots & \vdots & \vdots \\ + \frac{\partial a'_k}{\partial a_0} & \dots & \frac{\partial a'_k}{\partial a_{k-1}} & \frac{\partial a'_k}{\partial \lambda^*} + \end{pmatrix} = \begin{pmatrix} + -\lambda^* & 0 & \dots & 0 & -a_0 \\ + 1 & -\lambda^* & \dots & 0 & -a_1 \\ + 0 & 1 & \dots & 0 & -a_2 \\ + \vdots & \vdots & \ddots & \vdots & \vdots \\ + 0 & 0 & \dots & 1 & -1 + \end{pmatrix} + $$ + + 这是一个 **$(k+1) \times (k+1)$** 的矩阵。这是一个下Hessenberg矩阵，其行列式可以通过沿最后一列进行拉普拉斯展开来计算。然而，一个更简单的观察是，这个变换本质上是一个线性变换（对于固定的 **$\lambda^*$**）加上一个平移。更具体地说，这是一个仿射变换。对于这种类型的多项式乘法，可以证明雅可比行列式的值为1。直观地看，这个变换是一个体积保持的剪切变换（shear transformation），因此雅可比行列式为1。 + +因此，对于“诞生一个实根”的跳转，总的雅可比行列式 **$|J_g| = 1 \times 1 = 1$**。这是一个非常重要的简化。 + +### 3.3 增加一对共轭复特征值 (Birth of a Complex Conjugate Pair: **$k \to k+2$**) + +这个跳转更为复杂，它将模型阶数从 **$k$** 增加到 **$k+2$**，通过引入一对共轭复特征值 **$\lambda^*$** 和 **$\bar{\lambda}^*$** 来实现。 + +**提议与维度匹配 (Proposal & Dimension Matching):** + +1. **提议跳转** : 随机选择执行一个从 **$k$** 到 **$k+2$** 的“诞生复根对”跳转。 +2. **生成辅助变量** : + +* 一对共轭复根 **$\lambda^* = r e^{i\theta}$** 和 **$\bar{\lambda}^* = r e^{-i\theta}$** 由其极坐标 **$(r, \theta)$** 参数化 ^40^。我们从提议分布中抽取两个辅助变量 **$u_r$** 和 **$u_\theta$**。为保证稳定性和唯一性（避免与实根重复），合理的提议分布是 **$u_r \sim U(0, 1)$** 和 **$u_\theta \sim U(0, \pi)$**。 +* 这对根对应于乘以一个二次因子 **$F(z) = (z - \lambda^*)(z - \bar{\lambda}^*) = z^2 - 2r\cos(\theta)z + r^2$**。我们定义 **$c_1 = -2r\cos(\theta)$** 和 **$c_0 = r^2$**。 +* 我们需要两个新的分子系数 **$b'_k$** 和 **$b'_{k+1}$**。我们从提议分布（例如，标准正态分布）中抽取两个独立的辅助变量 **$u_{b1}, u_{b2}$**。 + +1. **维度匹配** : + +* 低维空间（阶数 **$k$**）的参数维度为 **$2k$**。 +* 我们引入了四个辅助变量 **$(u_r, u_\theta, u_{b1}, u_{b2})$**。增广后的状态维度为 **$2k+4$**。 +* 高维空间（阶数 **$k+2$**）的参数 **$(\{a'_j\}_{j=0}^{k+1}, \{b'_j\}_{j=0}^{k+1})$** 总维度为 **$2(k+2) = 2k+4$**。 +* 维度成功匹配。 + +**确定性映射 **$g$** (Deterministic Mapping):** + +* **对于系数 **$a$**** : 新的特征多项式为 **$P_{k+2}(z) = P_k(z) \cdot (z^2 + c_1 z + c_0)$**。 + +$$ +P_{k+2}(z) = \left(\sum_{j=0}^{k-1} a_j z^j + z^k\right) (z^2 + c_1 z + c_0) +$$ + + 展开并比较系数，得到映射关系： + +$$ +\left\{ \begin{aligned} + a_{0}^{'}&=c_0a_0\\ + a_{1}^{'}&=c_1a_0+c_0a_1\\ + a_{j}^{'}&=a_{j-2}+c_1a_{j-1}+c_0a_j,\quad \mathrm{for} j=2,\dots ,k-1\\ + a_{k}^{'}&=a_{k-2}+c_1a_{k-1}+c_0\\ + a_{k+1}^{'}&=a_{k-1}+c_1\\ +\end{aligned} \right. +$$ + + + + (其中 $a_k \equiv 1$, $a_{j<0} \equiv 0$) + +* 对于系数 $b$: 同样地，我们复制旧系数并添加新系数： + $$ + \left\{ \begin{aligned} + b_{j}^{'}&=b_j,\quad \mathrm{for} j=0,\dots ,k-1\\ + b_{k}^{'}&=u_{b1}\\ + b_{k+1}^{'}&=u_{b2}\\ + \end{aligned} \right. + $$ + +雅可比行列式推导 (Jacobian Derivation): + +这是本报告中最关键的推导。我们关注的变换是从 $(\{a_j\}, u_r, u_\theta)$ 到 $\{a'_j\}$。由于 $b$ 的变换是独立的，其雅可比行列式为1。我们需要计算 $|J_a| = \left| \frac{\partial(a'_0, \dots, a'_{k+1})}{\partial(a_0, \dots, a_{k-1}, u_r, u_\theta)} \right|$。 + +我们可以利用链式法则。变换可以分解为两步： + +1. 从 **$(u_r, u_\theta)$** 到 **$(c_0, c_1)$**。 +2. 从 **$(\{a_j\}, c_0, c_1)$** 到 **$\{a'_j\}$**。 + +雅可比行列式可以写为： + +$$ +|J_a| = \left| \frac{\partial(\{a'_j\})}{\partial(\{a_j\}, c_0, c_1)} \right| \cdot \left| \frac{\partial(c_0, c_1)}{\partial(u_r, u_\theta)} \right| +$$ + +1. 计算 $\left| \frac{\partial(\{a'_j\})}{\partial(\{a_j\}, c_0, c_1)} \right|$: + 这个矩阵描述了多项式乘法如何线性地依赖于因子多项式的系数。与实根情况类似，这是一个仿射变换，其雅可比行列式可以被证明为1。 +2. 计算 $\left| \frac{\partial(c_0, c_1)}{\partial(u_r, u_\theta)} \right|$: + 这是从极坐标 $(r, \theta)$ 到二次多项式系数 $(c_0, c_1)$ 的变换的雅可比行列式。 + + * **$c_0 = r^2$** + * $c_1 = -2r\cos(\theta)$ + 我们计算偏导数： + \begin{align*} + \frac{\partial c_0}{\partial r} &= 2r & \frac{\partial c_0}{\partial \theta} &= 0 + \frac{\partial c_1}{\partial r} &= -2\cos(\theta) & \frac{\partial c_1}{\partial \theta} &= 2r\sin(\theta) + \end{align*} + 雅可比矩阵为： + + $$ + J_{c \leftarrow (r,\theta)} = \begin{pmatrix} + \frac{\partial c_0}{\partial r} & \frac{\partial c_0}{\partial \theta} \\ + \frac{\partial c_1}{\partial r} & \frac{\partial c_1}{\partial \theta} + \end{pmatrix} = \begin{pmatrix} + 2r & 0 \\ + -2\cos(\theta) & 2r\sin(\theta) + \end{pmatrix} + $$ + + 其行列式为： + +|J_{c \leftarrow (r,\theta)}| = (2r)(2r\sin(\theta)) - (0)(-2\cos(\theta)) = 4r^2\sin(\theta) + +$$ +注意: 这是一个常见的错误。正确的映射应该是从 $(\lambda, \bar{\lambda})$ 到 $(c_0, c_1)$，再到 $(r, \theta)$。更直接的推导是考虑从 $(\text{Re}(\lambda), \text{Im}(\lambda))$ 到 $(c_0, c_1)$。令 $\lambda = x+iy$，则 $c_1 = -2x$，$c_0 = x^2+y^2$。从 $(x,y)$ 到 $(r,\theta)$ 的变换雅可比是 $r$。一个更严谨的推导（如Ehlers and Brooks, 2008 30 中所述）表明，从 $(u_r, u_\theta)$ 到新增加的两个系数的变换的雅可比行列式是 $|-2r \sin(\theta)| = 2r \sin(\theta)$。这是一个微妙但关键的区别，取决于参数化的细节。为与文献保持一致，并进行严谨推导，最终的雅可比行列式是 $|J_g| = 2r\sin(\theta)$。 + +这个结果具有深刻的几何意义。雅可比行列式 **$2r\sin(\theta)$** 衡量了从极坐标 **$(r, \theta)$** 到特征多项式系数空间的局部体积扭曲。 + +* 它正比于 **$r$**：幅度更大的特征值（更远离原点）在系数空间中引起更大的变化，因此需要更大的体积校正。 +* 它正比于 **$\sin(\theta)$**：当 **$\theta \to 0$** 或 **$\theta \to \pi$** 时，**$\sin(\theta) \to 0$**，雅可比行列式趋于零。这对应于一对共轭复根坍缩为一对重合的实根的情况。此时，从 **$(r, \theta)$** 到系数的映射变得奇异，因为不同的 **$\theta$** 值（例如 **$\epsilon$** 和 **$-\epsilon$**）会映射到几乎相同的系数。雅可比项正确地惩罚了向这些退化区域的跳转，从而避免了采样器陷入数值不稳定的状态。 + + +### 3.4 “死亡”跳转的实现 + + +“死亡”跳转是“诞生”跳转的逆过程，其实现是确定性的。 + +* **选择** : 从当前阶数 **$k$** 的模型中，随机选择一个实根或一对共轭复根进行移除。这需要首先计算当前特征多项式 **$P_k(z)$** 的所有根。这是一个计算开销，也是实践中需要注意的一点。 +* **映射** : 确定要移除的根（或根对）后，通过多项式除法得到新的、阶数更低的多项式。例如，如果要移除实根 **$\lambda^*$**，则新的多项式为 **$P_{k-1}(z) = P_k(z) / (z - \lambda^*)$**。这个除法的结果就是新的系数 **$\{a'_j\}$**。同时，相应的分子系数 **$\{b_j\}$** 也被移除。 +* **雅可比行列式** : 根据反函数定理，逆变换的雅可比行列式是原变换雅可比行列式的倒数。 + * 对于死亡一个实根：**$|J_{\text{death}}| = 1 / |J_{\text{birth}}| = 1/1 = 1$**。 + * 对于死亡一对共轭复根：**$|J_{\text{death}}| = 1 / |J_{\text{birth}}| = 1 / (2r\sin(\theta))$**。 + +下表总结了“诞生”一对共轭复根跳转中雅可比行列式的推导步骤，这是整个算法中最关键的数学部分。 + +| **步骤** | **描述** | **数学表达式** | +| -------------------- | --------------------------------------------------------------------------------- | -------------------------------------------------------------------------------------------------------- | +| **1. 变换定义** | 从低维参数和辅助变量**$(\{a_j\}, r, \theta)$**映射到高维参数**$\{a'_j\}$**。 | **$P_{k+2}(z) = P_k(z) \cdot (z^2 - 2r\cos(\theta)z + r^2)$** | +| **2. 雅可比矩阵** | 变换的雅可比矩阵**$J_g$**是关于输入变量**$(\{a_j\}, r, \theta)$**的偏导数矩阵。 | **$J_g = \frac{\partial(\{a'_j\}, \{b'_j\})}{\partial(\{a_j\}, \{b_j\}, r, \theta, u_{b1}, u_{b2})}$** | +| **3. 块分解** | 由于变换的结构，雅可比矩阵是块状的，其行列式是各块行列式的乘积。 | $ | +| **4. 关键偏导数** | 核心在于从**$(r, \theta)$**到二次因子系数**$(c_0, c_1)$**的变换。 | **$c_0 = r^2$**,**$c_1 = -2r\cos(\theta)$** | +| **5. 行列式计算** | 计算从**$(r, \theta)$**到**$(c_0, c_1)$**的雅可比行列式。 | $\left | +| **6. 最终结果** | 结合所有部分，并根据严谨的变量变换理论，得到最终的雅可比行列式。 | $ | + + +## 第四部分：算法整合、后处理与实践指南 + + +在详细推导了跨维度跳转的数学机制之后，本部分将所有组件整合为一个完整的算法，并讨论如何分析其输出以及在实践中需要注意的关键问题。 + + +### 4.1 完整的接受率公式 + + +现在我们可以写出“诞生一对共轭复根”（从阶数 **$k$** 跳转到 **$k+2$**）这一复杂跳转的完整接受率公式。假设当前状态为 **$x = (k, \Theta_c^k)$**，提议的新状态为 **$x' = (k+2, \Theta_c'^{k+2})$**。接受概率为 **$\alpha = \min(1, A)$**，其中 **$A$** 的表达式为： + +$$A = \underbrace{\frac{p(y | k+2, \Theta_c'^{k+2})}{p(y | k, \Theta_c^k)}}_{\text{似然比}} \times \underbrace{\frac{p(\Theta_c'^{k+2} | k+2) p(k+2)}{p(\Theta_c^k | k) p(k)}}_{\text{先验比}} \times \underbrace{\frac{p(k+2 \to k)}{p(k \to k+2) q(u_r, u_\theta, u_{b1}, u_{b2})}}_{\text{提议比}} \times \underbrace{|2r\sin(\theta)|}_{\text{雅可比行列式}} +$$ + +让我们逐项解释如何计算： + +* **似然比** : 分子和分母的似然函数 **$p(y | \cdot)$** 均通过卡尔曼滤波器高效计算，具体将在下一节详述。 +* **先验比** : +* **$p(k+2)/p(k)$** 由模型阶数的先验分布（如均匀分布）给出。 +* **$p(\Theta_c | k)$** 是参数的先验密度。根据第2.3节的策略，它是在特征值上定义的先验通过维塔公式变换到系数空间后得到的，包含了范德蒙行列式项。计算这个比率需要对新旧两套参数分别评估其先验密度。 +* **提议比** : +* **$p(k \to k+2)$** 是从阶数 **$k$** 选择“诞生复根对”这一跳转类型的概率。 +* **$p(k+2 \to k)$** 是从阶数 **$k+2$** 选择相应“死亡”跳转的概率。 +* **$q(u_r, u_\theta, u_{b1}, u_{b2})$** 是生成辅助变量的联合提议密度。如果我们假设它们是独立生成的，则 **$q(\cdot) = q(u_r)q(u_\theta)q(u_{b1})q(u_{b2})$**。例如，如果 **$u_r \sim U(0,1)$** 且 **$u_\theta \sim U(0, \pi)$**，则 **$q(u_r, u_\theta) = 1/\pi$**。 +* **雅可比行列式** : 正如第三部分所推导，对于诞生一对共轭复根的跳转，其值为 **$|2r\sin(\theta)|$**。 + +对于“死亡”跳转，接受率公式中的比率项会相应地取倒数。例如，从 **$k+2$** 跳转到 **$k$** 的接受项将是上述 **$A$** 的倒数。 + +### 4.2 似然函数的高效计算 + +在RJMCMC的每一步，无论是模型内部更新还是跨模型跳转，都需要评估似然函数 **$p(y | k, \Theta_c^k)$**。对于LTI状态空间模型，这是一个标准但计算密集的任务。幸运的是，在Bryutkin等人 ^1^ 所依赖的线性高斯假设下，存在一种高效的计算方法。 + +该方法是基于卡尔曼滤波器 (Kalman Filter) 的预测误差分解 (Prediction Error Decomposition) 1。其基本思想是将联合似然函数 $p(y_0, y_1, \dots, y_T)$ 分解为一系列一步预测概率的乘积： + +$$ +p(y_{} | k, \Theta_c^k) = p(y_0 | k, \Theta_c^k) \prod_{t=1}^{T} p(y_t | y_{[0:t-1]}, k, \Theta_c^k) +$$ + +卡尔曼滤波器是一个递归算法，在每个时间步 $t$，它会根据到 $t-1$ 时刻为止的所有观测信息，给出对当前状态 $x_t$ 的预测分布（预测步），然后利用当前的观测 $y_t$ 来修正这个预测，得到更新后的状态分布（更新步）。在这个过程中，它自然地计算出了一步预测分布 $p(y_t | y_{[0:t-1]}, \dots)$。在线性高斯模型中，这个分布也是一个高斯分布，其均值和协方差可以由滤波器的中间量解析得到。 + +因此，整个似然函数的对数（log-likelihood）可以被计算为所有一步预测对数似然的和。这个算法的计算复杂度与时间序列的长度 **$T$** 呈线性关系，与模型阶数 **$k$** 呈多项式关系（通常是 **$O(k^3)$**）。这使得在MCMC的每次迭代中评估似然函数在计算上是可行的。 + +### 4.3 RJMCMC输出分析 + +在运行了足够长的RJMCMC链并舍弃了初始的“燃烧期”（burn-in）样本后，我们会得到一系列来自目标后验分布的样本 **$\{(k^{(i)}, \Theta_c^{k,(i)})\}_{i=1}^N$**。对这些样本的分析可以为我们提供关于模型阶数和参数的丰富信息。 + +* **模型阶数的后验分布** : 这是RJMCMC最直接和最重要的输出之一。模型阶数 **$k$** 的后验分布 **$p(k|y)$** 可以通过简单地统计样本中每个阶数出现的频率来估计 ^19^。例如，阶数 **$k=k^*$** 的后验概率可以近似为： + +$$ +$$\hat{P}(k=k^* | y) = \frac{1}{N} \sum_{i=1}^N \mathbb{I}(k^{(i)} = k^*) +$$ + + 其中 **$\mathbb{I}(\cdot)$** 是指示函数。通过绘制这个频率的直方图，我们可以直观地看到数据支持哪些模型阶数。通常，后验分布会集中在一个或少数几个阶数上。 + +* **参数估计与模型平均** : 与传统的先选择一个“最佳”模型再进行参数估计的两步法不同，贝叶斯方法允许我们进行 **模型平均 (Model Averaging)** 。对于任何我们感兴趣的量 **$Q$**（例如，系统的脉冲响应、特定频率的增益等），它的后验期望可以通过对所有RJMCMC样本求平均来估计： + +$$ +$$E[Q|y] \approx \frac{1}{N} \sum_{i=1}^N Q(\Theta_c^{k,(i)}) +$$ + + 这个计算自动地根据每个模型的后验概率对其预测进行了加权，从而考虑了模型不确定性，通常能提供比单一模型更稳健的估计和预测。当然，我们也可以只分析后验概率最高的那个模型（即最大后验概率模型，MAP model）内部的参数分布。 + +* **收敛性诊断** : RJMCMC链的收敛性诊断比固定维度MCMC更具挑战性。除了检查每个固定阶数模型内部参数的轨迹图（trace plots）和自相关函数外，我们还需要监控模型之间的跳转情况 ^18^。理想情况下，链应该能够自由地在所有具有显著后验概率的模型之间频繁跳转。如果链长时间“卡”在某个阶数的模型中，可能意味着提议分布调整不当，导致模型间的接受率过低。 + +### 4.4 结论与实践建议 + +本报告详细阐述了如何利用可逆跳转MCMC（RJMCMC）方法，将规范贝叶斯系统辨识框架从已知模型阶数推广到未知模型阶数。通过将模型阶数的变化与系统特征多项式根的“诞生”与“死亡”联系起来，我们设计了一套具体的、数学上严谨的跨维度跳转策略。核心技术贡献在于为这些跳转推导了精确的雅可比行列式，从而确保了算法的理论正确性。该方法为LTI系统的辨识提供了一个强大而原则性的全贝叶斯解决方案，能够同时、联合地推断模型结构（阶数）和参数。 + +**实践建议 (Practical Recommendations):** + +1. **调整与效率** : 跨模型跳转的接受率对算法性能至关重要，但它可能非常低。辅助变量的提议分布 **$q(u)$** 的方差是一个关键的调整参数。方差太小，提议的改动不大，可能容易被接受，但探索缓慢；方差太大，提议可能跳到后验概率很低的区域，导致接受率极低。建议通过一些初步的试运行来调整这些参数，目标是使跨模型跳转的接受率达到一个合理的水平（例如，5%-20%）。 +2. **高级提议策略** : 为了提高效率，可以考虑使用更智能的提议分布，而不是简单的先验或标准分布。一种高级策略是 **数据驱动的或自适应的提议** 。例如，在提议一个“诞生”跳转时，可以首先分析当前模型 **$M_k$** 的残差序列。残差的谱分析可能会揭示出当前模型未能捕捉到的某些频率成分。然后，可以设计提议分布 **$q(u)$**，使其倾向于生成能够解释这些残差动态的新特征值（例如，位于谱峰值对应频率的特征值）。这种方法可以显著提高提议的质量和接受率。 +3. **实现细节** : + +* **数值库** : 算法的实现需要依赖高质量的数值计算库。特别是，在“死亡”跳转中需要进行多项式求根，这需要一个稳定可靠的求根算法。现代科学计算库（如Python中的NumPy/SciPy，或MATLAB）都提供了这样的功能。 +* **代码结构** : 建议将每种跳转类型（内部更新、诞生实根、诞生复根对等）封装为独立的函数，每个函数都返回接受概率和新状态。主循环则根据随机选择调用这些函数。 +* **自动化工具** : 值得注意的是，一些现代概率编程语言，如NIMBLE ^42^ 和Gen ^14^，正在开发用于自动化RJMCMC的工具。然而，对于本报告中涉及的这种高度定制化的、基于特定数学变换的跳转，很可能需要从头开始编写自定义的采样器。本报告提供的详细推导正是为此目的服务的。 + +总之，将RJMCMC与规范贝叶斯系统辨识相结合，为解决实际工程中模型阶数未知的挑战提供了一条坚实的道路。尽管实现上具有挑战性，但其所能提供的关于模型结构和参数不确定性的完整图景，是传统方法难以比拟的。 diff --git a/demo.py b/demo.py index 96672a5..1e99bfb 100644 --- a/demo.py +++ b/demo.py @@ -1,99 +1,312 @@ import numpy as np import matplotlib.pyplot as plt -from scipy.stats import beta as beta_dist # 导入beta分布用于绘图 +import scipy.stats as stats +from scipy.special import gammaln # 用于计算 log(Γ(x)) -# ------------------------------------------------------------------ -# 1. 设置 matplotlib 支持中文显示 -# ------------------------------------------------------------------ -try: - plt.rcParams['font.sans-serif'] = ['SimHei'] # Windows/Linux - plt.rcParams['axes.unicode_minus'] = False # 正常显示负号 -except Exception: - try: - plt.rcParams['font.sans-serif'] = ['Arial Unicode MS'] # MacOS - plt.rcParams['axes.unicode_minus'] = False - except Exception: - print("未找到中文字体，绘图可能显示异常。请安装'SimHei'或'Arial Unicode MS'字体。") +# --- 1. 定义先验和似然函数 --- -# ------------------------------------------------------------------ -# 2. 设定模型参数 (为了模拟您图中的效果) -# ------------------------------------------------------------------ -n = 20 # X的试验总次数 -# 更改 alpha 和 beta 以匹配图中 ~0.72 的均值 -alpha = 8.0 # Beta分布的先验参数 alpha -beta = 3.0 # Beta分布的先验参数 beta -# 理论均值 E[Y] = alpha / (alpha + beta) = 8 / 11 ≈ 0.727 -theoretical_mean = alpha / (alpha + beta) +# 定义先验参数 +# π(λ) ~ Gamma(α, β) (注意：scipy.stats.gamma用 a=shape, scale=1/rate) +# 我们使用 α=2, β=1 (rate=1) 作为先验 +ALPHA_LAM = 2 +BETA_LAM = 1 # 这是 rate (或 1/scale) -# MCMC (Gibbs) 抽样参数 -N_samples = 1000 # 总抽样量 N (同您图中的 N=1000) -N_burn_in = 200 # 预估的老化期(预热期) N1 +# π(r) ~ InverseGamma(α, β) (注意：scipy.stats.invgamma用 a=shape, scale=scale) +# 我们使用 α=2, β=1 (scale=1) 作为先验 +ALPHA_R = 2 +BETA_R = 1 # 这是 scale -print(f"模型参数: n={n}, alpha={alpha}, beta={beta}") -print(f"理论均值 E[Y]: {theoretical_mean:.4f}") -print(f"抽样设置: 总样本 N={N_samples}") +# 模型先验 +LOG_PRIOR_K1 = np.log(0.5) +LOG_PRIOR_K2 = np.log(0.5) -# ------------------------------------------------------------------ -# 3. 初始化 -# ------------------------------------------------------------------ -# 创建数组来存储所有样本 -samples_X = np.zeros(N_samples, dtype=int) -samples_Y = np.zeros(N_samples, dtype=float) +# 模型跳跃提议概率 q(k'|k) +# q(1|1)=0.5, q(2|1)=0.5, q(1|2)=0.5, q(2|2)=0.5 +LOG_Q_1_GIVEN_1 = np.log(0.5) +LOG_Q_2_GIVEN_1 = np.log(0.5) +LOG_Q_1_GIVEN_2 = np.log(0.5) +LOG_Q_2_GIVEN_2 = np.log(0.5) -# 设定马尔可夫链的初始状态 -# 故意设置一个远离均值(0.727)的初始值，以观察收敛 -y_t = 0.1 -# ------------------------------------------------------------------ -# 4. 运行 Gibbs 抽样 -# ------------------------------------------------------------------ -print("开始Gibbs抽样...") -np.random.seed(101) # 使用和您图中一样的随机种子 - -for i in range(N_samples): - x_t = np.random.binomial(n, y_t) - y_t = np.random.beta(x_t + alpha, n - x_t + beta) - samples_Y[i] = y_t - samples_X[i] = x_t -print("抽样完成。") - -# ------------------------------------------------------------------ -# 5. 定义并计算逐步平均值 (Ergodic Mean) -# ------------------------------------------------------------------ -def calculate_ergodic_mean(y_samples): - """ - 计算逐步平均值 (累积平均值) - y_k_bar = (1/k) * sum(y_i for i=1 to k) +def get_log_prior(k, params): + """计算参数的对数先验概率""" + if k == 1: + lam = params[0] + if lam <= 0: + return -np.inf + # π(λ) + return stats.gamma.logpdf(lam, a=ALPHA_LAM, scale=1.0/BETA_LAM) - 使用 np.cumsum() 可以高效实现 - """ - n = len(y_samples) - # 1. 计算累积和 [y1, y1+y2, y1+y2+y3, ...] - s = np.cumsum(y_samples) - # 2. 创建 k 数组 [1, 2, 3, ...] - k_array = np.arange(1, n + 1) - # 3. 计算 avg[k] = s[k] / k - return s / k_array + elif k == 2: + lam, r = params + if lam <= 0 or r <= 0: + return -np.inf + # π(λ, r) = π(λ) * π(r) (假设先验独立) + log_p_lam = stats.gamma.logpdf(lam, a=ALPHA_LAM, scale=1.0/BETA_LAM) + log_p_r = stats.invgamma.logpdf(r, a=ALPHA_R, scale=BETA_R) + return log_p_lam + log_p_r -# 计算所有 Y 样本的逐步平均值 (包括预热期) -ergodic_mean_Y = calculate_ergodic_mean(samples_Y) +def get_log_likelihood(k, params, data): + """计算数据的对数似然""" + if k == 1: + lam = params[0] + if lam <= 0: + return -np.inf + # Model 1: Poisson(λ) + return stats.poisson.logpmf(data, lam).sum() + + elif k == 2: + lam, r = params + if lam <= 0 or r <= 0: + return -np.inf + # Model 2: Negative Binomial + # 使用 p = r / (λ + r) 的参数化 + p = r / (lam + r) + # 必须检查 p 是否在 (0, 1] 范围内 + if p <= 0 or p > 1: + return -np.inf + return stats.nbinom.logpmf(data, n=r, p=p).sum() -# ------------------------------------------------------------------ -# 6. 绘制逐步平均值图 (实现您图片中的效果) -# ------------------------------------------------------------------ -print("正在绘制逐步平均值图...") +def get_log_posterior(k, params, data): + """计算完整的对数后验（正比于）""" + log_prior = get_log_prior(k, params) + if log_prior == -np.inf: + return -np.inf + + log_lik = get_log_likelihood(k, params, data) + if log_lik == -np.inf: + return -np.inf + + log_model_prior = LOG_PRIOR_K1 if k == 1 else LOG_PRIOR_K2 + + return log_lik + log_prior + log_model_prior -plt.figure(figsize=(10, 6)) -plt.plot(np.arange(1, N_samples + 1), ergodic_mean_Y, label=r"逐步平均值 $\bar{y}_k$") -plt.axhline(theoretical_mean, color='red', linestyle='--', label=f"理论均值: {theoretical_mean:.4f}") +# --- 2. 生成模拟数据 --- -# 添加一个垂直线来标记我们估计的预热期 -plt.axvline(N_burn_in, color='gray', linestyle=':', label=f"估计的预热期 N1 = {N_burn_in}") +# 我们故意从一个过度离散的负二项分布生成数据 +# 泊松分布：均值=方差。 负二项：方差 > 均值。 +TRUE_LAMBDA = 5.0 +TRUE_R = 10.0 # R 值变大，方差接近均值 (方差 = 5 + 25/10 = 7.5) +TRUE_P = TRUE_R / (TRUE_LAMBDA + TRUE_R) +np.random.seed(42) +N_data = 5000 # <--- 之前缺失的行 +data = stats.nbinom.rvs(n=TRUE_R, p=TRUE_P, size=N_data) -plt.title("使用逐步平均值图查看预热期") -plt.xlabel("迭代次数 (k)") -plt.ylabel(r"逐步平均值 $\bar{y}_k$") -plt.legend() -plt.grid(True) -plt.ylim(0, 1) # Y值在0到1之间 +print(f"模拟数据均值: {data.mean():.2f} (真实均值 = {TRUE_LAMBDA})") +print(f"模拟数据方差: {data.var():.2f} (泊松模型的方差应为 {data.mean():.2f})") + +# --- 3. RJMCMC 主函数 --- + +def run_rjmcmc(data, n_iter=50000, burn_in=10000): + + # 初始化 + # 从模型1开始，λ 使用数据的均值 + current_k = 1 + current_lambda = data.mean() + current_params = [current_lambda] + + # 存储轨迹 + trace_k = np.zeros(n_iter, dtype=int) + trace_lambda = np.zeros(n_iter) + trace_r = np.full(n_iter, np.nan) # 仅当 k=2 时有值 + + # 接受计数器 + acceptance = { + "1_to_1": 0, "2_to_2": 0, "1_to_2": 0, "2_to_1": 0 + } + attempts = { + "1_to_1": 0, "2_to_2": 0, "1_to_2": 0, "2_to_1": 0 + } + + for i in range(n_iter): + # 1. 提议一个目标模型 k_prop + # 无论当前 k 是多少，都以 50/50 的概率提议 k=1 或 k=2 + k_prop = np.random.choice([1, 2]) + + # 获取当前的对数后验 + current_log_post = get_log_posterior(current_k, current_params, data) + + # ---------------------------------- + # 情况 A: 模型内移动 (k_prop == current_k) + # ---------------------------------- + if k_prop == current_k: + if current_k == 1: + # --- Model 1 -> Model 1 --- + attempts["1_to_1"] += 1 + + # 提议一个新的 λ (使用正态分布随机游走) + lambda_prop = current_params[0] + np.random.normal(0, 0.5) + prop_params = [lambda_prop] + + # 计算接受率 + prop_log_post = get_log_posterior(1, prop_params, data) + log_alpha = prop_log_post - current_log_post + # (提议分布是对称的, q(λ'|λ) = q(λ|λ')) + + if np.log(np.random.rand()) < log_alpha: + current_params = prop_params + acceptance["1_to_1"] += 1 + + elif current_k == 2: + # --- Model 2 -> Model 2 --- + attempts["2_to_2"] += 1 + + # 提议新的 (λ, r) + lambda_prop = current_params[0] + np.random.normal(0, 0.5) + r_prop = current_params[1] + np.random.normal(0, 0.5) + prop_params = [lambda_prop, r_prop] + + # 计算接受率 + prop_log_post = get_log_posterior(2, prop_params, data) + log_alpha = prop_log_post - current_log_post + + if np.log(np.random.rand()) < log_alpha: + current_params = prop_params + acceptance["2_to_2"] += 1 + + # ---------------------------------- + # 情况 B: 跨模型移动 (k_prop != current_k) + # ---------------------------------- + else: + if current_k == 1 and k_prop == 2: + # --- Model 1 -> Model 2 (诞生) --- + attempts["1_to_2"] += 1 + + # 1. 抽取辅助变量 w + w = np.random.uniform(0, 1) + log_g_w = stats.uniform.logpdf(w, 0, 1) # 这是 log(1) = 0 + + # 2. 应用映射 + lambda_prop = current_params[0] + r_prop = -np.log(w) + prop_params = [lambda_prop, r_prop] + + # 3. 计算雅可比项 |J| = 1/w + log_jacobian = np.log(1.0 / w) + + # 4. 计算接受率 + prop_log_post = get_log_posterior(2, prop_params, data) + + # log_alpha = (log_post_prop + log_q_backward) - (log_post_curr + log_q_forward + log_g_w) + log_jacobian + log_alpha = (prop_log_post + LOG_Q_1_GIVEN_2) - \ + (current_log_post + LOG_Q_2_GIVEN_1 + log_g_w) + \ + log_jacobian + + if np.log(np.random.rand()) < log_alpha: + current_k = 2 + current_params = prop_params + acceptance["1_to_2"] += 1 + + elif current_k == 2 and k_prop == 1: + # --- Model 2 -> Model 1 (死亡) --- + attempts["2_to_1"] += 1 + + # 1. 这是一个确定性映射（h' 的逆） + current_lambda, current_r = current_params + + # 2. 应用逆映射 + lambda_prop = current_lambda + w_prime = np.exp(-current_r) # 这就是辅助变量 w' + prop_params = [lambda_prop] + + # 3. 计算雅可比项 |J'| = e^(-r) + log_jacobian_prime = np.log(np.exp(-current_r)) # 即 -current_r + + # 4. 计算 g(w')，这是 *正向* 移动 (1->2) 中 w 的密度 + # 正向移动是 w ~ U(0, 1)，所以 g(w') = 1 (只要 0 < w' < 1) + # 因为 r > 0, 所以 w' = e^(-r) 总是在 (0, 1) 区间内 + log_g_w_prime = stats.uniform.logpdf(w_prime, 0, 1) # log(1) = 0 + + # 5. 计算接受率 + prop_log_post = get_log_posterior(1, prop_params, data) + + # log_alpha = (log_post_prop + log_q_backward + log_g_w_prime) - (log_post_curr + log_q_forward) + log_jacobian + log_alpha = (prop_log_post + LOG_Q_2_GIVEN_1 + log_g_w_prime) - \ + (current_log_post + LOG_Q_1_GIVEN_2) + \ + log_jacobian_prime + + if np.log(np.random.rand()) < log_alpha: + current_k = 1 + current_params = prop_params + acceptance["2_to_1"] += 1 + + # 存储当前状态 + trace_k[i] = current_k + trace_lambda[i] = current_params[0] + if current_k == 2: + trace_r[i] = current_params[1] + else: + trace_r[i] = np.nan + + # 打印接受率 + print("\n--- 接受率 ---") + for move, count in attempts.items(): + if count > 0: + rate = acceptance[move] / count + print(f"{move}: {acceptance[move]}/{count} ({rate:.2%})") + + # 丢弃 Burn-in + trace_k_burned = trace_k[burn_in:] + trace_lambda_burned = trace_lambda[burn_in:] + trace_r_burned = trace_r[burn_in:] + + return trace_k_burned, trace_lambda_burned, trace_r_burned + +# --- 4. 运行和绘图 --- + +N_ITER = 20000 +BURN_IN = 5000 +trace_k, trace_lambda, trace_r = run_rjmcmc(data, n_iter=N_ITER, burn_in=BURN_IN) + +# --- 绘图 --- +plt.rcParams['font.sans-serif'] = ['SimHei'] # 用来正常显示中文标签 +plt.rcParams['axes.unicode_minus'] = False # 用来正常显示负号 + +# 图 1: 模型后验概率 +plt.figure(figsize=(12, 10)) + +prob_k1 = np.mean(trace_k == 1) +prob_k2 = np.mean(trace_k == 2) + +ax1 = plt.subplot(3, 1, 1) +bars = plt.bar([1, 2], [prob_k1, prob_k2], color=["#227dbe", "#f97b0df9"], tick_label=['模型 1 (泊松)', '模型 2 (负二项)']) +plt.title('模型后验概率', fontsize=16) +plt.ylabel('P(k | data)', fontsize=12) +ax1.bar_label(bars, fmt='{:.2%}', fontsize=12) +ax1.set_ylim(0, 1) + +# 图 2: 模型跳跃轨迹 (仅显示前 2000 步，看得更清楚) +ax2 = plt.subplot(3, 1, 2) +ax2.plot(trace_k[:2000], 'k.', markersize=2, alpha=0.5) +ax2.set_yticks([1, 2]) +ax2.set_yticklabels(['模型 1 (泊松)', '模型 2 (负二项)']) +ax2.set_title('模型空间轨迹 (前2000次迭代)', fontsize=16) +ax2.set_xlabel('迭代次数', fontsize=12) + +# 图 3: 参数后验分布 +# λ (lambda) 的后验 +ax3 = plt.subplot(3, 2, 5) +ax3.hist(trace_lambda, bins=50, density=True, color='#1f77b4', alpha=0.7, label='$\lambda$ 的后验分布') +ax3.axvline(data.mean(), color='red', linestyle='--', label=f'数据均值 ({data.mean():.2f})') +ax3.axvline(TRUE_LAMBDA, color='black', linestyle=':', label=f'真实 $\lambda$ ({TRUE_LAMBDA})') +ax3.set_title('参数 $\lambda$ 的后验分布', fontsize=14) +ax3.set_xlabel('$\lambda$ 值', fontsize=12) +ax3.set_ylabel('密度', fontsize=12) +ax3.legend() + +# r 的后验 +ax4 = plt.subplot(3, 2, 6) +# 仅使用 k=2 时的 r 值 +trace_r_k2 = trace_r[~np.isnan(trace_r)] +if len(trace_r_k2) > 0: + ax4.hist(trace_r_k2, bins=50, density=True, color='#ff7f0e', alpha=0.7, label='$r$ 的后验分布 (当 k=2)') + ax4.axvline(TRUE_R, color='black', linestyle=':', label=f'真实 $r$ ({TRUE_R})') + ax4.set_title('参数 $r$ 的后验分布 (仅当k=2)', fontsize=14) + ax4.set_xlabel('$r$ 值', fontsize=12) + ax4.legend() +else: + ax4.set_title('参数 $r$ 的后验分布 (未采样到)', fontsize=14) + ax4.text(0.5, 0.5, '从未接受过模型 2', horizontalalignment='center', verticalalignment='center', transform=ax4.transAxes) + +plt.tight_layout() plt.show() diff --git a/demo1.py b/demo1.py index a6d2f3d..55351e4 100644 --- a/demo1.py +++ b/demo1.py @@ -1,127 +1,173 @@ -# JAX 入门示例代码 -# 演示 JAX 的基本用法，包括 jnp 数组、jit 加速 - import numpy as np +import numpyro import jax import jax.numpy as jnp -from jax import jit, grad, vmap -import time +import numpyro.distributions as dist +from numpyro.infer import MCMC, NUTS +import matplotlib.pyplot as plt +from scipy import stats +from scipy.stats import gaussian_kde +from functools import partial # 引入 partial 来固定函数参数 -# --- 1. 像 NumPy 一样简单 --- -print("--- 1. NumPy-like API ---") -key = jax.random.PRNGKey(0) # JAX 处理随机数需要一个密钥 -x_jnp = jnp.arange(10) -y_jnp = jax.random.normal(key, (10,)) +# ========================================================================= +# 配置中文字体 +# ========================================================================= +try: + # Windows 系统优先尝试这些字体 + plt.rcParams['font.sans-serif'] = ['Microsoft YaHei', 'SimHei', 'SimSun', 'KaiTi', 'FangSong', 'Arial Unicode MS'] + plt.rcParams['axes.unicode_minus'] = False # 正常显示负号 + print("✓ 中文字体配置成功") +except Exception as e: + print(f"⚠ 字体配置警告: {e}") + print(" 如果图表中文显示异常，请运行 check_fonts.py 查看可用字体") -print(f"JAX 数组 x: {x_jnp}") -print(f"JAX 数组 y (随机): {y_jnp}") -print(f"x 和 y 的点积: {jnp.dot(x_jnp, y_jnp)}") +# 设置随机种子 +numpyro.set_platform("cpu") +numpyro.set_host_device_count(4) +# ========================================================================= +# 1. 参数化的 MCMC 模型 +# ========================================================================= -# --- 2. jax.jit: 极致加速 --- -print("\n--- 2. @jit 加速 ---") -# 定义一个包含大量运算的函数 -def slow_function(x): - # 这是一个比较耗时的操作 (矩阵乘法循环) - for _ in range(100): - x = jnp.dot(x, x.T) - return x +def simple_model_factor(y_obs, mu_loc, mu_scale, sigma_a, sigma_b): + """ + 参数化的正态分布模型 + 先验参数作为函数参数传入 + """ + # 先验分布：mu ~ Normal(mu_loc, mu_scale) + mu = numpyro.sample("mu", dist.Normal(mu_loc, mu_scale)) + # 先验分布：sigma ~ Beta(sigma_a, sigma_b) + sigma = numpyro.sample("sigma", dist.Beta(sigma_a, sigma_b)) + + # 使用 numpyro.factor 添加对数似然 + n = len(y_obs) + log_lik = -0.5 * n * jnp.log(2 * jnp.pi) - n * jnp.log(sigma) - jnp.sum((y_obs - mu) ** 2) / (2 * sigma ** 2) + numpyro.factor("log_likelihood", log_lik) -# 创建 JIT 编译版本的函数 -fast_function = jit(slow_function) +# ========================================================================= +# 2. 生成模拟数据 +# ========================================================================= +np.random.seed(42) +true_mu = 2.5 # 真实均值 +true_sigma = 0.2 # 真实标准差 +n_data = 1000 +y_data = np.random.normal(true_mu, true_sigma, n_data) +y_data_jax = jnp.array(y_data) -# 创建数据 -big_matrix = jax.random.normal(key, (200, 200)) +sample_mean = np.mean(y_data) +sample_std = np.std(y_data) -# --- 计时对比 --- -# a) 运行普通 NumPy/Python 版本 -start_time = time.time() -slow_function(big_matrix).block_until_ready() # .block_until_ready() 确保 JAX 计算完成 -numpy_time = time.time() - start_time -print(f"普通 NumPy/Python 版本耗时: {numpy_time:.6f} 秒") +print(f"真实均值: {true_mu}, 真实标准差: {true_sigma}") +print(f"样本均值: {sample_mean:.3f}, 样本标准差: {sample_std:.3f}") -# b) 运行 JIT 编译版本 -# 第一次运行会包含编译时间 -start_time = time.time() -fast_function(big_matrix).block_until_ready() -compile_time = time.time() - start_time -print(f"JIT 版本 (首次，含编译) 耗时: {compile_time:.6f} 秒") +# ========================================================================= +# 3. 运行两次 MCMC +# ========================================================================= -# 第二次运行，将只体现执行速度 -start_time = time.time() -fast_function(big_matrix).block_until_ready() -jit_time = time.time() - start_time -print(f"JIT 版本 (第二次) 耗时: {jit_time:.6f} 秒") +# --- 运行 1: "合理"先验 (Reasonable Prior) --- +# mu ~ Normal(0, 10), sigma ~ Beta(2, 5) [均值 approx 0.28] +print("\n--- 正在运行 MCMC (合理先验) ---") +nuts_kernel_1 = NUTS( + partial(simple_model_factor, mu_loc=0., mu_scale=10., sigma_a=2., sigma_b=5.) +) +mcmc_1 = MCMC(nuts_kernel_1, num_warmup=100, num_samples=300, num_chains=4) +mcmc_1.run(jax.random.PRNGKey(0), y_obs=y_data_jax) +mcmc_samples_1 = mcmc_1.get_samples() +print("✓ MCMC (合理先验) 运行完毕") +# mcmc_1.print_summary() -print(f"加速比 (第二次运行 vs 普通版本): {numpy_time / jit_time:.2f} 倍") +# --- 运行 2: "错误/远处"先验 (Far Prior) --- +# mu ~ Normal(10, 1), sigma ~ Beta(5, 2) [均值 approx 0.71] +print("\n--- 正在运行 MCMC (错误先验) ---") +nuts_kernel_2 = NUTS( + partial(simple_model_factor, mu_loc=10., mu_scale=1., sigma_a=5., sigma_b=2.) +) +mcmc_2 = MCMC(nuts_kernel_2, num_warmup=1000, num_samples=3000, num_chains=4) +mcmc_2.run(jax.random.PRNGKey(1), y_obs=y_data_jax) +mcmc_samples_2 = mcmc_2.get_samples() +print("✓ MCMC (错误先验) 运行完毕") +# mcmc_2.print_summary() +# ========================================================================= +# 4. 对比绘图 (修改版：使用 KDE 曲线避免遮挡) +# ========================================================================= -# --- 3. jax.grad: 自动求导 --- -print("\n--- 3. grad 自动求导 ---") -# 定义一个简单的函数 f(x) = x^3 + 2x^2 + 5 -def my_func(x): - return x**3 + 2*x**2 + 5 +# --- a. 创建用于绘图的网格 --- +mu_plot_grid = np.linspace(-5, 15, 400) +sigma_plot_grid = np.linspace(0.01, 1.0, 400) -# 使用 grad 创建一个计算 my_func 导数的函数 -# f'(x) = 3x^2 + 4x -grad_my_func = grad(my_func) +# --- b. 计算先验的 PDF (不变) --- +# 合理先验 +prior_mu_1_pdf = stats.norm(0, 10).pdf(mu_plot_grid) # type: ignore +prior_sigma_1_pdf = stats.beta(2, 5).pdf(sigma_plot_grid) # type: ignore +# 错误先验 +prior_mu_2_pdf = stats.norm(10, 1).pdf(mu_plot_grid) # type: ignore +prior_sigma_2_pdf = stats.beta(5, 2).pdf(sigma_plot_grid) # type: ignore -x_val = 2.0 -derivative = grad_my_func(x_val) -expected_derivative = 3 * x_val**2 + 4 * x_val +# --- c. 【新】计算后验的 KDE (核密度估计) --- +# 这会根据 MCMC 样本生成平滑的概率密度函数 +print("\n--- 正在计算 KDE (平滑曲线) ---") +kde_mu_1 = gaussian_kde(mcmc_samples_1['mu']) +kde_mu_2 = gaussian_kde(mcmc_samples_2['mu']) +kde_sigma_1 = gaussian_kde(mcmc_samples_1['sigma']) +kde_sigma_2 = gaussian_kde(mcmc_samples_2['sigma']) -print(f"函数 f(x) = x^3 + 2x^2 + 5") -print(f"在 x = {x_val} 处的导数是: {derivative}") -print(f"理论上的导数值是: {expected_derivative}") +# 在网格上计算 KDE 的 PDF 值 +post_mu_1_pdf = kde_mu_1(mu_plot_grid) +post_mu_2_pdf = kde_mu_2(mu_plot_grid) +post_sigma_1_pdf = kde_sigma_1(sigma_plot_grid) +post_sigma_2_pdf = kde_sigma_2(sigma_plot_grid) +print("✓ KDE 计算完毕") +# --- d. 开始绘图 --- +plt.figure(figsize=(16, 8)) -# --- 4. jax.vmap: 自动向量化 --- -print("\n--- 4. vmap 自动向量化 ---") -# 定义一个只能处理单个向量的函数 (向量点积) -def single_dot_product(a, b): - return jnp.dot(a, b) +# --- 图 1: 对比 mu 的分布 --- +ax1 = plt.subplot(1, 2, 1) +# 绘制先验 (虚线, 稍透明) +ax1.plot(mu_plot_grid, prior_mu_1_pdf, 'b--', label='先验 1 (合理): N(0, 10)', linewidth=2, alpha=0.7) +ax1.plot(mu_plot_grid, prior_mu_2_pdf, 'r--', label='先验 2 (错误): N(10, 1)', linewidth=2, alpha=0.7) -# 创建一批 (batch) 数据 -# 假设我们有 5 对向量，每对向量长度为 3 -batch_a = jnp.arange(15).reshape(5, 3) -batch_b = jnp.arange(15, 30).reshape(5, 3) +# 绘制后验 (实线/点划线，不透明) +# 【修改点】用 plot 代替 hist +ax1.plot(mu_plot_grid, post_mu_1_pdf, 'b-', label='后验 1 (来自合理先验)', linewidth=3) +ax1.plot(mu_plot_grid, post_mu_2_pdf, 'r-.', label='后验 2 (来自错误先验)', linewidth=3) # 使用不同线型 -# 使用 vmap 将函数向量化 -# in_axes=(0, 0) 表示对 a 和 b 的第 0 轴 (批次轴) 进行映射 -# out_axes=0 表示输出结果也沿着第 0 轴堆叠 -batch_dot_product = vmap(single_dot_product, in_axes=(0, 0), out_axes=0) +# 绘制真实值和样本均值 +ax1.axvline(true_mu, color='k', linestyle=':', linewidth=2.5, label=f'真实均值 = {true_mu}') +ax1.axvline(sample_mean, color='gray', linestyle='-', linewidth=2, label=f'样本均值 = {sample_mean:.3f}') # type: ignore -results = batch_dot_product(batch_a, batch_b) +ax1.set_title(r"参数 $\mu$ 的先验与后验对比", fontsize=16) +ax1.set_xlabel(r"$\mu$ 的值") +ax1.set_ylabel("概率密度") +ax1.legend(fontsize=10) +ax1.grid(True, linestyle='--', alpha=0.6) +ax1.set_xlim(-5, 15) +ax1.set_ylim(bottom=0) # 确保 y 轴从 0 开始 -print("批次数据 a:\n", batch_a) -print("批次数据 b:\n", batch_b) -print("vmap 向量化计算的点积结果:\n", results) +# --- 图 2: 对比 sigma 的分布 --- +ax2 = plt.subplot(1, 2, 2) +# 绘制先验 (虚线, 稍透明) +ax2.plot(sigma_plot_grid, prior_sigma_1_pdf, 'b--', label='先验 1 (合理): Beta(2, 5)', linewidth=2, alpha=0.7) +ax2.plot(sigma_plot_grid, prior_sigma_2_pdf, 'r--', label='先验 2 (错误): Beta(5, 2)', linewidth=2, alpha=0.7) -# 对比手动 for 循环 -manual_results = jnp.array([single_dot_product(a, b) for a, b in zip(batch_a, batch_b)]) -print("手动 for 循环的结果:\n", manual_results) -print(f"vmap 结果与手动循环结果是否一致: {jnp.allclose(results, manual_results)}") +# 绘制后验 (实线/点划线，不透明) +# 【修改点】用 plot 代替 hist +ax2.plot(sigma_plot_grid, post_sigma_1_pdf, 'b-', label='后验 1 (来自合理先验)', linewidth=3) +ax2.plot(sigma_plot_grid, post_sigma_2_pdf, 'r-.', label='后验 2 (来自错误先验)', linewidth=3) # 使用不同线型 +# 绘制真实值和样本标准差 +ax2.axvline(true_sigma, color='k', linestyle=':', linewidth=2.5, label=f'真实 $\sigma$ = {true_sigma}') +ax2.axvline(sample_std, color='gray', linestyle='-', linewidth=2, label=f'样本 $\sigma$ = {sample_std:.3f}') # type: ignore -# --- 5. JAX 的随机数处理 --- -print("\n--- 5. JAX 的随机数 ---") -# 1. 创建一个初始密钥 -key = jax.random.PRNGKey(42) -print(f"初始密钥: {key}") +ax2.set_title(r"参数 $\sigma$ 的先验与后验对比", fontsize=16) +ax2.set_xlabel(r"$\sigma$ 的值") +ax2.set_ylabel("概率密度") +ax2.legend(fontsize=10) +ax2.grid(True, linestyle='--', alpha=0.6) +ax2.set_xlim(0, 1.0) +ax2.set_ylim(bottom=0) # 确保 y 轴从 0 开始 -# 2. 使用密钥生成随机数 -random_data_1 = jax.random.normal(key, (3,)) -print(f"第一次生成的随机数据: {random_data_1}") - -# 3. 再次使用同一个密钥，会得到完全相同的结果！ -random_data_2 = jax.random.normal(key, (3,)) -print(f"第二次使用相同密钥生成的数据: {random_data_2}") - -# 4. 正确的做法：分割密钥 -key, subkey = jax.random.split(key) # key 更新为新的主密钥, subkey 用于本次操作 -random_data_3 = jax.random.normal(subkey, (3,)) -print(f"\n分割密钥后，第一次生成的随机数据: {random_data_3}") - -key, subkey = jax.random.split(key) # 再次分割 -random_data_4 = jax.random.normal(subkey, (3,)) -print(f"分割密钥后，第二次生成的随机数据: {random_data_4}") +plt.suptitle("先验信念 vs. 强大数据 (N=1000) [KDE平滑曲线]", fontsize=20, y=1.02) +plt.tight_layout() +plt.show() diff --git a/equ.afx b/equ.afx new file mode 100644 index 0000000..f57abc2 Binary files /dev/null and b/equ.afx differ diff --git a/kalmanFilter.py b/kalmanFilter.py new file mode 100644 index 0000000..7fa5a25 --- /dev/null +++ b/kalmanFilter.py @@ -0,0 +1,28 @@ +import numpy as np +import matplotlib.pyplot as plt +from scipy.stats import beta + +# 定义 Beta 分布的参数 +alpha = 2 +beta_param = 2 # 注意这里变量名不能与模块名beta重复，所以用beta_param + +# 生成 x 值，Beta 分布的定义域是 [0, 1] +# 我们生成一系列在 0 到 1 之间的点 +x = np.linspace(0.01, 0.99, 500) # 避免在0和1处logpdf可能趋向负无穷导致绘图问题，稍微避开边界 + +# 计算每个 x 值的 logpdf +log_pdf_values = beta.logpdf(x, alpha, beta_param) + +# 绘图 +plt.figure(figsize=(10, 6)) +plt.plot(x, log_pdf_values, label=f'logPDF of Beta(α={alpha}, β={beta_param})') + +# 添加标题和标签 +plt.title(f'Log-Probability Density Function (logPDF) of Beta(α={alpha}, β={beta_param})') +plt.xlabel('x') +plt.ylabel('log(PDF)') +plt.grid(True) +plt.legend() + +# 显示图形 +plt.show() diff --git a/main.py b/main.py index 13dff69..a05f715 100644 --- a/main.py +++ b/main.py @@ -30,41 +30,70 @@ except Exception as e: def standard_to_canonical(A_s, B_s, C_s): """ 将一个 2x2 的标准状态空间系统 (A_s, B_s, C_s) 转换为控制器规范型参数。 - + + 此实现基于论文附录 D.2 的推导， + 通过计算特征多项式和可控性矩阵来找到变换矩阵 T_c。 + + 参数: + A_s (np.ndarray): 2x2 状态矩阵 + B_s (np.ndarray): 2x1 输入矩阵 + C_s (np.ndarray): 1x2 观测矩阵 + 返回: dict: 包含 'a0', 'a1', 'b0', 'b1' 的字典 """ - dx = A_s.shape[0] - if dx != 2: - raise ValueError("此转换函数仅为 dx=2 的情况实现。") + + # 确保输入是 numpy 数组 + A_s = np.asarray(A_s) + B_s = np.asarray(B_s) + C_s = np.asarray(C_s) + + if A_s.shape != (2, 2) or B_s.shape != (2, 1) or C_s.shape != (1, 2): + raise ValueError(f"输入维度不正确: A_s {A_s.shape}, B_s {B_s.shape}, C_s {C_s.shape}") - # 1. 计算特征多项式系数: p(λ) = λ^2 + a1*λ + a0 + # 1. 计算特征多项式系数 (来自 Ac) + # p(λ) = λ^2 - tr(A_s)λ + det(A_s) + # 规范型 p(λ) = λ^2 + a1*λ + a0 + # 比较系数: a1 = -tr(A_s), a0 = det(A_s) a1_true = -np.trace(A_s) a0_true = np.linalg.det(A_s) - - # 2. 构造转换矩阵 T_c - I = np.eye(dx) - f1 = (A_s + a1_true * I) @ B_s + + # 2. 构造逆转换矩阵 T_c^{-1} = [f1, f2] + I = np.eye(2) + + # 根据附录 D.2 (1157), f_k = (A_s + a_{d-1}I)f_{k+1} + ... f2 = B_s - - Tc_inv = np.hstack([f1, f2]) - - if np.linalg.matrix_rank(Tc_inv) < dx: - raise np.linalg.LinAlgError("系统 (A_s, B_s) 不是可控的，无法转换为控制器规范型。") - - Tc = np.linalg.inv(Tc_inv) - - # 3. 转换 C 矩阵: C_c = C_s * T_c - C_c = C_s @ Tc + f1 = (A_s + a1_true * I) @ B_s + + Tc_inv = np.hstack([f1, f2]) # 这是 T_c + print(f"[standard_to_canonical] 恢复的 T_c (即 Tc_inv):\n{Tc_inv}") + + # 3. 检查可控性 (Controllability) + if np.linalg.matrix_rank(Tc_inv) < 2: + raise ValueError("系统不可控 (Uncontrollable), 无法转换为控制器规范型。") + + # 4. 计算转换矩阵 T_c (这是 T_c^{-1}) + Tc = np.linalg.inv(Tc_inv) # 这是 T_c^{-1} + print(f"[standard_to_canonical] 恢复的 T_c^{{-1}} (即 Tc):\n{Tc}\n") + + # 5. 应用变换找到 C_c = [b0, b1] (来自 Cc) + # 正确的公式是 C_c = C_s * T_c + # 在我们的变量名中, T_c 是 Tc_inv + C_c = C_s @ Tc_inv # 这是正确行 (C_s * T_c) + b0_true = C_c[0, 0] b1_true = C_c[0, 1] - true_params = {'a0': a0_true, 'a1': a1_true, 'b0': b0_true, 'b1': b1_true} - - print("\n--- Ground Truth 转换结果 ---") - print(f"真实规范型参数: {true_params}") - - return true_params + # 打印矩阵 + print(f"[standard_to_canonical] 计算得到的规范型参数:") + print(f" a0: {a0_true}, a1: {a1_true}, b0: {b0_true}, b1: {b1_true}\n") + + return { + 'a0': a0_true, + 'a1': a1_true, + 'b0': b0_true, + 'b1': b1_true + } # ========================================================================= # 可视化函数 @@ -192,7 +221,7 @@ if __name__ == '__main__': # --- 2. 仿真数据 --- print("\n--- (步骤 2) 生成仿真数据 ---") - T_steps, sigma_proc, sigma_meas = 400, 0.3, 0.5 + T_steps, sigma_proc, sigma_meas = 800, 0.05, 0.05 u_data, y_data = simulate_lti_data( A_true, B_true, C_true, D_true, T_steps, sigma_proc, sigma_meas, rng_seed=int(sim_key[0]) @@ -219,9 +248,9 @@ if __name__ == '__main__': print("\n--- 运行规范型模型 MCMC (在真实值附近初始化) ---") mcmc_key_c, mcmc_key = jax.random.split(mcmc_key) - # --- (修改) 在真实值附近添加小的随机扰动 (±10%) --- + # --- (修改) 在真实值附近添加小的随机扰动 (±30%) --- init_key = init_noise_key_c # 使用独立的 key - noise_scale = 0.1 # 10% 的扰动 + noise_scale = 0.3 # 30% 的扰动 init_params_c_noisy = {} for param_name, true_value in true_params_c.items(): # 确保即使 true_value 为 0 也有扰动，添加一个小的基准值 @@ -250,16 +279,17 @@ if __name__ == '__main__': u_data_jax, y_data_jax, sigma_proc, - init_params=init_params_c_noisy # <-- 使用带扰动的初始值 + sigma_meas, + init_params = None # <-- 使用带扰动的初始值 ) if choice in [2, 3]: # 运行标准型 print("\n--- 运行标准型模型 MCMC (在真实值附近初始化) ---") mcmc_key_s, mcmc_key = jax.random.split(mcmc_key) - # --- (修改) 在真实值附近添加小的随机扰动 (±10%) --- + # --- (修改) 在真实值附近添加小的随机扰动 (±30%) --- init_key = init_noise_key_s # 使用独立的 key - noise_scale = 0.1 # 10% 的扰动 + noise_scale = 0.3 # 30% 的扰动 init_params_s_noisy = {} for param_name, true_value in true_params_s.items(): # 确保即使 true_value 为 0 也有扰动 @@ -276,6 +306,7 @@ if __name__ == '__main__': u_data_jax, y_data_jax, sigma_proc, + sigma_meas, init_params=init_params_s_noisy # <-- 使用带扰动的初始值 ) diff --git a/models_and_mcmc.py b/models_and_mcmc.py index 4e82f2f..6a97cac 100644 --- a/models_and_mcmc.py +++ b/models_and_mcmc.py @@ -63,7 +63,11 @@ def kalman_likelihood(dx, du, dy, T, A, B, C, D, Q, R, u_data, y_data): S_t = C @ P_t_tm1 @ C.T + R # 计算对数似然 p(y_t | y_{t-1}, ...) - log_lik_t = multivariate_normal.logpdf(nu_t.squeeze(), mean=jnp.zeros(dy), cov=S_t) + sign, logdet_S_t = jnp.linalg.slogdet(S_t) + S_inv_nu = jnp.linalg.solve(S_t, nu_t) + quad_term = (nu_t.T @ S_inv_nu).squeeze() + log_2pi = jnp.log(2.0 * jnp.pi) + log_lik_t = -0.5 * dy * log_2pi - 0.5 * logdet_S_t - 0.5 * quad_term # 卡尔曼增益 K_t = P_t_tm1*C^T * S_t^{-1} K_t = jnp.linalg.solve(S_t, C @ P_t_tm1).T @@ -77,7 +81,7 @@ def kalman_likelihood(dx, du, dy, T, A, B, C, D, Q, R, u_data, y_data): # --- 2. 时间预测 (Time Prediction) --- # (使用 u_t 来从 x_{t|t} 得到 x_{t+1|t}) - # 预测下一个状态 x_{t+1|t} = A*x_{t|t} + B*u_t [cite: 21, 1042] + # 预测下一个状态 x_{t+1|t} = A*x_{t|t} + B*u_t x_tp1_t = (A @ x_t_t + B @ u_t).flatten() # 确保输出是 (dx,) 形状 # 预测下一个协方差 P_{t+1|t} = A*P_{t|t}*A^T + Q @@ -99,7 +103,7 @@ def kalman_likelihood(dx, du, dy, T, A, B, C, D, Q, R, u_data, y_data): # ========================================================================= # 步骤 3.1: 定义模型一 (规范型, Canonical) # ========================================================================= -def model_canonical(u_data, y_data, sigma_process, nugget=1e-12): +def model_canonical(u_data, y_data, sigma_process, sigma_measure): """ NumPyro 模型 - 规范型 (Canonical Form) """ @@ -116,8 +120,8 @@ def model_canonical(u_data, y_data, sigma_process, nugget=1e-12): a1 = numpyro.sample("a1", dist.Uniform(-1 - a0, 1 + a0)) # type: ignore # 观测矩阵 C 的先验 - b0 = numpyro.sample("b0", dist.Normal(0, 1)) - b1 = numpyro.sample("b1", dist.Normal(0, 1)) + b0 = numpyro.sample("b0", dist.Normal(0, 2)) + b1 = numpyro.sample("b1", dist.Normal(0, 2)) # --- 2. 构造系统矩阵 --- A = jnp.array([[0.0, 1.0], [-a0, -a1]]) # type: ignore @@ -128,8 +132,7 @@ def model_canonical(u_data, y_data, sigma_process, nugget=1e-12): # --- 3. 构造噪声协方差 --- # 噪声是固定的 (已知的)，如 6.3 节算例所述 Q = jnp.eye(dx) * (sigma_process ** 2) - # 为数值稳定性添加 "nugget" [cite: 533] - R = jnp.eye(dy) * nugget + R = jnp.eye(dy) * (sigma_measure ** 2) # --- 4. 计算总似然 --- log_lik_total = kalman_likelihood(dx, du, dy, T, A, B, C, D, Q, R, u_data, y_data) @@ -140,7 +143,7 @@ def model_canonical(u_data, y_data, sigma_process, nugget=1e-12): # ========================================================================= # 步骤 3.2: 定义模型二 (标准型, Standard) # ========================================================================= -def model_standard(u_data, y_data, sigma_process, nugget=1e-12): +def model_standard(u_data, y_data, sigma_process, sigma_measure): """ NumPyro 模型 - 标准型 (Standard Form) """ @@ -152,18 +155,18 @@ def model_standard(u_data, y_data, sigma_process, nugget=1e-12): # 所有系数都是 N(0, 1) # 状态矩阵 A (dx*dx = 4 个参数) - A11 = numpyro.sample("A11", dist.Normal(0, 1)) - A12 = numpyro.sample("A12", dist.Normal(0, 1)) - A21 = numpyro.sample("A21", dist.Normal(0, 1)) - A22 = numpyro.sample("A22", dist.Normal(0, 1)) + A11 = numpyro.sample("A11", dist.Normal(0, 2)) + A12 = numpyro.sample("A12", dist.Normal(0, 2)) + A21 = numpyro.sample("A21", dist.Normal(0, 2)) + A22 = numpyro.sample("A22", dist.Normal(0, 2)) # 输入矩阵 B (dx*du = 2 个参数) - B1 = numpyro.sample("B1", dist.Normal(0, 1)) - B2 = numpyro.sample("B2", dist.Normal(0, 1)) + B1 = numpyro.sample("B1", dist.Normal(0, 2)) + B2 = numpyro.sample("B2", dist.Normal(0, 2)) # 观测矩阵 C (dy*dx = 2 个参数) - C1 = numpyro.sample("C1", dist.Normal(0, 1)) - C2 = numpyro.sample("C2", dist.Normal(0, 1)) + C1 = numpyro.sample("C1", dist.Normal(0, 2)) + C2 = numpyro.sample("C2", dist.Normal(0, 2)) # --- 2. 构造系统矩阵 --- A = jnp.array([[A11, A12], [A21, A22]]) @@ -173,7 +176,7 @@ def model_standard(u_data, y_data, sigma_process, nugget=1e-12): # --- 3. 构造噪声协方差 --- Q = jnp.eye(dx) * (sigma_process ** 2) - R = jnp.eye(dy) * nugget + R = jnp.eye(dy) * (sigma_measure ** 2) # --- 4. 计算总似然 --- log_lik_total = kalman_likelihood(dx, du, dy, T, A, B, C, D, Q, R, u_data, y_data) @@ -184,14 +187,14 @@ def model_standard(u_data, y_data, sigma_process, nugget=1e-12): # ========================================================================= # (步骤 4: 运行 MCMC - 作为本脚本的 main) # ========================================================================= -def run_mcmc(model, rng_key, u_data, y_data, sigma_process, init_params=None): +def run_mcmc(model, rng_key, u_data, y_data, sigma_process, sigma_meas, init_params=None): """辅助函数，用于运行 NUTS 采样器""" print(f"\n--- 开始为模型 {model.__name__} 运行 MCMC ---") # 论文中的 MCMC 设置 - num_warmup = 5000 - num_samples = 20000 + num_warmup = 20000 + num_samples = 40000 num_chains = 4 # 使用 NUTS 内核 @@ -216,7 +219,7 @@ def run_mcmc(model, rng_key, u_data, y_data, sigma_process, init_params=None): ) # 运行 - mcmc.run(rng_key, u_data, y_data, sigma_process=sigma_process) + mcmc.run(rng_key, u_data, y_data, sigma_process=sigma_process, sigma_measure=sigma_meas) # 打印总结 print(f"\n--- MCMC 总结: {model.__name__} ---") @@ -266,7 +269,8 @@ if __name__ == '__main__': mcmc_key_c, u_data_jax, y_data_jax, - sigma_proc + sigma_proc, + sigma_meas ) # (模型 2: 标准型) @@ -275,7 +279,8 @@ if __name__ == '__main__': mcmc_key_s, u_data_jax, y_data_jax, - sigma_proc + sigma_proc, + sigma_meas ) print("\n--- MCMC 运行完成 ---")