实现方法来自论文:
Quantification of scaling exponents and crossover phenomena in nonstationary
heartbeat time series
C.‐K. Peng
,
Shlomo Havlin
,
H. Eugene Stanley
, and
Ary L. Goldberger
输入:一个数组B,长度为len
输出:斜率b(2)
第一步:计算y(k),其中Bave位B的均值。
第二步:计算yn(k),n是能被len整除的数,比如我这是3500个数,那么n就可以取1,2,5,10,35等。
根据n的大小划分,比如我取n=500,那么y就被划分成了7块,对应了1-500;501-1000;1001-1500;1501-2000;2001-2500;2501-3000;3001-3500。
然后计算每一块的regress(即利用回归分析得到预测值)。
第三步:根据如下公式计算F(n)
第四步:得到n和F(n),以log10为底,然后画图,再次用regress得到斜率b。
function out=dfa(a)
avg=mean(a);
sum=0.0;
len=length(a);
for i=1:3500
sum=sum+a(i)-avg;
y(i)=sum;
end
%x=1:3500;
%X=[ ones(length(y),1), x' ];
%b=regress(y',X)
index=1;
for n=10:3500
if(rem(3500,n)==0)
count=0;
k=3500/n;
for i=1:k
begin=(i-1)*n+1;
finish=i*n;
q=y(begin:finish);
p=begin:finish;
X=[ ones(length(q),1), p' ];
b=regress(q',X); %回归 y=b1+b2*x;
for j=begin:finish
cc=b(1)+b(2)*j;
op=y(j)-cc;
count=count+op*op;
end
end
f=(count/3500)^0.5;
K(index)=log10(n);
F(index)=log10(f);
index=index+1;
end
end
X=[ ones(length(F),1), K' ];
[b,bint,r,rint,stats]=regress(F',X);
out=b(2);
z=b(1)+b(2)*K;
figure,plot(K,F,'rp',K,z,'b');
title('散点图及回归线');
end