import numpy as np
from scipy.linalg import toeplitz, solve_toeplitz
def ar_gammas(phis, sigma2, max_lag=20):
"""
计算AR(n)过程的自协方差函数(0到max_lag)
Parameters:
phis (list): 自回归系数[𝜙_1, 𝜙_2, ..., 𝜙_n]
sigma2 (float): 噪声方差𝜎_𝜀^2
max_lag (int): 最大滞后阶数
Returns:
list: γ(0)到γ(max_lag)的自协方差值
"""
n = len(phis)
if n == 0: # 白噪声过程
return [sigma2] + [0] * max_lag
# 构建Yule-Walker方程的托普利茨矩阵
r = np.zeros(n)
r[0] = 1
for i in range(1, n):
r[i] = -sum(phis[j] * (r[i - 1 - j] if i - 1 - j < len(r) else 0) for j in range(min(i, n)))
# 构建右侧向量 [σ_ε^2, 0, 0, ...]
b = np.zeros(n)
b[0] = -sigma2
# 解Yule-Walker方程 (使用Levinson-Durbin算法优化)
try:
phi_est = solve_toeplitz((r, r), b)
except np.linalg.LinAlgError:
# 托普利茨解法失败时使用通用解法
R = toeplitz(r)
phi_est = np.linalg.solve(R, b)
# 计算γ(0)
if abs(sigma2) > 1e-10: # 检查sigma2是否为零
inner_term = -b / sigma2
denominator = 1 + np.dot(phi_est, inner_term)
# 检查分母是否为零
if abs(denominator) > 1e-10:
gamma0 = sigma2 / denominator
else:
gamma0 = np.nan # 分母为零时设为NaN
else:
gamma0 = np.nan # sigma2为零时设为NaN
# 计算初始自协方差γ(0)到γ(n)
gammas = [gamma0]
for k in range(1, n + 1):
gamma_k = sum(phis[j - 1] * gammas[k - j]
for j in range(1, min(k, n) + 1)
if k - j < len(gammas))
gammas.append(gamma_k)
# 递推高阶自协方差(k > n)
for k in range(n + 1, max_lag + 1):
gamma_k = sum(phis[j - 1] * gammas[k - j]
for j in range(1, n + 1)
if k - j < len(gammas))
gammas.append(gamma_k)
return gammas[:max_lag + 1]
def calculate_ar_properties(phi0, phis, sigma2, max_lag=20):
"""
计算AR过程的完整统计量
Parameters:
phi0 (float): 常数项𝜙_0
phis (list): 自回归系数[𝜙_1, 𝜙_2, ..., 𝜙_n]
sigma2 (float): 噪声方差𝜎_𝜀^2
max_lag (int): 最大滞后阶数
Returns:
dict: 包含期望、方差、自协方差和自相关系数的字典
"""
# 计算数学期望𝜇
mu = phi0 / (1 - sum(phis)) if len(phis) > 0 else phi0
# 计算自协方差函数
gammas = ar_gammas(phis, sigma2, max_lag)
variance = gammas[0] # γ(0)即方差
# 计算自相关系数𝜌(k)
rhos = [gamma_k / variance for gamma_k in gammas]
return {
'expectation': mu,
'variance': variance,
'autocovariance': gammas,
'autocorrelation': rhos
}
def calculate_mixed_process(props1, props2, p, max_lag=20):
"""
计算混合过程Y_t的统计量
Parameters:
props1 (dict): X_t1的统计量
props2 (dict): X_t2的统计量
p (float): 伯努利参数
max_lag (int): 最大滞后阶数
Returns:
dict: Y_t的期望、方差和自协方差
"""
mu1, mu2 = props1['expectation'], props2['expectation']
gamma1, gamma2 = props1['autocovariance'], props2['autocovariance']
# 混合过程的期望
mu_Y = p * mu1 + (1 - p) * mu2
# 混合过程的方差 (k=0)
var_Y = (p * gamma1[0] + (1 - p) * gamma2[0] +
p * (1 - p) * (mu1 - mu2) ** 2)
# 混合过程的自协方差函数
gamma_Y = [var_Y] # k=0
for k in range(1, max_lag + 1):
cov_k = p ** 2 * gamma1[k] + (1 - p) ** 2 * gamma2[k]
gamma_Y.append(cov_k)
return {
'expectation': mu_Y,
'variance': var_Y,
'autocovariance': gamma_Y
}
# ====================== 使用示例 ======================
if __name__ == "__main__":
# 示例参数 (AR(20)过程)
phi0_1 = 0.5
phis_1 = [0.3, 0.2, 0.3, 0.3, 0.4, 0.6, 0.4, 0.3, 0.2, 0.3, 0.3, 0.4, 0.6, 0.4, 0.5, 0.1, 0.4, 0.6, 0.7, 0.4] # 𝜙1=0.6, 𝜙2=-0.2
sigma2_1 = 1
# AR(20)过程
phi0_2 = 0.5
phis_2 = [0.1, 0.3, 0.2, 0.5, 0.7, 0.5, 0.3, 0.4, -0.5, -0.2, -0.1, 0.7, 0.5, 0.3, 0.4, -0.5, -0.2, -0.1, 0.7, 0.5] # 𝜙1=0.7
sigma2_2 = 1.0
# 伯努利混合参数
p = 0.6
# 计算两个AR过程的统计量
ar_props1 = calculate_ar_properties(phi0_1, phis_1, sigma2_1, max_lag=7)
ar_props2 = calculate_ar_properties(phi0_2, phis_2, sigma2_2, max_lag=7)
# 计算混合过程统计量
mixed_props = calculate_mixed_process(ar_props1, ar_props2, p, max_lag=7)
# 打印结果
print("AR过程1 (X_t1):")
print(f"期望 μ1 = {ar_props1['expectation']:.4f}")
print(f"方差 γ1(0) = {ar_props1['variance']:.4f}")
print("自相关系数 𝜌1(k):")
for k, rho in enumerate(ar_props1['autocorrelation'][:20]):
print(f"𝜌1({k}) = {rho:.4f}")
print("\nAR过程2 (X_t2):")
print(f"期望 μ2 = {ar_props2['expectation']:.4f}")
print(f"方差 γ2(0) = {ar_props2['variance']:.4f}")
print("自相关系数 𝜌2(k):")
for k, rho in enumerate(ar_props2['autocorrelation'][:20]):
print(f"𝜌2({k}) = {rho:.4f}")
print("\n混合过程 (Y_t):")
print(f"期望 E[Y_t] = {mixed_props['expectation']:.4f}")
print(f"方差 Var(Y_t) = {mixed_props['variance']:.4f}")
print("自协方差函数 𝛾_Y(k):")
for k, gamma in enumerate(mixed_props['autocovariance'][:20]):
print(f"𝛾_Y({k}) = {gamma:.4f}")
上述代码是原始代码，请修改原始代码，使其能计算并输出当 h=20时的的自相关系数 ρ_i(h)和自协方差函数 γ_i(h),也就是计算0到20的自相关系数 ρ_i(h)和自协方差函数 γ_i(h),请在原始代码中修改，谢谢
最新发布
本文介绍了如何在TensorFlow中调整输出格式，去除e表示并设置小数点后仅显示四位，以提高输出的可读性。通过np.set_printoptions方法实现，适合日常开发中的美化输出。
扫一扫
1017
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
