程序仍有问题
This commit is contained in:
+510
File diff suppressed because one or more lines are too long
@@ -17,29 +17,29 @@ $$
|
|||||||
### k <-> k+1
|
### k <-> k+1
|
||||||
|
|
||||||
1. Birth move
|
1. Birth move
|
||||||
$$
|
辅助变量的先验：$u~U\left( -1,1 \right) $
|
||||||
\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|
|
提议分布:$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
|
2. Death move
|
||||||
$$
|
辅助变量的先验：$\frac{1}{k+1} $
|
||||||
\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|
|
提议分布:$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
|
### k <-> k+2
|
||||||
|
|
||||||
1. Birth move
|
1. Birth move
|
||||||
$$
|
映射关系: $\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) $
|
||||||
\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|
|
辅助变量的先验: $\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
|
2. Death move
|
||||||
$$
|
辅助变量的先验: $\frac{1}{k+1}$
|
||||||
\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|
|
提议分布:$ 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 $$
|
||||||
$其中 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|$
|
|
||||||
|
|
||||||
## 同维度游动
|
## 同维度游动
|
||||||
|
|
||||||
@@ -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
|
r_1=a+b\text{、}r_2=a-b
|
||||||
$$
|
$$
|
||||||
|
|
||||||
|
|
||||||
再确定雅可比矩阵行列式
|
再确定雅可比矩阵行列式
|
||||||
|
|
||||||
$$
|
$$
|
||||||
J_{C\rightarrow R}=\left| \det \left( \begin{matrix}
|
J_{C\rightarrow R}=\left| \det \left( \begin{matrix}
|
||||||
\frac{\partial a}{\partial r_1}& \frac{\partial a}{\partial r_2}\\
|
\frac{\partial a}{\partial r_1}& \frac{\partial a}{\partial r_2}\\
|
||||||
|
|||||||
+468
-152
@@ -48,8 +48,10 @@ class RJMCMC_Sampler:
|
|||||||
self.lambda_b_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)
|
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):
|
def _update_current_eigenvalues(self):
|
||||||
@@ -59,6 +61,16 @@ class RJMCMC_Sampler:
|
|||||||
self.current_eigenvalues = np.roots(poly_coeffs)
|
self.current_eigenvalues = np.roots(poly_coeffs)
|
||||||
else:
|
else:
|
||||||
self.current_eigenvalues = np.array([])
|
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):
|
def _cal_current_eigenvalues(self, a_coeffs_temp):
|
||||||
"""一个辅助函数，根据给定的 a_coeffs 计算对应的特征值。"""
|
"""一个辅助函数，根据给定的 a_coeffs 计算对应的特征值。"""
|
||||||
@@ -206,7 +218,7 @@ class RJMCMC_Sampler:
|
|||||||
|
|
||||||
# --- 计算似然 ---
|
# --- 计算似然 ---
|
||||||
# 预测误差 nu_t
|
# 预测误差 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
|
nu_t = y_t - y_pred
|
||||||
|
|
||||||
# 预测误差协方差 S_t
|
# 预测误差协方差 S_t
|
||||||
@@ -228,7 +240,7 @@ class RJMCMC_Sampler:
|
|||||||
P_update = I_KC @ P_pred @ I_KC.T + K_t @ Gamma @ K_t.T
|
P_update = I_KC @ P_pred @ I_KC.T + K_t @ Gamma @ K_t.T
|
||||||
|
|
||||||
# --- 为下一次循环准备预测 (t+1) ---
|
# --- 为下一次循环准备预测 (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
|
P_pred = A @ P_update @ A.T + Sigma
|
||||||
|
|
||||||
return total_log_likelihood # 返回标量值
|
return total_log_likelihood # 返回标量值
|
||||||
@@ -241,6 +253,14 @@ class RJMCMC_Sampler:
|
|||||||
k = self.current_k
|
k = self.current_k
|
||||||
current_eigs = self.current_eigenvalues
|
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
|
can_add_complex = (k + 2) <= self.k_max
|
||||||
add_real = True
|
add_real = True
|
||||||
@@ -252,33 +272,51 @@ class RJMCMC_Sampler:
|
|||||||
if add_real:
|
if add_real:
|
||||||
k_new = k + 1
|
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)
|
# 从提议分布 q(u) 中采样辅助变量 u = (u_lambda, u_b)
|
||||||
u_lambda = np.random.uniform(-1.0, 1.0) # 新特征值
|
u_lambda = np.random.uniform(-1.0, 1.0) # 新特征值
|
||||||
u_b = np.random.normal(0, 1) # 新 b 系数
|
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|
|
# 计算对数雅可比行列式 log|J|
|
||||||
if k == 0:
|
if k == 0:
|
||||||
log_det_jacobian = 0.0
|
log_det_jacobian = 0.0
|
||||||
else:
|
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_eigs = np.append(current_eigs, u_lambda)
|
||||||
new_a = np.real(np.poly(new_eigs)[1:][::-1])
|
new_a = np.real(np.poly(new_eigs)[1:][::-1])
|
||||||
new_b = np.append(self.current_b, u_b)
|
new_b = np.append(self.current_b, u_b)
|
||||||
|
|
||||||
# 接受率中的项为 |J| / q(u)，在对数空间中为 log|J| - log(q(u))
|
# 接受率中的项为 q_backward(u) * |J| / q_forward(u)
|
||||||
log_ratio = log_det_jacobian - log_q_forward
|
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:
|
else:
|
||||||
# --- 诞生一对复共轭根 (k -> k+2) ---
|
# --- 诞生一对复共轭根 (k -> k+2) ---
|
||||||
k_new = 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)
|
# 采样辅助变量 u = (rho, theta, u_b1, u_b2)
|
||||||
rho = np.sqrt(np.random.uniform(0, 1.0))
|
rho = np.sqrt(np.random.uniform(0, 1.0))
|
||||||
theta = np.random.uniform(0, np.pi)
|
theta = np.random.uniform(0, np.pi)
|
||||||
@@ -286,25 +324,31 @@ class RJMCMC_Sampler:
|
|||||||
u_b1, u_b2 = np.random.normal(0, 1, 2)
|
u_b1, u_b2 = np.random.normal(0, 1, 2)
|
||||||
|
|
||||||
# 计算对数提议密度 log(q(u))
|
# 计算对数提议密度 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|
|
# 计算对数雅可比行列式 log|J|
|
||||||
# |J| = |product(|u_lambda - lambda_i|^2) * (2*Im(u_lambda))|
|
# |J| = 2*ρ
|
||||||
if k == 0:
|
log_jacobian = np.log(2 * rho)
|
||||||
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_eigs = np.append(current_eigs, [u_lambda, np.conjugate(u_lambda)])
|
||||||
new_a = np.real(np.poly(new_eigs)[1:][::-1])
|
new_a = np.real(np.poly(new_eigs)[1:][::-1])
|
||||||
new_b = np.append(self.current_b, [u_b1, u_b2])
|
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):
|
def _propose_death(self):
|
||||||
"""
|
"""
|
||||||
提议一个 "消亡" 转移 (k -> k-1 或 k -> k-2)。
|
提议一个 "消亡" 转移 (k -> k-1 或 k -> k-2)。
|
||||||
@@ -317,8 +361,8 @@ class RJMCMC_Sampler:
|
|||||||
real_eigs_indices = np.where(np.isreal(current_eigs))
|
real_eigs_indices = np.where(np.isreal(current_eigs))
|
||||||
complex_eigs_indices = np.where((np.iscomplex(current_eigs)) & (current_eigs.imag > 0)) #type: ignore
|
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_real = len(real_eigs_indices[0]) > 0
|
||||||
can_remove_complex = len(complex_eigs_indices) > 0
|
can_remove_complex = len(complex_eigs_indices[0]) > 0
|
||||||
|
|
||||||
if not can_remove_real and not can_remove_complex:
|
if not can_remove_real and not can_remove_complex:
|
||||||
return None # 无法执行消亡
|
return None # 无法执行消亡
|
||||||
@@ -335,66 +379,129 @@ class RJMCMC_Sampler:
|
|||||||
# --- 消亡一个实数根 (k -> k-1) ---
|
# --- 消亡一个实数根 (k -> k-1) ---
|
||||||
k_new = k - 1
|
k_new = k - 1
|
||||||
|
|
||||||
|
# 生成与消亡实根对应的诞生过程的对数概率之比
|
||||||
|
log_death_real_forward = np.log(0.5)
|
||||||
|
log_death_real_backward = np.log(0.5)
|
||||||
|
|
||||||
# 1. 随机选择一个实数根移除
|
# 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]
|
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
|
# 移除的根和b系数构成了逆向（诞生）提议的辅助变量 u
|
||||||
u_lambda = lambda_removed
|
u_lambda = lambda_removed
|
||||||
u_b = self.current_b[-1]
|
u_b = b_removed
|
||||||
|
|
||||||
# 2. 计算逆向提议的对数密度 log(q(u))
|
# 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|
|
# 3. 计算对应诞生过程的对数雅可比行列式 log|J|
|
||||||
remaining_eigs = np.delete(current_eigs, idx_to_remove)
|
remaining_eigs = np.delete(current_eigs, idx_to_remove)
|
||||||
if k_new == 0:
|
if k_new == 0:
|
||||||
log_jacobian_birth = 0.0
|
log_jacobian_birth = 0.0
|
||||||
else:
|
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])
|
if k_new == 0:
|
||||||
new_b = self.current_b[:-1]
|
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|
|
# 接受率中的项为 q_backward(u) / (|J| * q_forward(u))
|
||||||
log_ratio = log_q_reverse - log_jacobian_birth
|
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:
|
else:
|
||||||
# --- 消亡一对复共轭根 (k -> k-2) ---
|
# --- 消亡一对复共轭根 (k -> k-2) ---
|
||||||
k_new = 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]
|
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
|
||||||
u_lambda = lambda_removed
|
u_lambda = lambda_removed
|
||||||
u_b1, u_b2 = self.current_b[-2:]
|
u_b1, u_b2 = self.current_b[-2:]
|
||||||
|
|
||||||
# 计算逆向提议的对数密度 log(q(u))
|
# 计算逆向提议的对数密度 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])
|
remaining_eigs = np.delete(current_eigs, [complex_idx_to_remove, conjugate_idx])
|
||||||
if k_new == 0:
|
log_jacobian = np.log(2 * rho)
|
||||||
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])
|
if k_new == 0:
|
||||||
new_b = self.current_b[:-2]
|
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
|
log_ratio = -log_q_forward - log_jacobian + log_q_backward
|
||||||
|
|
||||||
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_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):
|
def _log_prior_b(self, b_coeffs):
|
||||||
"""
|
"""
|
||||||
@@ -404,6 +511,7 @@ class RJMCMC_Sampler:
|
|||||||
"""
|
"""
|
||||||
#
|
#
|
||||||
return stats.norm.logpdf(b_coeffs, 0, 1).sum()
|
return stats.norm.logpdf(b_coeffs, 0, 1).sum()
|
||||||
|
|
||||||
|
|
||||||
def _generate_stable_a(self, k):
|
def _generate_stable_a(self, k):
|
||||||
"""
|
"""
|
||||||
@@ -513,136 +621,344 @@ class RJMCMC_Sampler:
|
|||||||
|
|
||||||
return A_c, B_c, C_c, D_c
|
return A_c, B_c, C_c, D_c
|
||||||
|
|
||||||
def _log_prior_full(self, k, a_coeffs, b_coeffs, sigma2_w, sigma2_z):
|
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 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)
|
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)
|
if k > 0:
|
||||||
a_coeffs = self._cal_current_a_coeffs(eigenvalues)
|
log_prior_eigen = self._log_prior_eigenvalues(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:
|
else:
|
||||||
raise ValueError("Either a_coeffs or eigenvalues must be provided.")
|
log_prior_eigen = 0.0
|
||||||
|
|
||||||
if log_prior_a == -np.inf: return -np.inf
|
|
||||||
|
|
||||||
log_prior_b = self._log_prior_b(b_coeffs)
|
log_prior_b = self._log_prior_b(b_coeffs)
|
||||||
precision_w = 1.0 / sigma2_w
|
precision_w = 1.0 / sigma2_w
|
||||||
precision_z = 1.0 / sigma2_z
|
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_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_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
|
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
|
||||||
# 2. 计算似然
|
|
||||||
log_likelihood = self._log_likelihood(k, a_coeffs, b_coeffs, sigma2_w, sigma2_z)
|
|
||||||
|
|
||||||
return log_prior + 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
|
k = self.current_k
|
||||||
current_eigs = np.copy(self.current_eigenvalues)
|
|
||||||
|
|
||||||
# 随机选择两个不同的实数根
|
# 随机选择两个不同的实数根(可能是虚部为0的复数)，并删除他们
|
||||||
real_indices = np.where(np.isclose(current_eigs.imag, 0))
|
real_indices = np.where(np.isreal(current_eigs))
|
||||||
idx1, idx2 = np.random.choice(real_indices, 2, replace=False)
|
if len(real_indices[0]) < 2:
|
||||||
lambda1, lambda2 = current_eigs[idx1].real, current_eigs[idx2].real
|
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
|
||||||
|
|
||||||
# 确定性映射 -> 新的复共轭对
|
current_eigs = np.delete(current_eigs, [idx1, idx2])
|
||||||
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])
|
new_theta = np.random.uniform(0, np.pi)
|
||||||
proposal_eigs = np.append(proposal_eigs, new_complex_pair)
|
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)])
|
||||||
|
|
||||||
# 检查稳定性
|
# 对b进行游动
|
||||||
if np.any(np.abs(proposal_eigs) >= 1.0): return
|
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_acceptance_ratio = (log_post_unnormalized_proposal - log_post_unnormalized_current) + log_proposal_ratio
|
||||||
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. 计算提议比和雅可比项
|
acceptance_ratio = np.exp(log_acceptance_ratio)
|
||||||
# 正向 (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
|
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}
|
||||||
|
|
||||||
|
|
||||||
# 接受或拒绝
|
def _split_eigenvalues(self, current_eigs):
|
||||||
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
|
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))
|
complex_indices = np.where((np.iscomplex(current_eigs)) & (current_eigs.imag > 0)) #type: ignore
|
||||||
idx_c = np.random.choice(complex_indices)
|
if len(complex_indices[0]) < 1:
|
||||||
lambda_c = current_eigs[idx_c]
|
return # 没有复共轭对，无法分裂
|
||||||
a, b = lambda_c.real, lambda_c.imag
|
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)
|
current_eigs = np.delete(current_eigs, [idx_to_remove, conjugate_idx])
|
||||||
lambda1 = a + b * u
|
# 生成两个实数根
|
||||||
lambda2 = a - b * u
|
real_eig1 = np.random.uniform(-1.0, 1.0)
|
||||||
new_real_pair = [lambda1, lambda2]
|
real_eig2 = np.random.uniform(-1.0, 1.0)
|
||||||
|
current_eigs = np.append(current_eigs, [real_eig1, real_eig2])
|
||||||
|
|
||||||
# 构造提议的特征值集合
|
# 对b进行游动
|
||||||
idx_c_conj = np.where(np.isclose(current_eigs, np.conjugate(lambda_c)))
|
proposal_b, log_prior_b_current, log_prior_b_proposal = self._purpose_b_swimming()
|
||||||
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
|
# 正向过程是任选一个复共轭对，分裂为两个实根
|
||||||
|
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_acceptance_ratio = (log_post_unnormalized_proposal - log_post_unnormalized_current) + log_proposal_ratio
|
||||||
log_post_current = self._log_posterior(k, self.current_b, self.current_sigma2_w, self.current_sigma2_z, self.current_a, eigenvalues=current_eigs)
|
acceptance_ratio = np.exp(log_acceptance_ratio)
|
||||||
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)
|
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)
|
|
||||||
|
|||||||
+225
-19
@@ -31,25 +31,65 @@ def generate_ground_truth_system(dx=2, du=1, dy=1, rng_seed=None):
|
|||||||
# --- 1. 生成稳定的 A 矩阵 (dx x dx) ---
|
# --- 1. 生成稳定的 A 矩阵 (dx x dx) ---
|
||||||
# 论文图 2(b) 显示了复特征值，遵循 6.2 节的极坐标法
|
# 论文图 2(b) 显示了复特征值，遵循 6.2 节的极坐标法
|
||||||
|
|
||||||
# 在极坐标下采样一对共轭复特征值
|
all_eigenvalues = []
|
||||||
r_squared = rng.uniform(0, 1.0) # 采样 r^2，确保在单位圆内
|
|
||||||
r = np.sqrt(r_squared)
|
|
||||||
theta = rng.uniform(0, np.pi) # 仅在上半平面采样角度
|
|
||||||
|
|
||||||
lambda1 = r * (np.cos(theta) + 1j * np.sin(theta))
|
# 确保至少一个实特征值，如果 dx 为奇数
|
||||||
lambda2 = np.conjugate(lambda1)
|
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))
|
||||||
|
|
||||||
# 创建对应的实数块对角矩阵
|
rng.shuffle(all_eigenvalues) # 随机打乱特征值顺序
|
||||||
alpha, beta = lambda1.real, lambda1.imag
|
|
||||||
lambda_block = np.array([[alpha, beta],
|
# 构造一个 Real Schur Form 的块对角矩阵 (J_block)
|
||||||
[-beta, alpha]])
|
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))
|
Z = rng.standard_normal(size=(dx, dx))
|
||||||
V, _ = np.linalg.qr(Z)
|
V, _ = np.linalg.qr(Z)
|
||||||
|
|
||||||
# 组装 A = V * Lambda_block * V^T
|
# 组装 A = V * J_block * V^T
|
||||||
A = V @ lambda_block @ V.T
|
A = V @ J_block @ V.T
|
||||||
|
|
||||||
# --- 2. 生成 B (dx x du) 和 C (dy x dx) 矩阵 ---
|
# --- 2. 生成 B (dx x du) 和 C (dy x dx) 矩阵 ---
|
||||||
# 元素从 N(0, 1) 独立采样
|
# 元素从 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))
|
D = np.zeros((dy, du))
|
||||||
|
|
||||||
print(f"--- 成功生成 Ground Truth 系统 (尝试次数: {attempt + 1}) ---")
|
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(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
|
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} 次尝试后未能生成一个良态的系统。")
|
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__':
|
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 特征值匹配成功！")
|
||||||
|
|
||||||
|
|||||||
@@ -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 更稳健高效，尤其是需要大量退火步骤的复杂高维反问题。
|
||||||
|
$$
|
||||||
Reference in New Issue
Block a user