R/S分析源代码

sciren

是不是你给的连接里的那个代码是对的啊

sciren

诶 那个跑不起来

sciren

sciren 发表于 2018-11-5 21:34
酱紫 好的  那你还改吗 你改的话也上传吧  我就可以两个都保存啦 嘿嘿 （不要打我）

她的代码算不出来，，说未定义V

Jude

sciren 发表于 2018-11-5 21:55
她的代码算不出来，，说未定义V

喔，那个代码是有点问题的，只是针对那个人写的文章的数据，我改一下，更通用一些。

Jude

sciren 发表于 2018-11-5 21:55
她的代码算不出来，，说未定义V

function hurst=hurst_test(data)
[M,N]=size(data);
for n=3:N%子区间长度
    H=fix(N/n);
    Rs=zeros(1,H);
    for h=1:H%子区间个数
        xna=zeros(1,n);
        x=data(1,n*h-n+1:n*h);
        for r=1:n
            for t=1:r
                xna(r)=xna(r)+(x(t)-mean(x));%计算累积极差
            end
        end
        R=max(xna)-min(xna);
        s=std(x,1);
        Rs(h)=R/s;
    end
    Rsn(n-2)=mean(Rs);
end
n0=3:N;
lnRs=log(Rsn(:));
lnn=log(n0);
p2=polyfit(lnn,lnRs',1);
hurst=p2(1); % 赫斯特指数
v=Rsn./n0.^0.5;
figure(1)
plot(lnn,v);
figure(2)
scatter(lnn,lnRs,'o','k');hold on;
plot(lnn,lnn.*hurst+p2(2),'r--');
xlabel('ln(t)');ylabel('ln(R/S)');
r2=corrcoef(lnn,lnRs'); %拟合R的平方
if (p2(2)<0)
    str1=['y=',num2str(p2(1)),'*x',num2str(p2(2))];
else
    str1=['y=',num2str(p2(1)),'*x+',num2str(p2(2))];
end
a=get(gca);
text(a.XLim(1)+0.2,a.YLim(2)-0.2,str1);  %依据横纵坐标上下限，标注text
str2=['H=',num2str(p2(1))];str3=['R^2=',num2str(r2(1,2)^2)];
text(a.XLim(1)+0.2,a.YLim(2)-0.4,str2); text(a.XLim(1)+0.2,a.YLim(2)-0.6,str3);
end

Jude

sciren 发表于 2018-11-5 21:55
她的代码算不出来，，说未定义V

按照那篇博文，重新改过的代码，输入数据为行向量，您用您的数据验证一下是否合理？

sciren

Jude 发表于 2018-11-6 11:08
按照那篇博文，重新改过的代码，输入数据为行向量，您用您的数据验证一下是否合理？

哇咔咔  我去试试  上午去开会了  嘿嘿

sciren

不知道为啥  我也不知道我看的那个连接具体是怎么算的  你的和他的结果大小相差0.2左右 但是具体都是在量级范围内

sciren

%R/S分析Matlab程序  
function [logRS,logERS,V,p]=RSana(x,n,method,q)
% Alexandros Leontitsis
% Department of Education
% University of Ioannina
% 45110 - Dourouti
% Ioannina
% Greece
%
% University e-mail: me00743@cc.uoi.gr
% Lifetime e-mail: leoaleq@yahoo.com
% Homepage: http://www.geocities.com/CapeCanaveral/Lab/1421
%
% 1 Jan 2004.

if nargin<1 || isempty(x)==1
   error('You should provide a time series.');
else
   % x must be a vector
   if min(size(x))>1
      error('Invalid time series.');
   end
   x=x(:);
   % N is the time series length
   N=length(x);
end

if nargin<2 || isempty(n)==1
   n=1;
else
   % n must be either a scalar or a vector
   if min(size(n))>1
      error('n must be either a scalar or a vector.');
   end
   % n must be integer
   if n-round(n)~=0
       error('n must be integer.');
   end
   % n must be positive
   if n<=0
     error('n must be positive.');
   end
end

if nargin<4 || isempty(q)==1
   q=0;
else
    if q=='auto'
        t=autocorr(x,1);
        t=t(2);
        q=((3*N/2)^(1/3))*(2*t/(1-t^2))^(2/3);
    else
        % q must be a scalar
        if sum(size(q))>2
            error('q must be scalar.');
        end
        % q must be integer
        if q-round(q)~=0
            error('q must be integer.');
        end
        % q must be positive
        if q<0
            error('q must be positive.');
        end
    end
end

logn=log10(n);
for i=1:length(n)
   
    % Calculate the sub-periods
    a=floor(N/n(i));
   
    % Make the sub-periods matrix
    X=reshape(x(1:a*n(i)),n(i),a);
   
    % Estimate the mean of each sub-period
    ave=mean(X);
   
    % Remove the mean from each sub-period
    cumdev=X-ones(n(i),1)*ave;
   
    % Estimate the cumulative deviation from the mean
    cumdev=cumsum(cumdev);
   
    % Estimate the standard deviation
    switch method
    case 'Hurst'
        % Hurst-Mandelbrot variation
        stdev=std(X);
    case 'Lo'
        % Lo variation
        for j=1:a
            sq=0;
            for k=0:q
                v(k+1)=sum(X(k+1:n(i),j)'*X(1:n(i)-k,j))/(n(i)-1);
                if k>0
                    sq=sq+(1-k/(q+1))*v(k+1);
                end
            end
            stdev(j)=sqrt(v(1)+2*sq);
        end
    case 'MW'
        % Moody-Wu variation
        for j=1:a
            sq1=0;
            sq2=0;
            for k=0:q
                v(k+1)=sum(X(k+1:n(i),j)'*X(1:n(i)-k,j))/(n(i)-1);
                if k>0
                    sq1=sq1+(1-k/(q+1))*(n(i)-k)/n(i)/n(i);
                    sq2=sq2+(1-k/(q+1))*v(k+1);
                end
            end
            stdev(j)=sqrt((1+2*sq1)*v(1)+2*sq2);
        end
    case 'Parzen'
        % Parzen variation
        if mod(q,2)~=0
            error('For the "Parzen" variation q must be dived by 2.');
        end
        for j=1:a
            sq1=0;
            sq2=0;
            for k=0:q
                v(k+1)=sum(X(k+1:n(i),j)'*X(1:n(i)-k,j))/(n(i)-1);
                if k>0 && k<=q/2
                    sq1=sq1+(1-6*(k/q)^2+6*(k/q)^3)*v(k+1);
                elseif k>0 && k>q/2
                    sq2=sq2+(1-(k/q)^3)*v(k+1);
                end
            end
            stdev(j)=sqrt(v(1)+2*sq1+2*sq2);
        end
    otherwise
        error('You should provide another value for "method".');
    end
   
    % Estiamte the rescaled range
    rs=(max(cumdev)-min(cumdev))./stdev;
   
    clear stdev
   
    % Take the logarithm of the mean R/S
    logRS(i,1)=log10(mean(rs));
   
    if nargout>1
        
        % Initial calculations fro the log(E(R/S))
        j=1:n(i)-1;
        s=sqrt((n(i)-j)./j);
        s=sum(s);
        
        % The estimation of log(E(R/S))
        logERS(i,1)=log10(s/sqrt(n(i)*pi/2));
        
        % Other estimations of log(E(R/S))
        %logERS(i,1)=log10((n(i)-0.5)/n(i)*s/sqrt(n(i)*pi/2));
        %logERS(i,1)=log10(sqrt(n(i)*pi/2));
        
    end
   
    if nargout>2
        % Estimate V
       V(i,1)=mean(rs)/sqrt(n(i));
    end
end
    p=polyfit(logn,logRS,1);
  

%
%需要将该程序命名为RSana.m保存，以供调用。
%在matlab中导入所要分析的时间序列。
%然后在matlab命令行输入命令logRS=RSana(x,n,'Hurst')。其中x就是该时间序列，它本身就是向量数据;n为分析的时间周期，Hurst是分析方法。返回值logRS。
%这里指对一个周期进行分析，返回返回值logRS是一个标量。如果要作log-log图，需要一个循环，分别对n=5:200进行计算logRS。
我看有个人的代码这样的 是不是和分周期有关 所以结果不一样啊

sciren

诶不对 为啥我输入data=xlsread('D:\Desktop\rs\shuju.xlsx');
hurst=hurst_test(data)
会提示我未定义与 'double' 类型的输入参数相对应的函数 'hurst_test'。