From b1ba96b6919c51a693946ce781c9f9f718e24a66 Mon Sep 17 00:00:00 2001 From: Hongru Date: Wed, 5 Nov 2025 08:36:42 +0800 Subject: [PATCH] =?UTF-8?q?=E7=A8=8B=E5=BA=8F=E4=BB=8D=E6=9C=89=E9=97=AE?= =?UTF-8?q?=E9=A2=98?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- README.html | 510 +++++++++++++++++++++++++++++++++ README.md | 36 +-- RJMCMC_ModelSelector.py | 620 ++++++++++++++++++++++++++++++---------- equ.afx | Bin 26112 -> 131072 bytes generateGroudTruth.py | 244 ++++++++++++++-- temp.md | 173 +++++++++++ 6 files changed, 1394 insertions(+), 189 deletions(-) create mode 100644 README.html create mode 100644 temp.md diff --git a/README.html b/README.html new file mode 100644 index 0000000..13f8de3 --- /dev/null +++ b/README.html @@ -0,0 +1,510 @@ + +README.md + + + + + + + + + + + + +
+

MOTIVATION

+

论文中的辨识是认为系统的阶数已知的辨识，而且 MCMC 的抽样的变量是特征多项式矩阵，我想使用特征值来抽样，并同时考虑系统阶数未知的情况。

+

采用RJMMC的方法

+

细致平稳条件

+

α(x,y)=min⁡{1,R}π(x)p(x,y)=π(y)p(y,x)\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⁡P(xk+1)P(y∣xk+1)P(xk)P(y∣xk)+log⁡P(dead)P(birth)+log⁡1k+1q(u)+log⁡∣J1∣\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. +
  2. Death move

    log⁡R=log⁡P(xk)P(y∣xk)P(xk+1)P(y∣xk+1)+log⁡P(birth)P(dead)+log⁡q(u)1k+1−log⁡∣J1∣\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| +

    +
    3. +
+

k <-> k+2

+
    +
  1. Birth move

    log⁡R=log⁡P(xk+2)P(y∣xk+2)P(xk)P(y∣xk)+log⁡P(dead)P(birth)+log⁡1m+1q(θ)q(r)+log⁡∣J2∣\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. +
  2. Death move

    log⁡R=log⁡P(xk)P(y∣xk)P(xk+2)P(y∣xk+2)+log⁡P(birth)P(dead)+log⁡q(θ)q(r)1m+1−log⁡∣J2∣\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| +

    +
    3. +
+

其中m其中 m 是当前系统中复共轭极点对的数量 +J1J_1 和 J2J_2 是雅可比矩阵行列式 +∣J1∣=∣∏i=1k(λi−ui)∣\left| J_1 \right|=\left| \prod_{\boldsymbol{i}=1}^{\boldsymbol{k}}{\left( \lambda _i-u_i \right)} \right| +∣J2∣=∣∏i=1k(λk+1−λi)(λk+2−λi)∣∣λk+1−λk+1∣\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=r1+r22、b=r1−r22r1=a+b、r2=a−ba=\frac{r_1+r_2}{2}\text{、}b=\frac{r_1-r_2}{2} +\\ +r_1=a+b\text{、}r_2=a-b +

+

再确定雅可比矩阵行列式

+

JC→R=∣det⁡(∂a∂r1∂a∂r2∂b∂r1∂b∂r2)∣=12JR→C=JC→R−1=2J_{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. 实根 <-> 复共轭根对 +在nrn_r个实根中任选两个实根r1r_1和r2r_2，转换为复共轭根对a±jba\pm jb，其中a=r1+r22a=\frac{r_1+r_2}{2}，b=r1−r22b=\frac{r_1-r_2}{2}。
    2. +
+

q(x∣x′)=pm×1C(nr,2)q\left( x|x' \right) =p_m\times \frac{1}{C\left( n_r,2 \right)} +

+

此时的接受率写作:Merge：α(x′,x)=min⁡{1,π(x)π(x′)⋅1q(x∣x′)⋅12}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\}

+
    +
  1. 复共轭根对 <-> 实根 +在ncn_c个复共轭根对中任选一个复共轭根对a±jba\pm jb，转换为两个实根r1=a+br_1=a+b和r2=a−br_2=a-b。
    2. +
+

q(x′∣x)=pc×1ncq\left( x'|x \right) =p_c\times \frac{1}{n_c} +

+

此时的接受率写作:Split：α(x,x′)=min⁡{1,π(x′)π(x)⋅1q(x′∣x)⋅2}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. 实根调整 +选择一个实根rr，通过添加噪声ϵ∼N(0,σ2)\epsilon \sim N\left( 0,\sigma ^2 \right)调整该实根的位置为r′=r+ϵr' = r + \epsilon。
    2. +
+

α(x,x′)=min⁡{1,π(x′)π(x)}\alpha \left( x,x' \right) =\min \left\{ 1,\frac{\pi \left( x' \right)}{\pi \left( x \right)} \right\} +

+
    +
  1. 复共轭根对调整 +选择一个复共轭根对a±jba\pm jb，通过添加噪声ϵa∼N(0,σa2)\epsilon _a \sim N\left( 0,\sigma _a^2 \right)和ϵb∼N(0,σb2)\epsilon _b \sim N\left( 0,\sigma _b^2 \right)调整该复共轭根对的位置为a′=a+ϵaa' = a + \epsilon _a和b′=b+ϵbb' = b + \epsilon _b。
    2. +
+

α(x,x′)=min⁡{1,π(x′)π(x)}\alpha \left( x,x' \right) =\min \left\{ 1,\frac{\pi \left( x' \right)}{\pi \left( x \right)} \right\} +

