DFA指数计算的问题

大海绵1990

有没有人计算过DFA指数，下面是我总结的计算DFA指数的步骤，还有我自己编的程序，想拿来计算极端降水阈值，但是运行出来结果有问题，不知道哪里出问题，望各位大神指教！

1.      对某站点湿日（日降水大于等于0.1mm）降水序列xi，计算该降水序列最大值xmax及序列平均值xave。Xave记为XR

2.      从降水序列{xi}的最大值xmax开始，依次舍去{xi,xi>=xmax-d*k}数据区间的数据点，直到xi=R，其实就是说，第一次舍去这个序列的最大的数，第二次舍去倒数第二个，直到xmax-d*k=xave，每舍去一个就得到一个新序列Yj，j=xmax-d*k，d是区间间隔，一般取1，k=1,2,…,(xmax-XR)/d。

3.      计算每个新生序列Yj的长程相关性指数DFAj，分析长程相关性指数随区间j变化的变化过程。对某一长度为m的降水序列Yj={xi，j=1,2,…,m}，长程相关性指数计算过程如下：

（1）      计算降水时间序列{xj,j=1,2,…,m}的累积离差，公式就不写了

（2）      讲累积离差Y(j)等分为N个时间长度为n的区间，这里离差序列不一定能被整除，我为了方便就只能把剩下的数据去掉，这样不太好的，改进的方法可以看珠江的那篇文献。对于分出来的每个区间V(V=1,2,…,N)，利用多项式回归拟合，得到降水序列局部趋势，滤去该趋势后的序列记为Ys(j)，表示原序列与拟和序列之差，公式略。。。

（3）      计算每个子区间滤去局部趋势后的方差，下式中的根号内的部分。

（4）      计算标准DFA的波动函数，y(k)就是上面的累计离差序列Y(j)，yn(k)是拟合出来的序列，得到F(n)与n的关系

（5） 对F(n)与n取对数，得到散点图，并进行线性拟合，其斜率就是DFA指数。

for n=1

    for m=4

        data=mon_pre{n,m};

        xmax=max(data);

        [x]=find(data<0.1);

        data(x,:)=[];

        mean_s=mean(data);

        me=floor(mean_s);

        k=me:1:xmax;

        for o=1:length(k)

            data1=data;

            [z]=find(data1>=k(o));%得到不同长度的新序列

            data1(z,:)=[];

            b=length(data1);

            cc(o)=b;

            mean_l=mean(data1);

            %s1=0;

            for i=1:b

           licha(i)=sum(data1(1:i)-mean_l);

            %for i=1:b

                %s1=s1+(data1(b,1)-mean_l);

                %licha(i)=s1;

              

                %计算累积离差

            end

            

            h=length(licha);

              

            for d=8:30

               

       w=floor(h/d);    %几种分段数w

         licha1=licha;

       N1=w*d;   %分段后的总的数据长度

         if h>N1

         licha1(:,N1+1:h)=[];

         end

       A=reshape(licha1,w,[]);

           ys=[];

            for q=1:w

                e=1:1:d;

            pv=polyfit(e,A(q,:),1);%得到所有段的线性拟合

            

            yp=polyval(pv,e);

            ys=[ys,A(q,:)-yp];  %得到一种分段的滤趋势

            end

            

            

            j=0;

            

            for g=1:N1

             j=j+(licha1(g)-ys(g)).^2;

             end

             ff=sqrt(j/N1);

             F_w(d)=ff;

            

           

            end

           d=8:1:30;

           p=polyfit(log10(d),log10(F_w(:,8:30)),1);

           pk(o)=p(1,1);

        end

    end

end

刚接触matlab不久，程序写的比较繁琐，大神们有好的建议可以告诉我哈！

山_脉

我也不懂，不过我觉得你的结论和你已经不在遥远了

山_脉

哦，对了，楼主是干嘛的呢，怎么会习惯称降水日为湿日呢。

山_脉

你去看看论坛的这个帖子，看看对你有帮助不http://bbs.06climate.com/forum.php?mod=viewthread&tid=19701，祝你早日成功。

大海绵1990

山_脉 发表于 2014-9-15 15:32
你去看看论坛的这个帖子，看看对你有帮助不http://bbs.06climate.com/forum.php?mod=viewthread&tid=19701 ...

谢谢，其实我看到了，但是那个我不知道怎么用到自己的数据上去

大海绵1990

山_脉 发表于 2014-9-15 15:07
我也不懂，不过我觉得你的结论和你已经不在遥远了

谢谢只是觉得这一步跨得好久

大海绵1990

山_脉 发表于 2014-9-15 15:08
哦，对了，楼主是干嘛的呢，怎么会习惯称降水日为湿日呢。

文献里面的，就是这样定义的，还有好的方法可以推荐给我哦

山_脉

大海绵，最近你的这个问题弄清楚了吗，还是去研究其他东西去了

大海绵1990

山_脉 发表于 2015-4-9 10:28
大海绵，最近你的这个问题弄清楚了吗，还是去研究其他东西去了

之前是基本弄出来了，但是好长时间没弄了，现在在研究其他的东西

scwjx520

谢谢楼主,正要用的.