gp_algorithm.m

GP算法计算时间序列想空间重构的最佳嵌入维数和延迟时间

% GP算法求关联维和嵌入维

%clc
clear
close all

IsPlot = 1;
%---------------------------------------------------
% 产生 Lorenz 时间序列

sigma = 10;          % Lorenz方程参数
r = 28;               
b = 8/3;          

y = [-1;0;1];        % 起始点 (3x1 的列向量)
h = 0.01;            % 积分时间步长

k1 = 30000;           % 前面的迭代点数
k2 = 3000;           % 后面的迭代点数

z = LorenzData(y,h,k1+k2,sigma,r,b);
X = z(k1+1:end,1);
X = normalize_1(X);  % 归一化到均值为 0，振幅为 1

%---------------------------------------------------'

disp('----- GP算法求关联维和嵌入维 -----');

t = 10;
m_vector = 2:5;
r_vector = exp(-5:0.5:0);

num_m = length(m_vector);
num_r = length(r_vector);
ln_Cr = zeros(num_m,num_r);

%------------------------------------------------------
tic
type_norm = 2       % 使用范数类型 (缺省值为2)
                    % type_norm = 0,1,2时，分别对应无穷范数、1范数和2范数
block = 1           % 分块计计算关联积分 - 分块数 (缺省值为1)
                    % t越大速度越快，但有误差
for i = 1:num_m
    for j = 1:num_r
        % 计算关联积分S(m,N,r,t), 参见 <<混沌时间序列分析及应用>> P35 式(2.29)
        m = m_vector(i);
        r = r_vector(j);
        %ln_Cr(i,j) = log(CorrelationIntegral(m,X,r,t)); % 缺省用法
        ln_Cr(i,j) = log(CorrelationIntegral(m,X,r,t,type_norm,block)); 
    end
end
t = toc

ln_r = log(r_vector);
plot(ln_r,ln_Cr','+:');grid;
xlabel('ln(r)'); ylabel('ln(C(r))');
title(['norm = ',num2str(type_norm),', block = ',num2str(block),', t = ',num2str(t)]);
legend('m=2','m=3','m=4','m=5',4)