+ +
+ + + \ No newline at end of file diff --git a/README.md b/README.md index ef0d906..3190470 100644 --- a/README.md +++ b/README.md @@ -17,29 +17,29 @@ $$ ### 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| - $$ + 辅助变量的先验：$u~U\left( -1,1 \right) $ + 提议分布:$q\left( \lambda _k,\lambda _{k+1} \right) =q\left( u \right) \cdot P\left( birth \right) $ + 接受率:$\alpha \left( \lambda _k,\lambda _{k+1} \right) =\min \left\{ 1,\frac{P\left( \lambda _{k+1} \right) P\left( y|\lambda _{k+1} \right)}{P\left( \lambda _k \right) P\left( y|\lambda _k \right)}\cdot \frac{P\left( dead \right)}{P\left( birth \right)}\cdot \frac{\small{\frac{1}{k+1}}}{q\left( u \right)}\cdot \left| J_1 \right| \right\} $ + 缩放因子:$\left| J_1 \right|=1$ 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| - $$ + 辅助变量的先验：$\frac{1}{k+1} $ + 提议分布:$q\left( \lambda _{k+1},\lambda _k \right) =\frac{1}{k+1}\cdot P\left( dead \right) $ + 接受率:$\alpha \left( \lambda _{k+1},\lambda _k \right) =\min \left\{ 1,\frac{P\left( \lambda _k \right) P\left( y|\lambda _k \right)}{P\left( \lambda _{k+1} \right) P\left( y|\lambda _{k+1} \right)}\cdot \frac{P\left( birth \right)}{P\left( dead \right)}\cdot \frac{q\left( u \right)}{\frac{1}{k+1}}\cdot \frac{1}{\left| J_1 \right|} \right\} $ + 缩放因子:$\left| J_1 \right|=1$ ### 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| - $$ + 映射关系: $\lambda _1,\lambda _2,\cdots \lambda _k,\theta ,\rho \longrightarrow \lambda _1,\lambda _2,\cdots \lambda _k,\lambda _{k+1},\lambda _{k+2}$ 其中: $\lambda _{k+1}=\rho \left( \cos \theta +i\sin \theta \right) $, $\lambda _{k+2}=\rho \left( \cos \theta -i\sin \theta \right) $ + 辅助变量的先验: $\theta ~U\left( 0,\pi \right) ;\rho ~U\left( 0,1 \right) $ + 提议分布:$ q\left( \lambda _k,\lambda _{k+1} \right) =\frac{\small{1}}{q\left( \theta \right) q\left( \rho \right)}\cdot P\left( birth \right) $ + 接受率:$\alpha \left( \lambda _k,\lambda _{k+1} \right) =\min \left\{ 1,\frac{P\left( \lambda _{k+2} \right) P\left( y|\lambda _{k+2} \right)}{P\left( \lambda _k \right) P\left( y|\lambda _k \right)}\cdot \frac{P\left( dead \right)}{P\left( birth \right)}\cdot \frac{\small{\frac{1}{k+1}}}{q\left( \theta \right) q\left( \rho \right)}\cdot \left| J_2 \right| \right\} $ + 缩放因子: $$ \left| J_2 \right|=\left| \det \left( \frac{\partial \left( \lambda _1,\lambda _2,\cdots \lambda _k,\lambda _{k+1},\lambda _{k+2} \right)}{\partial \left( \lambda _1,\lambda _2,\cdots \lambda _k,\theta ,\rho \right)} \right) \right|=\left| \det \left( \begin{matrix}\cos \theta +i\sin \theta&-\rho \sin \theta +i\rho \cos \theta\\\cos \theta -i\sin \theta&-\rho \sin \theta -i\rho \cos \theta\\\end{matrix} \right) \right|=2\rho $$ 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|$ + 辅助变量的先验: $\frac{1}{k+1}$ + 提议分布:$ q\left( \lambda _{k+1},\lambda _k \right) =\frac{1}{k+1}\cdot P\left( dead \right) $ + 接受率:$\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{q\left( \theta \right) q\left( \rho \right)}{\frac{1}{k+1}}-\log \left| J_1 \right| $ + 缩放因子: $$ \left| J_2 \right|=\left| \det \left( \frac{\partial \left( \lambda _1,\lambda _2,\cdots \lambda _k,\lambda _{k+1},\lambda _{k+2} \right)}{\partial \left( \lambda _1,\lambda _2,\cdots \lambda _k,\theta ,\rho \right)} \right) \right|=\left| \det \left( \begin{matrix}\cos \theta +i\sin \theta&-\rho \sin \theta +i\rho \cos \theta\\\cos \theta -i\sin \theta&-\rho \sin \theta -i\rho \cos \theta\\\end{matrix} \right) \right|=2\rho $$ ## 同维度游动 @@ -53,8 +53,8 @@ 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}\\ diff --git a/RJMCMC_ModelSelector.py b/RJMCMC_ModelSelector.py index 03ccd45..e4be382 100644 --- a/RJMCMC_ModelSelector.py +++ b/RJMCMC_ModelSelector.py @@ -48,8 +48,10 @@ class RJMCMC_Sampler: 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) + # 为避免除零错误，在分母上增加一个极小值 + epsilon = 1e-9 + self.current_sigma2_w = 1.0 / (np.random.gamma(self.lambda_a_w, 1.0/self.lambda_b_w) + epsilon) + self.current_sigma2_z = 1.0 / (np.random.gamma(self.lambda_a_z, 1.0/self.lambda_b_z) + epsilon) def _update_current_eigenvalues(self): @@ -59,6 +61,16 @@ class RJMCMC_Sampler: self.current_eigenvalues = np.roots(poly_coeffs) else: self.current_eigenvalues = np.array([]) + + def _update_current_a_coeffs(self): + """一个辅助函数，根据 current_eigenvalues 更新 current_a。""" + k = len(self.current_eigenvalues) + if k > 0: + # np.poly 返回 [1, a_{k-1}, ..., a_0]，去掉首项“1”，反转剩下的 + poly_coeffs = np.poly(self.current_eigenvalues) + self.current_a = np.real(poly_coeffs[1:][::-1]) + else: + self.current_a = np.array([]) def _cal_current_eigenvalues(self, a_coeffs_temp): """一个辅助函数，根据给定的 a_coeffs 计算对应的特征值。""" @@ -206,7 +218,7 @@ class RJMCMC_Sampler: # --- 计算似然 --- # 预测误差 nu_t - y_pred = C @ x_pred + D @ u_t + y_pred = C @ x_pred + D @ np.array([[u_t]]) nu_t = y_t - y_pred # 预测误差协方差 S_t @@ -228,7 +240,7 @@ class RJMCMC_Sampler: 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 + x_pred = A @ x_update + B @ np.array([[u_t]]) P_pred = A @ P_update @ A.T + Sigma return total_log_likelihood # 返回标量值 @@ -241,6 +253,14 @@ class RJMCMC_Sampler: k = self.current_k current_eigs = self.current_eigenvalues + # 获取实数特征值数目 + real_eigs_indices = np.where(np.isreal(current_eigs)) + num_real_eigs = len(real_eigs_indices[0]) + + # 获取复数特征值数目，只保留 imag > 0 的部分 (共轭对只算一次) + complex_eigs_indices = np.where((np.iscomplex(current_eigs)) & (current_eigs.imag > 0)) #type: ignore + num_complex_eigs = len(complex_eigs_indices[0]) + # 决定是诞生一个实数根还是一对复共轭根 can_add_complex = (k + 2) <= self.k_max add_real = True @@ -252,33 +272,51 @@ class RJMCMC_Sampler: if add_real: k_new = k + 1 + # 诞生实根的对数概率之比 + log_birth_real_forward = np.log(0.5) + log_birth_real_backward = np.log(0.5) + # 从提议分布 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_q_forward = log_birth_real_forward + stats.norm.logpdf(u_b, 0, 1) - np.log(2.0) + log_q_backward = log_birth_real_backward - np.log(num_real_eigs + 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))) + log_det_jacobian = 1 # 构造新状态 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 + # 接受率中的项为 q_backward(u) * |J| / q_forward(u) + log_ratio = log_q_backward + 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} + # 计算未归一化的后验之比 + log_post_unnormalized_current = self._log_posterior_unnormalized(k, self.current_a, self.current_b, + self.current_sigma2_w, self.current_sigma2_z) + log_post_unnormalized_new = self._log_posterior_unnormalized(k_new, new_a, new_b, self.current_sigma2_w, self.current_sigma2_z) + + # 计算接受率 + log_acceptance_ratio = log_post_unnormalized_new - log_post_unnormalized_current + log_ratio + acceptance_ratio = np.exp(log_acceptance_ratio) + + return {"k_new": k_new, "a_new": new_a, "b_new": new_b, "log_ratio": log_ratio, "type": "birth_real", "new_eigs": new_eigs, "acceptance_ratio": acceptance_ratio} else: # --- 诞生一对复共轭根 (k -> k+2) --- k_new = k + 2 + # 诞生实根的对数概率之比 + log_birth_complex_forward = np.log(0.5) + log_birth_complex_backward = np.log(0.5) + # 采样辅助变量 u = (rho, theta, u_b1, u_b2) rho = np.sqrt(np.random.uniform(0, 1.0)) theta = np.random.uniform(0, np.pi) @@ -286,25 +324,31 @@ class RJMCMC_Sampler: 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_q_forward = log_birth_complex_forward + stats.norm.logpdf(u_b1, 0, 1) + stats.norm.logpdf(u_b2, 0, 1) - np.log(np.pi) - np.log(1.0) + log_q_backward = log_birth_complex_backward - np.log(num_complex_eigs + 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)) + # |J| = 2*ρ + log_jacobian = np.log(2 * rho) # 构造新状态 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 + log_ratio = log_q_backward + log_jacobian - log_q_forward + + # 计算未归一化的后验之比 + log_post_unnormalized_current = self._log_posterior_unnormalized(k, self.current_a, self.current_b, + self.current_sigma2_w, self.current_sigma2_z) + log_post_unnormalized_new = self._log_posterior_unnormalized(k_new, new_a, new_b, self.current_sigma2_w, self.current_sigma2_z) + + # 计算接受率 + log_acceptance_ratio = log_post_unnormalized_new - log_post_unnormalized_current + log_ratio + acceptance_ratio = np.exp(log_acceptance_ratio) + + return {"k_new": k_new, "a_new": new_a, "b_new": new_b, "log_ratio": log_ratio, "type": "birth_complex", "new_eigs": new_eigs, "acceptance_ratio": acceptance_ratio} - 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)。 @@ -317,8 +361,8 @@ class RJMCMC_Sampler: 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 + can_remove_real = len(real_eigs_indices[0]) > 0 + can_remove_complex = len(complex_eigs_indices[0]) > 0 if not can_remove_real and not can_remove_complex: return None # 无法执行消亡 @@ -335,66 +379,129 @@ class RJMCMC_Sampler: # --- 消亡一个实数根 (k -> k-1) --- k_new = k - 1 + # 生成与消亡实根对应的诞生过程的对数概率之比 + log_death_real_forward = np.log(0.5) + log_death_real_backward = np.log(0.5) + # 1. 随机选择一个实数根移除 - idx_to_remove = np.random.choice(real_eigs_indices) + idx_to_remove = np.random.choice(real_eigs_indices[0]) lambda_removed = current_eigs[idx_to_remove] + # 2. 随机选择一个实数 b 系数移除 + idx_b_to_remove = np.random.choice(len(self.current_b)) + b_removed = self.current_b[idx_b_to_remove] + # 移除的根和b系数构成了逆向（诞生）提议的辅助变量 u u_lambda = lambda_removed - u_b = self.current_b[-1] + u_b = b_removed # 2. 计算逆向提议的对数密度 log(q(u)) - log_q_reverse = -np.log(2.0) + stats.norm.logpdf(u_b, 0, 1) - + log_q_forward = log_death_real_forward - np.log(len(real_eigs_indices[0])) - np.log(len(self.current_b)) + log_q_backward = log_death_real_backward + stats.norm.logpdf(u_b, 0, 1) - np.log(2.0) + # 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))) + log_jacobian_birth = 1 # 构造新状态 - new_a = np.real(np.poly(remaining_eigs)[1:][::-1]) - new_b = self.current_b[:-1] + if k_new == 0: + new_a = np.array([]) + else: + new_a = np.real(np.poly(remaining_eigs)[1:][::-1]) + new_b = np.delete(self.current_b, idx_b_to_remove) - # 接受率中的项为 q_reverse(u) / |J_birth|，在对数空间中为 log(q_reverse) - log|J_birth| - log_ratio = log_q_reverse - log_jacobian_birth + # 接受率中的项为 q_backward(u) / (|J| * q_forward(u)) + log_ratio = -log_q_forward - log_jacobian_birth + log_q_backward + + # 计算未归一化的后验之比 + log_post_unnormalized_current = self._log_posterior_unnormalized(k, current_eigs, self.current_b, self.current_sigma2_w, self.current_sigma2_z) + log_post_unnormalized_new = self._log_posterior_unnormalized(k_new, remaining_eigs, new_b, self.current_sigma2_w, self.current_sigma2_z) + + # 计算接受率 + log_acceptance_ratio = log_post_unnormalized_new - log_post_unnormalized_current + log_ratio + acceptance_ratio = np.exp(log_acceptance_ratio) + + return {"k_new": k_new, "a_new": new_a, "b_new": new_b, "log_ratio": log_ratio, "type": "death_real", "new_eigs": remaining_eigs, "acceptance_ratio": acceptance_ratio} - 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 + rho = np.sqrt(np.random.uniform(0, 1.0)) + + # 生成与消亡复共轭对对应的诞生过程的对数概率之比 + log_death_complex_forward = np.log(0.5) + log_death_complex_backward = np.log(0.5) # 随机选择一对共轭根移除 - complex_idx_to_remove = np.random.choice(complex_eigs_indices) + complex_idx_to_remove = np.random.choice(complex_eigs_indices[0]) lambda_removed = current_eigs[complex_idx_to_remove] - + + # 随机选择一对 b 系数移除，需要确保两次选的不一样 + idx_b1_to_remove, idx_b2_to_remove = np.random.choice(len(self.current_b), size=2, replace=False) + b1_removed = self.current_b[idx_b1_to_remove] + b2_removed = self.current_b[idx_b2_to_remove] + # 找到其共轭对 - conjugate_idx_to_remove = np.where(current_eigs == np.conjugate(lambda_removed)) + conjugate_idx_to_remove_array = np.where(current_eigs == np.conjugate(lambda_removed))[0] + if len(conjugate_idx_to_remove_array) == 0: + return None # 找不到共轭对，无法执行消亡 + conjugate_idx = conjugate_idx_to_remove_array[0] # 逆向提议的辅助变量 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) + n_b = len(self.current_b) + pair_count = n_b * (n_b - 1) / 2.0 + log_q_forward = log_death_complex_forward - np.log(len(complex_eigs_indices[0])) - np.log(pair_count) + log_q_backward = log_death_complex_backward + stats.norm.logpdf(u_b1, 0, 1) + stats.norm.logpdf(u_b2, 0, 1) - np.log(np.pi) - np.log(1.0) # 计算对应诞生过程的对数雅可比行列式 - 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 + remaining_eigs = np.delete(current_eigs, [complex_idx_to_remove, conjugate_idx]) + log_jacobian = np.log(2 * rho) # 构造新状态 - new_a = np.real(np.poly(remaining_eigs)[1:][::-1]) - new_b = self.current_b[:-2] + if k_new == 0: + new_a = np.array([]) + else: + new_a = np.real(np.poly(remaining_eigs)[1:][::-1]) + new_b = np.delete(self.current_b, [idx_b1_to_remove, idx_b2_to_remove]) - 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} + log_ratio = -log_q_forward - log_jacobian + log_q_backward + + # 计算未归一化的后验之比 + log_post_unnormalized_current = self._log_posterior_unnormalized(k, current_eigs, self.current_b, self.current_sigma2_w, self.current_sigma2_z) + log_post_unnormalized_new = self._log_posterior_unnormalized(k_new, remaining_eigs, new_b, self.current_sigma2_w, self.current_sigma2_z) + + # 计算接受率 + log_acceptance_ratio = log_post_unnormalized_new - log_post_unnormalized_current + log_ratio + acceptance_ratio = np.exp(log_acceptance_ratio) + + return {"k_new": k_new, "a_new": new_a, "b_new": new_b, "log_ratio": log_ratio, "type": "death_complex", "new_eigs": remaining_eigs, "acceptance_ratio": acceptance_ratio} + + def _propose_within_model(self): + """ + 随机选择实根变复根，复根变实根，或微调现有根，调用现有函数_merge_eigenvalues，_split_eigenvalues，_disturbe_swimming。 + """ + k = self.current_k + current_eigs = self.current_eigenvalues.copy() + + if k == 0: + return None # 无法在 k=0 时进行模型内提议 + + proposal_type = np.random.choice(['merge', 'split', 'disturb'], p=[0.3, 0.3, 0.4]) + + if proposal_type == 'merge': + return self._merge_eigenvalues(current_eigs) # type: ignore + elif proposal_type == 'split': + return self._split_eigenvalues(current_eigs) # type: ignore + else: + return self._disturbe_swimming(current_eigs) # type: ignore + def _log_prior_b(self, b_coeffs): """ @@ -404,6 +511,7 @@ class RJMCMC_Sampler: """ # return stats.norm.logpdf(b_coeffs, 0, 1).sum() + def _generate_stable_a(self, k): """ @@ -513,136 +621,344 @@ class RJMCMC_Sampler: 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. 计算先验 + def _log_posterior_unnormalized(self, k, eigenvalues, b_coeffs, sigma2_w, sigma2_z): + """计算所有参数的非归一化对数后验。""" 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 + + if k > 0: + log_prior_eigen = self._log_prior_eigenvalues(eigenvalues) else: - raise ValueError("Either a_coeffs or eigenvalues must be provided.") - - if log_prior_a == -np.inf: return -np.inf - + log_prior_eigen = 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) - 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 + log_likelihood = self._log_likelihood(k, self._cal_current_a_coeffs(eigenvalues), b_coeffs, sigma2_w, sigma2_z) + return log_prior_k + log_prior_eigen + log_prior_b + log_prior_w + log_prior_z + log_likelihood + def _purpose_b_swimming(self): + """ + 对 b 系数进行游动：对每个 b_i 添加一个小的高斯扰动, 返回b游动前后的先验 + """ + k = self.current_k + current_b = self.current_b + proposal_b = np.copy(current_b) + + # 对每个 b_i 添加高斯扰动 + for i in range(k): + b_proposal = proposal_b[i] + np.random.normal(0, 0.1) # 小的高斯扰动 + proposal_b[i] = b_proposal + + # 计算新b的先验 + log_prior_current = self._log_prior_b(current_b) + log_prior_proposal = self._log_prior_b(proposal_b) + + return proposal_b, log_prior_current, log_prior_proposal + + def _purpose_lumbda_swimming(self): + """ + 对 lumbda 进行游动：对每个 lumbda 添加一个小的高斯扰动，虚部也要有扰动 + """ + k = self.current_k + current_eigs = np.copy(self.current_eigenvalues) + proposal_eigs = np.copy(current_eigs) + + # 对每个 lumbda 添加高斯扰动 + for i in range(k): + eig_proposal = proposal_eigs[i] + np.random.normal(0, 0.1) # 小的高斯扰动 + proposal_eigs[i] = eig_proposal + + # 虚部也要有扰动，共轭的虚部记得相等 + for i in range(k): + if np.iscomplex(proposal_eigs[i]) and proposal_eigs[i].imag != 0: + imag_perturbation = np.random.normal(0, 0.1) + proposal_eigs[i] = proposal_eigs[i].real + 1j * (proposal_eigs[i].imag + imag_perturbation) + # 找到共轭根并更新 + conjugate_idx = np.where(proposal_eigs == np.conjugate(current_eigs[i]))[0] + if len(conjugate_idx) > 0: + proposal_eigs[conjugate_idx[0]] = proposal_eigs[conjugate_idx[0]].real - 1j * (current_eigs[i].imag + imag_perturbation) + + # 计算新lumbda的先验 + log_prior_current = self._log_prior_eigenvalues(current_eigs) + log_prior_proposal = self._log_prior_eigenvalues(proposal_eigs) + + return proposal_eigs, log_prior_current, log_prior_proposal - def _merge_eigenvalues(self): + + def _merge_eigenvalues(self, current_eigs): """合并游动：选择两个实数根，确定性地合并为一个复共轭对。""" 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 + # 随机选择两个不同的实数根(可能是虚部为0的复数)，并删除他们 + real_indices = np.where(np.isreal(current_eigs)) + if len(real_indices[0]) < 2: + return # 不足两个实数根，无法合并 + idx_pair = np.random.choice(real_indices[0], size=2, replace=False) + idx1, idx2 = idx_pair + lambda1 = current_eigs[idx1].real # type: ignore + lambda2 = current_eigs[idx2].real # type: ignore - # 确定性映射 -> 新的复共轭对 - a = (lambda1 + lambda2) / 2.0 - b = abs(lambda2 - lambda1) / 2.0 - new_complex_pair = [a + 1j*b, a - 1j*b] + current_eigs = np.delete(current_eigs, [idx1, idx2]) - # 构造提议的特征值集合 - proposal_eigs = np.delete(current_eigs, [idx1, idx2]) - proposal_eigs = np.append(proposal_eigs, new_complex_pair) + # 生成一对共轭复根 + new_theta = np.random.uniform(0, np.pi) + new_rho = np.random.uniform(0, 1.0) + new_eigenvalue = new_rho * (np.cos(new_theta) + 1j * np.sin(new_theta)) + current_eigs = np.append(current_eigs, [new_eigenvalue, np.conjugate(new_eigenvalue)]) - # 检查稳定性 - if np.any(np.abs(proposal_eigs) >= 1.0): return + # 对b进行游动 + proposal_b, log_prior_b_current, log_prior_b_proposal = self._purpose_b_swimming() + + # 计算提议密度 + # 正向过程是任选两个实根，合并为复共轭对 + log_q_merge_forward = -np.log(len(real_indices[0]) * (len(real_indices[0]) - 1) / 2.0) - np.log(np.pi) - np.log(1.0) # 提议密度 q_forward(u) 在 (0, pi) x (0,1) + # 逆向过程是从复共轭对中任选一个删去，然后抽样两个新实根,提议分布就是均匀分布 + num_complex_eigs = len(np.where((np.iscomplex(current_eigs)) & (current_eigs.imag > 0))[0]) #type: ignore + log_q_merge_backward = -np.log(num_complex_eigs + 1) - np.log(2) - np.log(2) + + # 计算对数提议比 + log_proposal_ratio = log_q_merge_backward - log_q_merge_forward + + # 计算后验比 + log_post_unnormalized_current = self._log_posterior_unnormalized(k, self.current_eigenvalues, self.current_b, self.current_sigma2_w, self.current_sigma2_z) + log_post_unnormalized_proposal = self._log_posterior_unnormalized(k, current_eigs, proposal_b, self.current_sigma2_w, self.current_sigma2_z) # 计算接受率 - # 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) + log_acceptance_ratio = (log_post_unnormalized_proposal - log_post_unnormalized_current) + log_proposal_ratio - # 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 + acceptance_ratio = np.exp(log_acceptance_ratio) - log_acceptance_ratio = (log_post_proposal - log_post_current) + log_proposal_ratio + return {"a_new": np.real(np.poly(current_eigs)[1:][::-1]), "b_new": proposal_b, "log_ratio": log_proposal_ratio, "type": "merge", "new_eigs": current_eigs, "acceptance_ratio": acceptance_ratio} + - # 接受或拒绝 - if np.log(np.random.rand()) < log_acceptance_ratio: - self.current_a = a_proposal - self.current_eigenvalues = proposal_eigs - - - def _split_eigenvalues(self): - """分裂游动：选择一个复共轭对，随机地分裂为两个实数根。""" + def _split_eigenvalues(self, current_eigs): + """ + 分裂游动：选择一个复共轭对，随机地分裂为两个实数根。 + """ 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 + # 随机选择一个复共轭对，并删除它们 + complex_indices = np.where((np.iscomplex(current_eigs)) & (current_eigs.imag > 0)) #type: ignore + if len(complex_indices[0]) < 1: + return # 没有复共轭对，无法分裂 + idx_to_remove = np.random.choice(complex_indices[0]) + lambda_removed = current_eigs[idx_to_remove] + conjugate_idx_to_remove = np.where(current_eigs == np.conjugate(lambda_removed))[0] + + if len(conjugate_idx_to_remove) == 0: + return # 找不到共轭对，无法分裂 + conjugate_idx = conjugate_idx_to_remove[0] - # 映射到两个实数根 - u = np.random.beta(2, 2) - lambda1 = a + b * u - lambda2 = a - b * u - new_real_pair = [lambda1, lambda2] + # 删除选中的复共轭对 + current_eigs = np.delete(current_eigs, [idx_to_remove, conjugate_idx]) + # 生成两个实数根 + real_eig1 = np.random.uniform(-1.0, 1.0) + real_eig2 = np.random.uniform(-1.0, 1.0) + current_eigs = np.append(current_eigs, [real_eig1, real_eig2]) - # 构造提议的特征值集合 - 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) + # 对b进行游动 + proposal_b, log_prior_b_current, log_prior_b_proposal = self._purpose_b_swimming() - # 检查稳定性 - if np.any(np.abs(proposal_eigs) >= 1.0): return + # 计算提议密度 + # 正向过程是任选一个复共轭对，分裂为两个实根 + log_q_split_forward = -np.log(len(complex_indices[0])) - np.log(2) - np.log(2) # 提议密度 q_forward(u) 在 (-1,1) x (-1,1) + # 逆向过程是从实根中任选两个删去，然后抽样一个复共轭对,提议分布就是均匀分布 + num_real_eigs = len(np.where(np.isreal(current_eigs))[0]) + log_q_split_backward = -np.log(num_real_eigs * (num_real_eigs - 1) / 2.0) - np.log(np.pi) - np.log(1.0) + + # 计算对数提议比 + log_proposal_ratio = log_q_split_backward - log_q_split_forward + + # 计算后验比 + log_post_unnormalized_current = self._log_posterior_unnormalized(k, self.current_eigenvalues, self.current_b, self.current_sigma2_w, self.current_sigma2_z) + log_post_unnormalized_proposal = self._log_posterior_unnormalized(k, current_eigs, proposal_b, self.current_sigma2_w, self.current_sigma2_z) # 计算接受率 - # 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) + log_acceptance_ratio = (log_post_unnormalized_proposal - log_post_unnormalized_current) + log_proposal_ratio + acceptance_ratio = np.exp(log_acceptance_ratio) + + return {"a_new": np.real(np.poly(current_eigs)[1:][::-1]), "b_new": proposal_b, "log_ratio": log_proposal_ratio, "type": "split", "new_eigs": current_eigs, "acceptance_ratio": acceptance_ratio} + + def _disturbe_swimming(self, current_eigs): + """ + 扰动游动：保持维度不变，保持实数根和复共轭对的数量不变，进行扰动。 + """ + k = self.current_k + + # 对特征值进行游动 + proposal_eigs, log_prior_eigs_current, log_prior_eigs_proposal = self._purpose_lumbda_swimming() + + # 对b进行游动 + proposal_b, log_prior_b_current, log_prior_b_proposal = self._purpose_b_swimming() + + # 计算提议密度比 (对称提议，密度相等) + log_proposal_ratio = 0.0 + + # 计算后验比 + log_post_unnormalized_current = self._log_posterior_unnormalized(k, self.current_eigenvalues, self.current_b, self.current_sigma2_w, self.current_sigma2_z) + log_post_unnormalized_proposal = self._log_posterior_unnormalized(k, proposal_eigs, proposal_b, self.current_sigma2_w, self.current_sigma2_z) + + # 计算接受率 + log_acceptance_ratio = (log_post_unnormalized_proposal - log_post_unnormalized_current) + log_proposal_ratio + acceptance_ratio = np.exp(log_acceptance_ratio) + + return {"a_new": np.real(np.poly(proposal_eigs)[1:][::-1]), "b_new": proposal_b, "log_ratio": log_proposal_ratio, "type": "disturbe", "new_eigs": proposal_eigs, "acceptance_ratio": acceptance_ratio} + + + def run_MCMC(self, num_iterations): + """ + 运行 RJMCMC 采样器指定次数的迭代。 + """ + # 存储采样历史 + k_history = [] + a_history = [] + b_history = [] + sigma2_w_history = [] + sigma2_z_history = [] + eigenvalues_history = [] + + print(f"Starting MCMC with initial k={self.current_k}") + + for i in range(num_iterations): + # 随机选择一种转移类型 + # 这里的概率可以根据需要调整 + move_type = np.random.choice(['birth_death', 'within_model'], p=[0.5, 0.5]) + + proposal = None + accepted = False + move_name = 'None' + + if move_type == 'birth_death': + # 决定是 birth 还是 death + if self.current_k == self.k_min: + proposal = self._propose_birth() + elif self.current_k == self.k_max: + proposal = self._propose_death() + else: + if np.random.rand() < 0.5: + proposal = self._propose_birth() + else: + proposal = self._propose_death() + + elif move_type == 'within_model': + proposal = self._propose_within_model() + + # 处理提议 + if proposal and 'acceptance_ratio' in proposal: + move_name = proposal.get('type', 'N/A') + if np.random.rand() < proposal['acceptance_ratio']: + # 接受提议，更新状态 + if 'k_new' in proposal: self.current_k = proposal['k_new'] + if 'new_eigs' in proposal: self.current_eigenvalues = proposal['new_eigs'] + self.current_a = self._cal_current_a_coeffs(self.current_eigenvalues) + if 'b_new' in proposal: self.current_b = proposal['b_new'] + # sigma2_w 和 sigma2_z 在这些提议中没有更新，保持不变 + accepted = True + + # 存储当前状态 (无论是否接受，都存储当前链的状态) + k_history.append(self.current_k) + a_history.append(self.current_a) + b_history.append(self.current_b) + sigma2_w_history.append(self.current_sigma2_w) + sigma2_z_history.append(self.current_sigma2_z) + eigenvalues_history.append(self.current_eigenvalues) + + if (i + 1) % 100 == 0: + status = "Accepted" if accepted else "Rejected" + print(f"Iteration {i+1}/{num_iterations}, k: {self.current_k}, Move: {move_name}, Status: {status}") + + return { + "k": np.array(k_history), + "a": a_history, + "b": b_history, + "sigma2_w": np.array(sigma2_w_history), + "sigma2_z": np.array(sigma2_z_history), + "eigenvalues": eigenvalues_history + } + + +if __name__ == "__main__": + # 导入生成数据的模块 + from generateSimData import simulate_lti_data + from generateGroudTruth import generate_ground_truth_system + import matplotlib.pyplot as plt + + # 配置 matplotlib 支持中文显示 + plt.rcParams['font.sans-serif'] = ['SimHei', 'Microsoft YaHei', 'Arial Unicode MS'] # 用来正常显示中文标签 + plt.rcParams['axes.unicode_minus'] = False # 用来正常显示负号 + + # 1. 生成 Ground Truth 系统和仿真数据 + print("--- 1. 生成仿真数据 ---") + # 真实系统阶数为 2 + A_true, B_true, C_true, D_true = generate_ground_truth_system(dx=2, rng_seed=42) + + T_steps = 400 + sigma_proc = 0.1 # 过程噪声标准差 + sigma_meas = 0.1 # 测量噪声标准差 + + u_data, y_data = simulate_lti_data( + A_true, B_true, C_true, D_true, + T=T_steps, + sigma_process=sigma_proc, + sigma_measurement=sigma_meas, + rng_seed=123 + ) + # RJMCMC_Sampler 需要一维的 y 和 u + y_data_flat = y_data.flatten() + u_data_flat = u_data.flatten() + + print(f"仿真数据生成完毕。 y shape: {y_data_flat.shape}, u shape: {u_data_flat.shape}") + + # 2. 初始化 RJMCMC 采样器 + print("\n--- 2. 初始化 RJMCMC 采样器 ---") + sampler = RJMCMC_Sampler( + y=y_data_flat, + u=u_data_flat, + k_min=1, + k_max=4, # 探索的最大阶数 + initial_k=3 # 从一个不等于真实阶数的阶数开始 + ) + + # 3. 运行 MCMC + print("\n--- 3. 开始运行 MCMC 采样 ---") + num_iterations = 20000 + results = sampler.run_MCMC(num_iterations) + print("MCMC 采样完成。") + + # 4. 分析和可视化结果 + print("\n--- 4. 分析结果 ---") + + # 丢弃早期样本 (burn-in) + burn_in = 1000 + k_samples = results['k'][burn_in:] + + # 计算模型阶数的后验分布 + # 使用 bincount 统计每个 k 出现的次数 + k_posterior_counts = np.bincount(k_samples, minlength=sampler.k_max + 1) + k_posterior_prob = k_posterior_counts / len(k_samples) + + # 打印后验概率 + print("模型阶数的后验概率分布:") + for k_val in range(sampler.k_min, sampler.k_max + 1): + print(f" P(k={k_val} | y) ≈ {k_posterior_prob[k_val]:.4f}") + + # 可视化 k 的后验分布 + plt.figure(figsize=(10, 6)) + plt.bar(range(sampler.k_min, sampler.k_max + 1), + k_posterior_prob[sampler.k_min:sampler.k_max + 1], + color='skyblue', alpha=0.8, label='Posterior Probability') + plt.axvline(x=2, color='red', linestyle='--', label='True Model Order (k=2)') + plt.xlabel("模型阶数 (k)") + plt.ylabel("后验概率 P(k|y)") + plt.title("模型阶数的后验分布 (After Burn-in)") + plt.xticks(range(sampler.k_min, sampler.k_max + 1)) + plt.legend() + plt.grid(axis='y', linestyle='--', alpha=0.7) + plt.show() + + + - # 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/equ.afx b/equ.afx index f57abc28247afac5552183558a11b76e2fd60fb5..6e04110a3d8c938399ea0aa10a56c42c3aa6b957 100644 GIT binary patch literal 131072 zcmeHw37AyHweW2&7!(r+f*1+vJoF_R7aR!csFmQdo8zEv%2esd5d~(bd%)7DaUWY|E6Pqv*;)1PCR}7AO7^oOW(Wa zXwuM;KdZmsob|D+yohI~L+_&BarpmhNp+n}H>X1^pVT_<^n*(N4xO5Mb;9aR3h>0dP9NV1OY2 zX8@cDFcjb{fbRhOFMzWF&H*?V;5>j~0K)-B0E`3}1#mvVcL6Q{xDenXfbRi(A7C`V z4*-4$a52D-0LB2+0*nRtFMx3X;{kpQFah8afQbN?0$c`gIlvVFR{~rGa5ca+077#n z0ZazC4qytvQ~;sP*8}_%;0A#I4e&F7p93TTN~!P@LfTIBZ9pD=P-vkiO;~0Q%0UQf(9Ki7aCjj&TI1!)*;0b{2 zIHKQuivo*)_fG=+H^Bb@{2Aac0E+?Y0iFU_0`N4zQh;Ru%K;hyQUE?c8sHg#X8~3K zJO{86;CX<*0=xjw2+#!ZBETwu)c|V%)&jf)@G`(F0Dl8m2k@7w*WQ(Yy@~4U=zSQ0RI5k4Dc?%djRhP{1adcz*c|{0JZ^a2lxW?448QwX2ph&;ibYL7`%KP8!*{@++qD=r%7xMCbSuF|WqR$*)+x`0W z!}1rQT+|bNi{IPQI{}BCx%!SXH!j{4_zOG^c=brAC)zxF;4!0b8#H_N>_Owl|8dZZ zd1Dgoa-9QpPoOxYE* zY@{Adu@*&!rAv~uE<}s_ljIvl9+ii?At;YyF~IpJbht=S5#zS}l})NzP-sb24v#?7 z+a2wJj6rK zj|@cxb1}>K;W#RmZ!qN8`G{hU6g1aTjD-G9E+3qbQdtRpUt17ZKZpai7lSyF-?vUt zY_L7~b6#tV?YK25)z+AH=+E`X7WsWnPNkGmSIgA{my)FAsHU_?tWTdK44dlLL)FJc4R2Q&Cfc99DS=NmPR< zOo=Rbc`kesVHBR`jk6^w1MQ!Z`IRi2cWz59;3%@x;#>pUyDYzpOdnS=Wk056sOckU z85bv-^pPaZ2d?0o8l8RJ;8QAN&e%{|qiKi`(^L*BInW6vLxK8RBHb#s!z_mE(m+D0 zLZLxOx%j3T44}$b^tB8-tS{0c&cs1qBk(nwz5=oSDcFyDo|i zT6203>v6G^rp?BXOBMXqnkfxifr`XVRFDV;GjE%QxAk~4v++t;o71*e<&!lOqG)X{ zCq_t*TrQX3lqkD!VW0^m4@GkWB8g}kkeLv3XD_W;du?j~f@N?{(9}BDGNNH~?Ta)o zGsf~~)1w&EzPNFAMJ?=^&`LLasi4T45nU*f&Rh2k!Iv+}Z$-Xqa)zLuZbl3v&ZpAGU)%>_a`xgoO%|Msy^+ z5E*HC3mY4gl`FdXGxf-a65Dw+=lg zvHRs)eYbk)(8R3c>!4)D8Ht)ReNpF(MB|3%{Fq40TG#|7Zx2m0&Rq?aN1qENkQjB| zNLffRSP%o;-dK)*1HEO@r{Q>eL)4 zfp+Geb!!}3UVBs@5BrL3{e$ODn}1QF{{v5WZgsB<6B9-}4JCkKZ0ico6*$MXz6h0J z6tS&qp$?2oChi+56ZL-{lDN0iYN!l}{zv>(CMMjm(nCVUwx;0qK11QrpT(nLiIcmG zgOb%_0iW%0WI^4ho#V(9k%$R4Rj1Ysas>y6ee?^N5UgADJG;Iqiee}sg-L6qiDvX;v|T%(wQcc2;i6bZAxG+;iCl%kMOm$3SAVcU#Y1YI z_(K^VOGYP2&u#QwOY32`ijn3%LBx8>d;t4Z+QTxgp7OCF(=S0Rn~8XoF@Fk=xP|Vr z@u)96w!zA~HRs>P5}36M64gL?!w*2``dtp_ZX8p!7!PA9D0W zdgOAsgcMxeX%4mF*>y()Eh_mriD0Qx$Ce{eUP4mKV7%5mMSG+@5Rx^;$vg@}63rnk z5;-6`7MQ8iW>Q2~{*F_b9yd-X-1Rw8#)01IQ6I6BtL|Kgp}HTS%mPq6ZllIh7W| z1%CU8l%fV9@2$L~u`#yhsDlFz)nFOpTtUH{%J+*&@>-y) zG1aYF(v4)1O{YwnmlBieE&M~bkj=FGO5L;dc&+sjY*AIG8N8OhPpxZ}gYu?fgKj8RhOxCMPfN-;EL9xQ_@39A)zOwpsfFADI%~K9UAPlIC?1jD&q-370U@EwE){)AwvBbrbh_wN&sEG595c@EQx@svavAL3Lz8K6&tns%U$PjAq zJ~bz?#@!IHfTH(JD*H^oUxYLqEjByz^Ix4HG)8a;FT zitSNM)`6olc2fM5$by&W!Y2_%;VC1Npf%_Rqw{ak`8UTYxwB;@(%ev5b4XF3gB%Jb zL($JsIZk_=f*Ig!GYBac-!y{(jN|Oh;#p&2d;Je_rmMrjR907g*w1Xy{}!l66%k2e$4RKnVvw!LE4;E~MM0xS?1Zhbde zo1+qh^q@r|$0<>^)IbxWX+Ys$bDLVBlG9$>8h~I$_kV}lk+T}go6^3xadu^e*#@*i zOREvSR8ZvI+ApAQz>(;lA&Bn(b{u;1iVS}2oqIrboJe@Z##fpGnDPCE1!7zuEHVbc z9`>Q02*yI?0BSqD$u1V%{~g``?L_x~OKXNoO$-}m*vW})D}x5A&{YJ^8W~u;P`D_T zu{HU!^SHt+d?;L$aEG0dV}pu^)VydI$!za#Ms~u2=qabGuwQ%kQ*taa^1^#FkGNBw zW#dtaVs%xeSv}Ik1!Y&6KVsQ9B8JkODmU_Js%4V2)b^Kq(f!|)Shn@PqK-?J)kmvJ zq|hC+62V9&723#v%UZn9hD_3n#De98#tQ=73b%>T{olZ%K;KQ-DcNo>zS55CxKPw` zjHzfAX61M-#<78PE`1zZHlFL!pfo&ZtrLnCbxYFHQgvWO2GKTHc`t=Ax7G)%z(_r1 zg2{VNO&Z%>sfcGC@sW3avTv z8qFas5;YZPRE|QejTYWc2M6wWk5H zhITT#*C;1Op~7(~X$I4ol83Rmagv1bIILO3;22Exmc}@joc8)-E9=_fK0ggG_#Q}e z?GTNOD$qnDWv24KlPXboN`p>@Qb9#Bh%0a&fN{>=v2JGVGHeJVOv682!)p;|PCSgq zVGTch$6$r@%cz$2@d$&!)8-T`n8gL!-6U04v^TU!JIB7E;uQeXE z(^$S-@rLQV$7amM2$r#2)sjYKbJTjuvTK-nI>&s=#W&4GnIsuj;0YF87mCuNEUb4! ze@^<*-6GBEszD%Y`|LTD13)lvcm>V_G-D~Gbde8(?alGQFd@=ia{*~CLF6bBYmKJC3aG{<#yFRp_WEP9Z-&nyWBr^KQtRMO*@3s| zxQZPglqpt;)}T-=4c4IO>RPZ3Qda+1VK&8u&bGj%XkN}@+OIARlrW+sGyp1u(XG8i z5_CDDNM+U(lHX|VB15wh#Mvr4jAOzIinG$F8v~;U zQTCI=1fw?R0gPZ7%T+CDR5nMwtt_KKwJeL)HVN^vnpNmDKxquKlF(O^N3k&SSM13Z z8|vVMqd>i^q`3=@<7`TNEE~#{jihDVE;AaF6xq|3G04&MU94x%;>$jcp?z?NIkDdn zv>uMP$oSzs&|7l-ej8F#cAu0^!tKbIsNQ@^gZoHT|MZ4IiyGW118dSx9JsB)9kqS` zgh<5d-rVYaa^9mUchV!PycZ_*^xdiebk7{mILEfW=;0`0Th~Gz7?(_(GG8VpU3#B)Z>QCw z@_k7xWY3BPyN~>IFc?n6=eaM2M_v(grRFDD1B99)Kr0(h? zwo`I{5mM4Lq@XHVT!>^j?Uv~w}{GhOc_^?Ymd9Ln7&II z&t{~)jmT#4m))&Pk}84c&N$J~6nlKt;*J7qN(#pX;>R#FFi?Gty_A*G9F$F548~R= zX044VPQ%GlRluASHiO$7lqojagX2g?vJ1N0&)a8wSy}kRvP)=#Sat<18;OOvGMZ~C z#&X4UVi|`sH7~hF@)NAGhS^u9^vkSUT$SHNxP5Bco>Vw)Y&!hljAEFvseNF324kb` z1&@u=fY36Vk$j`y0?GylWUJt>EK^o4QDE+-E3jb`_5$i-NK`z+h6^091GK)q!@Y7<-qBh1+ zn6)epO$)=NyfoF2dn6Vo&x$T!JP2Q! z8nC)mIW#aE8;)0lE}WA>8YuJ&BL@~E^?X~HO-tKZ5prtAqrme(3hftWM=pPnwX{g& zBw6gmFZ}sRl#MZ@cimXewvFIK7ipu$SDZYTho-?t<5Zkk1IWa{Ts;FGA^ntioG?3+ zsURH;xg(=t6OjnIzM>*@Ri6gB@`5x_V#hE(z@{Cl9R1~OA{YOafX?A`bqrBs8S@+P zmu`{qI{*fD2w$ar=-W=tF>o>zMxG0eY3$N~@GzT^wpigW7h5E=Z^&OnjTVWVY);w! z(d}Yz76{5+Zu1RBqQRqSQf^Jhm!-YB^BAZL6THyd-vI3Y%2W%#o3^bLg_~- zAam9QhLc@Rj{J;xL846%C$|ZL z5;B5F2=CQ3zc5CnQ=&jDQ}z$BO!;g~b6Z;b-;__%Vc}3|k;qBk(4lw!6VA8hKhLIv z-!-r?iX+JOMSf#Rk_~i9NahL}n}!3@KI#-WwpvOXCf9!XR;IBrd}DoI7_r>5oJ%$_ zeXA^$`fA`KN01gHA1S`liatZj{|t>rQM`{t4`U#Ow2w#%BPmorhNO1NiD|B9bJGuW zpNb7;Z8}!mvk{gySXnR$Gf%mO*kCpz`?#7s3v~k!xAo3_ZKKpf#r`Jgwhaq#6{#3% zpM)Nqbz%p1ByEXIsL)U0Z>z?#>pqI*(Tv2AEcQ4q-Mw_}VV+n)xsGRaLi8xcik_}#xKa9*gNa3dC6WK}9W8ka> zho*do6-Fc#s?bhlPIAp!AF?^cLhv8tY0?0>rmMb{VG&C*-jMdd9%Taca?k#Q*Gp{B zz*mf=a$PcCp_qs;E?z;I4j_}0sx*>3lVu96kbajj^%T2?^^~Qi64-IFqe8L;aCNMZ;S0xaMO&3)~axZCZf zg_| zD(Ll{BmMq2ragDb%Gue(T1><~7_d_&M$Y?GKAPVBbHD#HzUPjd_qo5S7hDXqg;$(Yam*vxuJY}^GUPl}Qq{;JUnJ-7O?_x)A({mFBuzw(xU??dxE_k#hiK^?&7 z=I$>+gZGHaFZi*lyW!FG5BqB`J48IX-FNzSS?<;)Ztzd|T`%8FjQ*~_`^tMoVxZr6 z>U$zl<2PT>GwoKNko5aD9tw$V5BN>b_DH)ozp~h$J-ln${b0#zzxmZJkofFPzj;=N zv^)LIfB0j2Z1DlH&HmVR@4_Wbue=Y5O}@M2OWXZf`@iJ_h95&0uR~(VXa4TT{|2`{ zog7Oy-munp-#Bu=^tijkJy5d_NRL~$QY5~TZanol-+g#yw{*<{xEZSIU)|HY`>cSH zH3z3_HbI?-t9n4N>qgc9m@qPl;?S9`bRe=0} zzr(l=fpOir%3*}{^x)BHWTjJoY}&0$9}1OENW1l0o=@E~{KT}o zX5g@tb0H))EPOF860xO6p6=nJZYzH6H9dHI+U?fqaVUXys)oGmIiK`~#EK8SrjB1v zyW5_36U|cqjZu9LO7z{*DecZk_ej)?-2oJv+%r+L0LZfK=))7^)-8g>gu@c!di>Uh z%2Il7TtEDu-d&#_;5$_-KK5##y9h=)WS2MVSGU3KT5oRic3-*BcLz^>*PC^gplYYr zyv9?38g+A?fy5{9=#*vI04+lwb zw(nSuA~E-XV^Sg!TQhKcs%9)m%X^1En%cb`nGJzCQV4 zf8u^$1@T)l)1UB{-@}bp=`%19tKYIS)&ChNSu=2OYQnCWo*SRGc+jfRH+!H~%iA`Me`9O{de~`I zr;XlmAI%WDEa|i0x+xnI{rdeJ0n}eyiv=LB$8wU6GZ+fskBt(N_$c?kAh(3&NdwC_0jJBA}8Wc@QVIRb> z<6sG0CS2Ivcf^{0DYz|7$!BYxrqvOgp5#{@D04%E^3|IcM{<1@farPF^D* zQphf|u`=Pu#sZqSoJ>Q*fECwRM=?UWbiDFX?#4;%bUybGWrkF+@8z*hXjhKxk%9(J#+kONuzSO)OyOYYnXaE3qR%Jn`WGtB*O}viR-#hlp`vq$>4OT!^#Eo$-IqmhwX5S7a^VuD4<0Xtd z*OG1|i){K|(!9*Dau5b*n=wr=iwm^7NjhS*eY8mAVw&**zOK06Fr?CPzhT-_-EZrW z!Pt_tG>8Of^Y)b56Tp)uyKSYF3}p4kQ}tsO(FT+o7XnOILs2n*@|fJ zYcgSQuAn$u9UGCcOvX5#w8kSa8)y#<=DH=PJrFV2{}eX{Hl>x#c>p6=#=NX0jmqY1 z&aKMDH_f?Ak_^9{tZjy|lf|A~v7ru5I11FuN}9XiIL;;_B4dW3Oqp_bHMh%*#zeUl z3H#vwonm_yU*dcW#hE*_i1;b%v8BKd4~K5a;rDGwP1$`?THL6GiR#U#G`Npc^-pgY zw5Y+IGO#B7#DUuy+)>;2Pl!aU?#->i7 zU%u6M<5TB*vyQKWlCE>SnlpVUIpuz@al>xvy^OpcDYPM|L~GS|8YWV(jV-Al(;R1lCG;ISLdWv zB;>nA)#8xj;69eHM$^oADXj=WG8svieaM2M_v(gr6y<^=3W!URmgnlzVi}E<++T#0 zGz}@h@s)fWyp)=J#em}|x#nlaQB;jC63Vy(!GQW0*4TUAjcV z-=tU`4-!(IsV8K4iG;49gee9uIXaJ{-6AU2F=bpmuDxsb_%5lC%}9M4kSyIYqe zRRY&$IML7)dwkX6jsj~+3dbd6plG1_9D6x2D4Vz#jIBb zWjcuV;5gEelxMo!&)a8wSvmiJWtY$fvFr+3X8B*6M6yuxLi)>QSY-`!D;cF8>IoEWi})EMn9gH4Gzdwk+eryQ`Tg# z>~8W^QK*k*#*!{FI<549w4y;xs6G@L$gwH01&)Ig4ueA$idy0@6eoaXG}g<|fWxJh z5Gxc7jBrtC6|F8utIJMqb-6uSfWYqul$ivMjY>aJ!AwS4P-CxtXn=odx~ zEJo`2wlJHPwzDGS)Qm@g=YbU3FU*cy{vvB>k;qB1*o$BI^OYzYV@U71v7T)k!HF)? zMvbpHc`gr4gOA3kII{+jiGjI#20TLgDe*XAb|zCnIv8?CM#Cl|5p;b;Md+$N4RqxN zX`sZ8VSIp1J5)LP%iBaQ{wo2U!|Cc6qQ)}jH{dVbBI9=e4D1lTO8d~aot$IfWGIY0 z7aG&pr2*k#HY07Z!e1`7NM_%Vzla(w5;@tNvi+w?Dilc_SB8BX^wGC9-fBrk#A&fC zj1dnf2V6@kvU7tZsnYiDv%%O@`UQ%!K_7+Ek5WM9tP2b$yPO^c+E|5`D3u(~?Gccb z7Y^q@W>TDk>8$JwYbqw9ayn@>`scZ$f1W$(bm*P`g!8TW&$H>^cMWWe;s~;R zk>B<@HVp@)eI_2ofa6;&Ar+TtzkDmx*ciUCzAs6s-KNUzSF$-&7MjwRjRDRfLzIJ4rXGkj6*gH<`8O`Hp7JhY)3Z6h6Zx*a z)7Olwv4p{TR8D!BM4YK$$p$aUA4X;#q;S*niR>ilvF)War^yQ)q?syNs!)*fp%DEH#zDj$;>!GecCUbz&_?a2&P|+i{jsg#AQ1j)E#}OUsk0dG@)`}mmuHltzqPI=Ot`hjP-0OpFF~mZw>s<0g zwnxufkP4ZMtY^S=1?52f9YcE84ezlzMGjbw;72M0bA1S-H>%Iq5_f%ihVP7eyjx=Q zH{MOTqxy7C%sBLrl)HZT4sYt_&T&Y*Xm}G&iulkuxlxljQQM;&8&yS#$8bIq$t_p zuNuA3bE_YF-(Pj#pFDT^D{uMtJ~YpBKN#>D)B${M?*0-qc#o+3f*-588y;Q%u)p@Q zL&T%oeW!1iiG4vyX|>5(L4pv7}e*XMBgo)((a6O zk3`Mb9YC?kJrgwxfGpdNK0GmQ-6BX#I4m))$8UY8ET#9x^~3+^-Sz1KzEid0W3TqP zi(sTfc6qaYbsOBS_2xEj_mvBMcktAAy;*1B?OLySji&-N>gGHHiBI6sDUW&M4+m-t zo_YtM`ei_4*^j*D#of~G(gnRdB+C=m&xy|-4wB$(-?1D;V(tOQq(maNX5jc#%~+6@ z_YQwFwR=0r&$RmH)E93Z?z>-k=omlVRKMJPee%Wr#QnYs;0MG3< zNXTfz-HYWP@`#A!r z{~!rT@#{EA#~F+XKB|Pr`F}xf3Cop$N^^=88K}rWMFz^P>J%9$VWhy?1Q;kbyg|cL z`cTXKqkIl$E||6#9-8B95tn?Yl1SXlr5|NsBr=8GKffjpi`Og=0?{{QKU z`HU))SXe-`#dL#wMv+Mz9Fu*Rw}9B1AVL}-Epe JTtM#B1OPk;WwHPO diff --git a/generateGroudTruth.py b/generateGroudTruth.py index 9a1959d..47b85e7 100644 --- a/generateGroudTruth.py +++ b/generateGroudTruth.py @@ -31,25 +31,65 @@ def generate_ground_truth_system(dx=2, du=1, dy=1, rng_seed=None): # --- 1. 生成稳定的 A 矩阵 (dx x dx) --- # 论文图 2(b) 显示了复特征值，遵循 6.2 节的极坐标法 - # 在极坐标下采样一对共轭复特征值 - r_squared = rng.uniform(0, 1.0) # 采样 r^2，确保在单位圆内 - r = np.sqrt(r_squared) - theta = rng.uniform(0, np.pi) # 仅在上半平面采样角度 + all_eigenvalues = [] - lambda1 = r * (np.cos(theta) + 1j * np.sin(theta)) - lambda2 = np.conjugate(lambda1) + # 确保至少一个实特征值，如果 dx 为奇数 + num_real_needed = dx % 2 + for _ in range(num_real_needed): + all_eigenvalues.append(rng.uniform(-0.99, 0.99)) # 稳定的实数特征值 + + # 生成复共轭特征值对 + num_complex_pairs_needed = (dx - num_real_needed) // 2 + for _ in range(num_complex_pairs_needed): + r_squared = rng.uniform(0, 1.0) + r = np.sqrt(r_squared) * rng.uniform(0.1, 0.99) # 确保半径在稳定范围内 + theta = rng.uniform(0, np.pi) # 仅在上半平面采样角度 + lambda1 = r * (np.cos(theta) + 1j * np.sin(theta)) + all_eigenvalues.append(lambda1) + all_eigenvalues.append(np.conjugate(lambda1)) - # 创建对应的实数块对角矩阵 - alpha, beta = lambda1.real, lambda1.imag - lambda_block = np.array([[alpha, beta], - [-beta, alpha]]) + rng.shuffle(all_eigenvalues) # 随机打乱特征值顺序 + + # 构造一个 Real Schur Form 的块对角矩阵 (J_block) + J_blocks = [] + current_eigenvalues = list(all_eigenvalues) # 创建一个可变副本 - # 生成一个随机正交矩阵 V (dx x dx) [cite: 511, 584] + while current_eigenvalues: + e = current_eigenvalues.pop(0) + if np.isreal(e): + J_blocks.append(np.array([[e.real]])) + else: + # 寻找其共轭 + conjugate_found = False + for i, ce in enumerate(current_eigenvalues): + if np.isclose(ce, np.conjugate(e)): # 使用 np.isclose 进行浮点比较 + alpha = e.real + beta = e.imag + J_blocks.append(np.array([[alpha, beta], [-beta, alpha]])) + current_eigenvalues.pop(i) # 移除已使用的共轭特征值 + conjugate_found = True + break + if not conjugate_found: + # 如果没有找到共轭（理论上不应发生，但作为健壮性处理），退化为实部 + J_blocks.append(np.array([[e.real]])) + + # 组装 J_block + if not J_blocks: # 处理 dx=0 或空列表的情况 + J_block = np.zeros((0,0)) + else: + # 使用 np.block 构造块对角矩阵 + J_block = np.block([ + [J_blocks[i] if i == j else np.zeros((J_blocks[i].shape[0], J_blocks[j].shape[1])) + for j in range(len(J_blocks))] + for i in range(len(J_blocks)) + ]) + + # 生成一个随机正交矩阵 V (dx x dx) Z = rng.standard_normal(size=(dx, dx)) V, _ = np.linalg.qr(Z) - # 组装 A = V * Lambda_block * V^T - A = V @ lambda_block @ V.T + # 组装 A = V * J_block * V^T + A = V @ J_block @ V.T # --- 2. 生成 B (dx x du) 和 C (dy x dx) 矩阵 --- # 元素从 N(0, 1) 独立采样 @@ -77,12 +117,14 @@ def generate_ground_truth_system(dx=2, du=1, dy=1, rng_seed=None): D = np.zeros((dy, du)) print(f"--- 成功生成 Ground Truth 系统 (尝试次数: {attempt + 1}) ---") - print(f"特征值: {lambda1:.4f}, {lambda2:.4f}") + + # 打印部分特征值，保持输出简洁 + print_eigs = [f"{e:.4f}" for e in all_eigenvalues[:min(dx, 4)]] + if dx > 4: + print_eigs.append("...") + print(f"特征值: {', '.join(print_eigs)}") + print(f"Gramian 能量占比: Wc={cond_c:.4f}, Wo={cond_o:.4f}") - print("A = \n", A) - print("B = \n", B) - print("C = \n", C) - print("D = \n", D) return A, B, C, D @@ -94,7 +136,171 @@ def generate_ground_truth_system(dx=2, du=1, dy=1, rng_seed=None): # 如果循环结束仍未成功 raise RuntimeError(f"在 {max_retries} 次尝试后未能生成一个良态的系统。") + +def transform_to_canonical_form(A, B, C, D): + """ + 将给定的 LTI 系统转换为可控规范形 (controllable canonical form)。 + + 此实现特指论文 Definition 3.1 中的 SISO 控制器规范形，并确保马尔可夫参数不变。 + + 参数: + A (np.ndarray): 状态矩阵 (dx x dx) + B (np.ndarray): 输入矩阵 (dx x du) + C (np.ndarray): 输出矩阵 (dy x dx) + D (np.ndarray): 直通矩阵 (dy x du) + + 返回: + tuple: 转换后的 (A_canon, B_canon, C_canon, D_canon) + + 抛出: + ValueError: 如果系统不是 SISO 或不可控。 + """ + dx = A.shape[0] # 状态维度 + du = B.shape[1] # 输入维度 + dy = C.shape[0] # 输出维度 + + # 1. 验证 SISO 系统，因为论文中描述的控制器规范形是针对 SISO 的 + if du != 1 or dy != 1: + raise ValueError("此转换仅针对单输入单输出 (SISO) 系统实现。") + + # 2. D 矩阵在状态空间变换中是不变的 + D_canon = D + + # 3. 计算 A 的特征多项式系数 + # np.poly(A) 返回 [1, p_{dx-1}, p_{dx-2}, ..., p_0] + # 其中特征多项式为 s^dx + p_{dx-1}s^{dx-1} + ... + p_1 s + p_0 + poly_coeffs_with_leading_one = np.poly(A) + # 论文 Definition 3.1 中 A_c 的最后一行是 [-a_0, -a_1, ..., -a_{dx-1}] + # 这里的 a_k 与 np.poly(A) 返回的 p_k 是一致的 (a_k = p_k)。 + # 需要将 np.poly(A)[1:] 的系数反序以匹配 A_c 的最后一行结构 [-p_0, -p_1, ..., -p_{dx-1}] + a_coeffs_for_Ac_row = -poly_coeffs_with_leading_one[1:][::-1] + + # 4. 构造 A_canon (控制器规范形) + A_canon = np.zeros((dx, dx)) + # 填充超对角线为 1 + for i in range(dx - 1): + A_canon[i, i+1] = 1 + # 填充最后一行 + A_canon[dx-1, :] = a_coeffs_for_Ac_row + + # 5. 构造 B_canon (控制器规范形) + B_canon = np.zeros((dx, du)) + B_canon[-1, 0] = 1 # SISO 情况下，最后一个元素为 1 + + # 6. 计算原始系统 (A, B) 的可控性矩阵 Q + # Q = [B, A@B, A^2@B, ..., A^(dx-1)@B] + Q = np.zeros((dx, dx)) + current_B_col = B + for i in range(dx): + Q[:, i] = current_B_col.flatten() # B 是 dx x 1 向量，需要展平 + current_B_col = A @ current_B_col + + # 7. 检查可控性 + if np.linalg.matrix_rank(Q) != dx: + raise ValueError("系统不可控，无法转换为可控规范形。") + + # 8. 计算规范形系统 (A_canon, B_canon) 的可控性矩阵 Q_canon + # Q_canon = [B_canon, A_canon@B_canon, ..., A_canon^(dx-1)@B_canon] + Q_canon = np.zeros((dx, dx)) + current_B_c_col = B_canon + for i in range(dx): + Q_canon[:, i] = current_B_c_col.flatten() + current_B_c_col = A_canon @ current_B_c_col + + # 9. 计算 C_canon + # 根据控制器规范形的定义，C_canon = [b_0, b_1, ..., b_{dx-1}] + # 其中 b_i 是传递函数分子系数。这可以通过将原始系统的马尔可夫参数 M_k + # 与规范系统的可控性矩阵 Q_canon 关联起来得到。 + # 我们知道 [M_1, M_2, ..., M_dx] = C_canon @ Q_canon + # 因此 C_canon = [M_1, M_2, ..., M_dx] @ Q_canon_inv + + # 计算原始系统的马尔可夫参数 M_k (k=1...dx) + # M_k = C @ A^(k-1) @ B + markov_params_row_matrix = np.zeros((1, dx)) + for k in range(dx): + # markov_params_row_matrix[0, k] 存储的是 M_{k+1}，对应 C A^k B + markov_params_row_matrix[0, k] = (C @ np.linalg.matrix_power(A, k) @ B).item() + + C_canon = markov_params_row_matrix @ np.linalg.inv(Q_canon) + + return A_canon, B_canon, C_canon, D_canon + if __name__ == '__main__': # --- 运行示例 --- # 设置一个随机种子，以便每次运行时都能得到相同的结果 - A_true, B_true, C_true, D_true = generate_ground_truth_system(rng_seed=42) + A_true, B_true, C_true, D_true = generate_ground_truth_system(dx=2, du=1, dy=1, rng_seed=42) + + print("\n--- 原始系统 ---") + print("A_true = \n", A_true) + print("B_true = \n", B_true) + print("C_true = \n", C_true) + print("D_true = \n", D_true) + + # 转换为可控规范形 + A_canon, B_canon, C_canon, D_canon = transform_to_canonical_form(A_true, B_true, C_true, D_true) + + print("\n--- 可控规范形系统 ---") + print("A_canon = \n", A_canon) + print("B_canon = \n", B_canon) + print("C_canon = \n", C_canon) + print("D_canon = \n", D_canon) + + # 验证转换结果 (不变的系统特性) + print("\n--- 验证转换结果 ---") + + # 1. D矩阵应相同 + print(f"D_true (原始): \n{D_true}") + print(f"D_canon (规范形): \n{D_canon}") + assert np.allclose(D_true, D_canon), "D 矩阵不匹配！" + print("D 矩阵匹配成功！") + + # 2. 马尔可夫参数应相同 (输入-输出行为不变) + # M_t = C A^(t-1) B (对于 t >= 1), M_0 = D + def get_markov_params(A, B, C, D, steps=5): + params = [D] # M_0 = D + for t in range(1, steps): + params.append(C @ np.linalg.matrix_power(A, t - 1) @ B) + return params + + markov_true = get_markov_params(A_true, B_true, C_true, D_true, steps=5) + markov_canon = get_markov_params(A_canon, B_canon, C_canon, D_canon, steps=5) + + print("\n马尔可夫参数验证 (前5步):") + for i in range(len(markov_true)): + print(f" Step {i}: True={markov_true[i].flatten()}, Canon={markov_canon[i].flatten()}") + assert np.allclose(markov_true[i], markov_canon[i]), f"马尔可夫参数在第 {i} 步不匹配！" + print("马尔可夫参数匹配成功！") + + # 3. 特征值应相同 + eig_true = np.linalg.eigvals(A_true) + eig_canon = np.linalg.eigvals(A_canon) + # 对特征值进行排序以进行比较，因为它们的顺序可能不同 + eig_true_sorted = np.sort(eig_true) + eig_canon_sorted = np.sort(eig_canon) + + print(f"\n特征值: True={eig_true_sorted}, Canon={eig_canon_sorted}") + assert np.allclose(eig_true_sorted, eig_canon_sorted), "特征值不匹配！" + print("特征值匹配成功！") + + print("\n所有验证成功！") + + # 尝试一个更高维度的系统 + print("\n--- 测试 dx=3 的系统 ---") + A_true_3, B_true_3, C_true_3, D_true_3 = generate_ground_truth_system(dx=3, du=1, dy=1, rng_seed=100) + A_canon_3, B_canon_3, C_canon_3, D_canon_3 = transform_to_canonical_form(A_true_3, B_true_3, C_true_3, D_true_3) + + print("\nA_canon (dx=3) = \n", A_canon_3) + print("B_canon (dx=3) = \n", B_canon_3) + print("C_canon (dx=3) = \n", C_canon_3) + + markov_true_3 = get_markov_params(A_true_3, B_true_3, C_true_3, D_true_3) + markov_canon_3 = get_markov_params(A_canon_3, B_canon_3, C_canon_3, D_canon_3) + for i in range(len(markov_true_3)): + assert np.allclose(markov_true_3[i], markov_canon_3[i]), f"dx=3 马尔可夫参数在第 {i} 步不匹配！" + print("dx=3 马尔可夫参数匹配成功！") + + eig_true_3 = np.linalg.eigvals(A_true_3) + eig_canon_3 = np.linalg.eigvals(A_canon_3) + assert np.allclose(np.sort(eig_true_3), np.sort(eig_canon_3)), "dx=3 特征值不匹配！" + print("dx=3 特征值匹配成功！") + diff --git a/temp.md b/temp.md new file mode 100644 index 0000000..5efac8d --- /dev/null +++ b/temp.md @@ -0,0 +1,173 @@ +好的，让我们从头开始，结合关键公式，详细讲解一下这篇文章。 + +### 1. 问题的起点：贝叶斯模型更新 (Bayesian Model Updating) + +**在科学和工程中，我们经常有一个数学模型 **$G$** 来描述一个物理系统，但模型中包含未知的参数 **$\theta$** **^1^^1^^1^^1^。我们通过实验测量得到一组数据 **$d$** ^2^。**贝叶斯更新** 的目的就是利用这些测量数据 **$d$** 来反推参数 **$\theta$** 的概率分布 ^3^^3^^3^^3^^3^^3^^3^^3^^3^。 + +这是通过贝叶斯定理实现的，它是本文所有工作的基础： + +$$ +p_{d}(\theta)=c_{E}^{-1}L(\theta)p_{0}(\theta) +$$ + +我们来拆解这个公式： + +* **$p_{d}(\theta)$：后验概率密度函数 (Posterior PDF)**。这是我们**想要的结果**，即在“看到”测量数据 $d$ 之后，参数 $\theta$ 的概率分布 。 +* **$p_{0}(\theta)$：先验概率密度函数 (Prior PDF)**。这是我们**开始的地方**，即在“看到”数据 *之前*，我们对 $\theta$ 的已有认知或假设。 +* **$L(\theta)$：似然函数 (Likelihood Function)**。它描述的是：**如果** 参数真的是 $\theta$，那么我们“看到”测量数据 $d$ 的概率有多大 [cite: 102]。 +* **$c_{E}$：模型证据 (Model Evidence)**。这是一个归一化常数，它保证 $p_{d}(\theta)$ 的总概率积分为1 [cite: 102]。它的计算公式为： +$$ + + +$$c\_{E}=\\int\_{\\mathbb{R}^{n}}L(\\theta)p\_{0}(\\theta)d\\theta +$$ + +$c_{E}$ 本身也很有用，它可以用来比较不同模型的好坏 [cite: 37, 105]。 + + + **挑战** **：在参数 **$\theta$** 维度很高（高维问题）时，这个后验分布 **$p_{d}(\theta)$** 非常复杂，几乎不可能直接求解，连 **$c_{E}$** 的积分也算不出来 **^4^^4^^4^^4^。 + +### 2. 求解方法：序列蒙特卡洛 (SMC) + +**为了解决这个难题，论文使用了 ****序列蒙特卡洛 (SMC)** 方法 ^5^。 + +**SMC 的核心思想是“逐步逼近”：它不试图一步到位从 **$p_{0}(\theta)$**（先验）跳到 **$p_{d}(\theta)$**（后验），而是构建一系列平滑过渡的中间分布 **$f_{l}(\theta)$** **^6^^6^^6^^6^。 + +最常用的方法是“退火”或“回火”（tempering），通过一个指数 **$\beta_l$** 来实现： + +$$ +f_{l}(\theta)\propto L(\theta)^{\beta_{l}}p_{0}(\theta), \quad \text{其中 } 0 = \beta_{0} < \beta_{1} < \dots < \beta_{L} = 1 +$$[cite\_start][cite: 174-175] + +本文重排并清理了数学公式与格式，确保在常见 Markdown 渲染器中正确显示（提示：需要支持数学渲染的环境，如 VS Code 预览的数学渲染或 GitHub 的内置 KaTeX/MathJax）。 + +### 1. 问题的起点：贝叶斯模型更新 (Bayesian Model Updating) + +在科学和工程中，我们经常有一个数学模型 $G$ 来描述物理系统，但模型中包含未知参数 $\theta$。我们通过实验测量得到一组数据 $d$。贝叶斯更新的目标是利用测量数据 $d$ 反推参数 $\theta$ 的后验分布。 + +这是通过贝叶斯定理实现的： + +$$ +p_{d}(\theta) = c_{E}^{-1} \, L(\theta) \, p_{0}(\theta) +$$ + +其中 $c_E$（模型证据）保证后验归一： + +$$ +c_{E} = \int_{\mathbb{R}^{n}} L(\theta) \, p_{0}(\theta) \, d\theta. +$$ + +当参数维度很高时，$p_d(\theta)$ 的形状可能非常复杂，直接计算甚至连证据 $c_E$ 的积分都很难。 + +### 2. 求解方法：序列蒙特卡洛 (SMC) + +SMC 的核心思想是“逐步逼近”：不从先验 $p_0(\theta)$ 一步跳到后验 $p_d(\theta)$，而是构建一系列中间分布 $f_l(\theta)$，常见做法是退火/回火（tempering）： + +$$ +f_{l}(\theta) \propto L(\theta)^{\beta_{l}} \, p_{0}(\theta), \quad 0=\beta_0 < \beta_1 < \cdots < \beta_L = 1. +$$ + +- 当 $\beta_0 = 0$ 时，$f_0(\theta) \propto p_0(\theta)$（先验）。 +- 当 $\beta_L = 1$ 时，$f_L(\theta) \propto L(\theta) p_0(\theta)$（后验）。 + +从 $f_{l-1}$ 到 $f_l$，SMC 通常包含三步： + +1) 重加权（Reweighting）：根据新的 $\beta_l$ 重新计算粒子权重； +2) 重采样（Resampling）：复制高权重、淘汰低权重粒子（会导致样本贫化）； +3) 移动（Moving）：在保持 $f_l$ 不变的前提下，对粒子做若干步 MCMC 以恢复多样性。 + +### 3. “移动”步骤的经典算法：pCN + +pCN（预条件 Crank–Nicolson）在高维问题中表现稳健，尤其适用于高斯先验 $p_0(\theta)=\mathcal{N}(0,I)$ 的情形。其提议为： + +$$ +v = \sqrt{1-s^{2}}\, \theta_{0} + s\, \xi, \quad \xi \sim \mathcal{N}(0, I), \; s\in[0,1]. +$$ + +对应的 Metropolis–Hastings 接受率在退火分布 $f_l \propto L^{\beta_l} p_0$ 下可简化为： + +$$ +\alpha(\theta_{0}, v) = \min\left\{1, \frac{L(v)^{\beta_{l}}}{L(\theta_{0})^{\beta_{l}}} \right\}. +$$ + +关键点：先验项在接受率中相互抵消，接受与否仅由似然的相对变化决定，这使 pCN 的性能对参数维度不敏感。 + +然而，pCN 的提议协方差等同于先验（单位阵），在后期 $f_l$ 已经很“窄”且相关性强时，方向上不自适应，效率会下降。 + +--- + +### 4. 新算法一：cov-pCN（协方差信息 pCN） + +思想：在第 $l$ 步，利用带权样本估计目标分布 $f_l$ 的协方差 $\hat{\Sigma}_{f_l}$，并让提议分布适应该协方差结构。 + +令 $\hat{\Sigma}_{f_l} = W L W^{\top}$ 为特征分解，则广义 pCN 提议为： + +$$ +v = A\,\theta_0 + s\, W L^{1/2} \xi, \quad \xi \sim \mathcal{N}(0,I), +$$ + +其中 + +$$ +A = \sqrt{I - s^{2} \hat{\Sigma}_{f_l}} = W\, \sqrt{I - s^{2} L} \, W^{\top}. +$$ + +这样构造可保持接受率与 pCN 同型： + +$$ +\alpha(\theta_{0}, v) = \min\left\{1, \frac{L(v)^{\beta_{l}}}{L(\theta_{0})^{\beta_{l}}} \right\}, +$$ + +但提议方向对齐于 $f_l$ 的主协方差方向，混合更高效。 + +--- + +### 5. 新算法二：pc-M（主成分 M–H） + +思想：高维协方差的主要方差集中在少数主方向上。仅沿主成分方向做随机游走即可。 + +做特征分解 $\hat{\Sigma}_{f_l} = W L W^{\top}$，取前 $n_r$ 个特征对 $(\rho_k, c_k)$，按权重（如与 $\rho_k$ 相关）随机选方向 $k$，然后： + +$$ +v = \theta_0 + s\, \rho_k\, c_k\, \xi, \quad \xi \sim \mathcal{N}(0,1). +$$ + +这是普通随机游走（RWM），接受率为： + +$$ +\alpha(\theta_0, v) = \min\left\{1, \frac{f_l(v)}{f_l(\theta_0)} \right\} += \min\left\{1, \frac{L(v)^{\beta_l} \, p_0(v)}{L(\theta_0)^{\beta_l} \, p_0(\theta_0)} \right\}. +$$ + +--- + +### 6. 协方差矩阵的估计与递归更新 + +两种新算法都依赖于“好的”协方差估计 $\hat{\Sigma}_{f_l}$。 + +1) 初始估计：用重采样前的带权样本估计均值 $\hat{\mu}_{f_l,0}$ 和协方差 + +$$ +\hat{\Sigma}_{f_l,0} = \sum_{j=1}^{J} w_j\, (\theta_j - \hat{\mu}_{f_l,0})(\theta_j - \hat{\mu}_{f_l,0})^{\top}. +$$ + +2) 递归更新：每获得 $J_a$ 个新样本就批量更新一次，降低特征分解开销（$O(n^3)$）。令步长为 $\gamma_{\mathrm{iter}}$，则 + +$$ +\hat{\Sigma}_{f_l,\mathrm{iter}} = \hat{\Sigma}_{f_l,\mathrm{iter-1}} \\ +\quad + \; \gamma_{\mathrm{iter}} \Bigg[ \frac{1}{J_a} \sum_{j=1}^{J_a} +(\theta_j - \hat{\mu}_{f_l,\mathrm{iter-1}})(\theta_j - \hat{\mu}_{f_l,\mathrm{iter-1}})^{\top} +\; - \hat{\Sigma}_{f_l,\mathrm{iter-1}} \Bigg]. +$$ + +这是标准的递归平均思想，用 $\gamma$ 平衡新旧信息。 + +### 7. 实验结论（摘要） + +- 高维（如 236 维水文层析）问题中，需要较多退火步数（如约 40）才能从先验到后验； +- 标准 pCN 在较早阶段接受率就显著下降，移动步骤失效； +- cov-pCN 的接受率下降更慢，多数退火步骤中仍能有效“移动”粒子； +- 在难问题上，cov-pCN 的均值误差、失配与不确定性更优。 + +总之：在 SMC 中引入自适应协方差的 cov-pCN，通常比标准 pCN 更稳健高效，尤其是需要大量退火步骤的复杂高维反问题。 + $$