G_P算法中参数r的确定

xiaocheng_2007

在G_P算法中r是如何确定的？我看了版主发的总结帖，还是心存疑惑，请问大家有没有这方面的理论知识？

xiaocheng_2007

大家是怎么选择的呀？

yuling

首先计算出重构的相空间中相点之间的最大距离dmax和最小距离dmin,选取一定的步长令r从
dmin到dmax逐步变化，这么做比合理。

xiaocheng_2007

这两天我做了一个实验，发现有这样一种现象：如果r的长度少于20，会增大lnC出现NaN的机率。大家有没有这样的情况？

sandman

http://www.studa.net/yingyong/080505/08422984.html
这个文章不知道看过没？‘一种简单的GP算法无标度区识别方法’

板凳 yuling  你的g-p算法的程序能贴一下吗？我的那个运算量特别大5000个数据我要运行40多分钟，我用的是matlab6.5.1，不知道你运行时多少时间？ 学习下你的程序。。。。:@)

xiaocheng_2007

楼上，你运行时间挺快的呀，我要运行在一个小时以上:@L

yuling

好的，这是我的程序，我在编程的时候尽量减少子函数的调用，这可以显著提高运行的速度。
function [ln_r,ln_C]=G_P(data,tau,min_m,max_m,ss)
% 本函数是利用G-P 方法计算混沌吸引子关联维
% data::待计算的时间序列
% tau:  时间延迟
% min_m:最小嵌入维
% max_m:最大嵌入维
% ss:半径搜索次数
N=length(data);   %待计算的时间序列长度
ln_C=zeros(max_m-min_m+1,ss);
ln_r=zeros(max_m-min_m+1,ss);
for m=min_m:max_m
    Y=reconstitution(data,m,tau);%重构相空间
    M=N-(m-1)*tau;%相空间点的数目
    d=zeros(M-1,M);
    for i=1:M-1
        for j=i+1:M
            d(i,j)=max(abs(Y(:,i)-Y(:,j)));%计算相空间中每两点之间的距离           
        end                                
    end
    max_d=max(max(d));%相空间中两点之间的最大距离
     for i=1:M-1      %计算相空间中两点之间的最小距离
        for j=1:i
            d(i,j)=max_d;   
        end                              
     end
    min_d=min(min(d));%相空间中两点之间的最小距离
    delt=(max_d-min_d)/ss;%搜索半径增加的步长
    for k=1:ss
        r=min_d+k*delt;
        C(k)=correlation_integral(Y,M,r);  %计算关联积分
        ln_C(m-min_m+1,k)=log(C(k));  %求lnC(r)
        ln_r(m-min_m+1,k)=log(r);     %求lnr
    end
    plot(ln_r(m-min_m+1,:),ln_C(m-min_m+1,:));%画出双对数图
    hold on;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function Data=reconstitution(data,m,tau)
%该函数用来重构相空间
% m:嵌入空间维数
% tau:时间延迟
% data:输入时间序列
% Data:输出,是m*n维矩阵

N=length(data); % N为时间序列长度
M=N-(m-1)*tau; %相空间中点的个数
Data=zeros(m,M);
for j=1:M
  for i=1:m           %相空间重构
    Data(i,j)=data((i-1)*tau+j);
  end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function C_I=correlation_integral(X,M,r)
%该函数用来计算计算关联积分
%C_I:关联积分的返回值
%X:重构的相空间矢量，是一个m*M的矩阵
%M::M是重构的m维相空间中的总点数
%r:Heaviside 函数中的搜索半径
sum_H=0;
for i=1:M-1
    for j=i+1:M
        d=norm((X(:,i)-X(:,j)),inf);%计算相空间中每两点之间的距离，其中NORM(V,inf) = max(abs(V)).
        if r>d   
        %sita=heaviside(r,d);%计算Heaviside 函数之值n
           sum_H=sum_H+1;
        end
    end
end
C_I=2*sum_H/(M*(M-1));%关联积分的值

heqiangyong

不对、，算出来的都是等于零！

yuling

哪地方不对啊，能具体点吗

heqiangyong

我用数据运行了一次，得出的结果都是inf，不知道错在哪？

yuling

我用如下Lorenz方程跑了一下本程序，结果很正常。你的数据和采用的参数怎么得到的，可以告诉我吗，大家可以一起分析。
function dy = Lorenz(t,y)
a=16;
b=4.0;
c=45.92;
dy=zeros(3,1);
dy(1)=-a*(y(1)-y(2));
dy(2)=-y(1)*y(3)+c*y(1)-y(2);
dy(3)=y(1)*y(2)-b*y(3);

%%%%%%%%%%%%%%%%%%%%
clear
t0=0;
tf=130;
[t,x]=ode45(@Lorenz,[t0:0.01:tf],[-1,0,1]);
Lorenz_data=x(10002:end,1);
tic
[lnr,lnC]=G_P(Lorenz_data,10,9,10,40);
toc

李锦鹏2008268

给我发一个编好的程序吗？我是初学者，谢谢，421069750@163.com