0-1 test
0-1测试是一个能够衡量时间序列是否有混沌的一种测试算法,与李亚普洛夫指数不同的是,它不需要进行相空间重构,通过输出结果是否接近于1来判别混沌现象的产生。
算法步骤
1.对于一个常数
c
∈
(
0
,
π
)
c\in(0,\pi)
c∈(0,π),计算:
这里
n
=
1
,
2
,
.
.
.
,
N
n=1,2,...,N
n=1,2,...,N。我们可以画出分别以
p
c
p_c
pc和
q
c
q_c
qc为横轴和纵轴的轨迹图,其图形产生的布朗运动可以表征混沌现象。
2.为了分析
p
c
p_c
pc和
q
c
q_c
qc的扩散行为,我们可以计算位移均方差MSD,计算式为:
这里,一般
n
=
N
/
10
n=N/10
n=N/10,但是这里的
N
N
N如果采用原始序列的长度,那么横线处的
p
c
p_c
pc将超出本身大小,所以我认为
N
=
N
−
n
N=N-n
N=N−n,这块我查阅了一些参考文献,都没有给出具体说明,我认为还可以讨论。
3.本身算出
M
c
(
n
)
M_c(n)
Mc(n)后可以直接计算
k
c
k_c
kc,但是一般我们还需对于结果进行修正,首先我们给出原始算法被称为回归算法。
4.这里我们再给出修正后的算法
这里需要根据(7.4)计算得出
D
c
(
n
)
D_c(n)
Dc(n),然后进行计算即可;
这里
a
a
a一般取1.1。接着计算
k
c
k_c
kc
可以得到
k
c
k_c
kc的修正值。
5.文献中还给出了一种相关方法计算
k
c
k_c
kc
这种方法的好处是计算结果能够比较好的接近于1或者0,符合0-1测试的内涵。
下面分享下自己编写的MATLAB代码,希望交流指正。
function [Kc,Kcreg,KcCorr] =ZeroOneTest(dataSet,cont)
%函数功能:完成0-1测试
%输入参数说明
%dataSet:待测试序列
%cont:规定常数值,范围0-pi自选
%------------------------
%输出参数说明
%Kc:利用Mc计算
%Kcreg:利用修正D得到
%------------------------
%测试样例
%dataset = logisticmap(0.3,3.97,4999);
%cont = 2;
% [Kc,Kcreg,KcCorr] =
ZeroOneTest(dataSet,cont);
%参考文献:The 0-1 Test
for Chaos: A Review
%计算开始
Ndata = size(dataSet,1);
p(1)=dataSet(1)*cos(cont);
s(1)=dataSet(1)*sin(cont);
for n = 1:Ndata
p(n+1)=p(n)+dataSet(n)*cos(n*cont);
s(n+1)=s(n)+dataSet(n)*sin(n*cont);
end
% %画出p和s的轨迹图
% figure
% plot(p,s);
% xlabel('p_c');ylabel('q_c');
% %画图结束;
%计算均方位移
Numn = Ndata/10;%建议不超过数据集大小的十分十一
meanpower = mean(dataSet,1)^2;
for n = 1:Numn
sumMean = 0;
for j = 1:Ndata-Numn%此处公式中N的大小取剩余的值
sumMean = sumMean + (p(j+n)-p(j))^2 + (s(j+n)-s(j))^2;
end
Mc(n) = 1/Ndata*sumMean;
Vosc(n) =
meanpower*((1-cos(n*cont))/(1-cos(cont)));
end
D = Mc -Vosc;
x = 1:Numn;
KcCorr = corr(x',D');
% %画图:画出Mc和D的变化
% figure
% n = 1:Numn;
% plot(n,Mc);
% hold on
% plot(n,D);
% xlabel('n');ylabel('Mc,D');
% %结束画图
Dreg = D(Numn) + 1.1*min(abs(D));
Kc = log(Mc(Numn))/log(Numn);
Kcreg = log(Dreg)/log(Numn);
end
测试代码
function [medianKcreg,medianKcCorr] =calcuKcreg(dataSet)
%函数功能:01测试进行100次,随机选择常数cont,得到Kcreg修正值
%输入参数:dataSet测试数据集
%输出参数:medianKcreg计算出来的Kcreg的中位数
for Numcalcu = 1:100
cont = pi/5+3*pi/5.*rand(1);
[Kc,Kcreg,KcCorr] = ZeroOneTest(dataSet,cont);
Kcreg1(Numcalcu) = Kcreg;
KcCorr1(Numcalcu)
= KcCorr;
end
medianKcreg = median(Kcreg1,2);
medianKcCorr = median(KcCorr1,2);
end
共同学习,欢迎留言指正。