互信息法计算延迟时间及其在高维相空间重构中的应用

互信息法是计算时间序列延迟时间的重要方法，用于构建高维相空间矩阵（相空间重构）。

互信息法原理

互信息法基于信息论，通过计算原始时间序列与其延迟版本之间的信息相关性来确定最优延迟时间τ。

核心公式

互信息I(τ)定义为：
$$I(\tau )=\sum _{x(t),x(t+\tau )}p(x(t),x(t+\tau )){\log }_2\frac{p(x(t),x(t+\tau ))}{p(x(t))p(x(t+\tau ))}$$

其中：

• $p(x(t))$是时间序列在t时刻的概率分布
• $p(x(t+\tau ))$是延迟τ后的概率分布
• $p(x(t),x(t+\tau ))$是联合概率分布

算法流程

1. 对时间序列进行离散化（分箱）
2. 计算原始序列的概率分布
3. 对每个τ值计算延迟序列的概率分布
4. 计算联合概率分布
5. 计算互信息值I(τ)
6. 选择互信息函数第一个极小值对应的τ作为最优延迟时间

Mt代码

function [tau_opt, mutual_info] = mutual_information_delay(x, max_tau, n_bins)
    % 互信息法计算最优延迟时间
    % 输入:
    %   x - 时间序列 (1xN)
    %   max_tau - 最大延迟时间
    %   n_bins - 直方图分箱数
    % 输出:
    %   tau_opt - 最优延迟时间
    %   mutual_info - 互信息值向量
    
    N = length(x);
    mutual_info = zeros(1, max_tau);
    
    % 计算原始序列的直方图
    [counts, edges] = histcounts(x, n_bins);
    p_x = counts / N;
    
    % 计算不同τ值的互信息
    for tau = 1:max_tau
        % 创建延迟序列
        x_delayed = x(1+tau:end);
        x_original = x(1:end-tau);
        M = length(x_original);
        
        % 计算延迟序列的直方图
        [counts_delayed, ~] = histcounts(x_delayed, edges);
        p_y = counts_delayed / M;
        
        % 计算联合直方图
        joint_counts = histcounts2(x_original, x_delayed, edges, edges);
        p_xy = joint_counts / M;
        
        % 计算互信息
        I = 0;
        for i = 1:n_bins
            for j = 1:n_bins
                if p_xy(i,j) > 0 && p_x(i) > 0 && p_y(j) > 0
                    I = I + p_xy(i,j) * log2(p_xy(i,j) / (p_x(i) * p_y(j)));
                end
            end
        end
        
        mutual_info(tau) = I;
    end
    
    % 寻找第一个局部极小值
    for tau = 2:max_tau-1
        if mutual_info(tau) < mutual_info(tau-1) && mutual_info(tau) < mutual_info(tau+1)
            tau_opt = tau;
            return;
        end
    end
    
    % 如果没有找到局部极小值，选择互信息下降至初始值1/e的τ
    threshold = mutual_info(1)/exp(1);
    for tau = 2:max_tau
        if mutual_info(tau) < threshold
            tau_opt = tau;
            return;
        end
    end
    
    tau_opt = round(max_tau/2); % 默认值
end

构建高维相空间矩阵

获得最优延迟时间τ后，我们可以构建m维相空间矩阵：

function phase_space = build_phase_space(x, tau, m)
    % 构建高维相空间矩阵
    % 输入:
    %   x - 时间序列
    %   tau - 延迟时间
    %   m - 嵌入维数
    % 输出:
    %   phase_space - 相空间矩阵 (N x m)
    
    N = length(x);
    n_points = N - (m-1)*tau;
    
    phase_space = zeros(n_points, m);
    
    for i = 1:n_points
        for j = 1:m
            idx = i + (j-1)*tau;
            phase_space(i, j) = x(idx);
        end
    end
end

参考代码 互信息法计算延迟时间的一种方法，可用于构建高维相空间矩阵 youwenfan.com/contentcsa/66168.html

应用 Lorenz系统分析

% 生成Lorenz系统时间序列
sigma = 10; rho = 28; beta = 8/3;
dt = 0.01; T = 100; 
t = 0:dt:T;

x = zeros(1, length(t));
y = zeros(1, length(t));
z = zeros(1, length(t));

x(1) = 1; y(1) = 1; z(1) = 1;

for i = 2:length(t)
    dx = sigma*(y(i-1) - x(i-1));
    dy = x(i-1)*(rho - z(i-1)) - y(i-1);
    dz = x(i-1)*y(i-1) - beta*z(i-1);
    
    x(i) = x(i-1) + dx*dt;
    y(i) = y(i-1) + dy*dt;
    z(i) = z(i-1) + dz*dt;
end

% 计算x分量的最优延迟时间
max_tau = 50;
n_bins = 30;
[tau_opt, mutual_info] = mutual_information_delay(x, max_tau, n_bins);

% 绘制互信息函数
figure;
plot(1:max_tau, mutual_info, 'b-o', 'LineWidth', 1.5);
hold on;
plot(tau_opt, mutual_info(tau_opt), 'ro', 'MarkerSize', 10, 'LineWidth', 2);
xlabel('延迟时间 \tau');
ylabel('互信息 I(\tau)');
title('互信息函数');
grid on;
legend('互信息值', '最优延迟时间');
set(gca, 'FontSize', 12);

% 构建相空间
m = 3; % 嵌入维数
phase_space = build_phase_space(x, tau_opt, m);

% 绘制重构的相空间
figure;
plot3(phase_space(:,1), phase_space(:,2), phase_space(:,3), 'b', 'LineWidth', 0.5);
xlabel('x(t)');
ylabel(['x(t+', num2str(tau_opt), ')']);
zlabel(['x(t+', num2str(2*tau_opt), ')']);
title(['Lorenz系统相空间重构 (τ=', num2str(tau_opt), ', m=3)']);
grid on;
rotate3d on;
set(gca, 'FontSize', 12);

% 计算并显示延迟时间
disp(['最优延迟时间 τ = ', num2str(tau_opt)]);

算法优化技巧

1. 自适应分箱方法：

   % 使用Freedman-Diaconis规则确定最优分箱数
   IQR = iqr(x);
   bin_width = 2*IQR*length(x)^(-1/3);
   n_bins = ceil((max(x)-min(x))/bin_width);
   

2. 核密度估计：

   % 使用核密度估计替代直方图
   [density, points] = ksdensity(x);
   

3. 快速互信息计算：

   % 使用快速互信息计算函数
   function I = fast_mutual_info(x, y, bins)
       [joint, ~, ~] = histcounts2(x, y, bins, bins);
       joint = joint/sum(joint(:));
       
       px = sum(joint, 2);
       py = sum(joint, 1);
       
       I = 0;
       for i = 1:bins
           for j = 1:bins
               if joint(i,j) > 0 && px(i) > 0 && py(j) > 0
                   I = I + joint(i,j)*log2(joint(i,j)/(px(i)*py(j)));
               end
           end
       end
   end
   

应用场景

1. 混沌系统分析：Lorenz、Rössler系统等
2. 生物信号处理：EEG、ECG时间序列分析
3. 金融时间序列：股票价格预测
4. 气候数据分析：温度、降水序列分析
5. 机械故障诊断：振动信号特征提取

互信息法的主要优势在于它同时考虑了线性和非线性相关性，能够捕捉时间序列中复杂的依赖关系。通过合理选择延迟时间τ和嵌入维数m，我们可以准确重建原始系统的动力学特性，为后续的非线性时间序列分析（如关联维数、Lyapunov指数计算等）奠定基础。