hillplot.m

极值理论中各种函数及图像的程序。matlab实现。

function [c,uband,lband]=hillplot(data,param,opt,ci,plotlen)
%Plots the Hill estimate of the tail index of heavy-tailed data
%
%  USAGE: [out,uband,lband]=hillplot(data,param,opt,ci,plotlen)
%
%   data: Data vector
%  param: whether 'alpha' or 'xi' (1/alpha) should be plotted
%    opt: Option for plotting with respect to threshold ('t') or number of exceedances 
%         ('n').
%     ci: Confidence interval
%plotlen: Optional argument that determines plotting range of hill estimate starting from
%         first order statistics. It can be entered as an integer (no quotes) or as a 
%         percentage as '%5' including the quotes.
%
%    out: Parameter plotted
%  uband: Upper confidence interval
%  lband: Lower confidence interval
	

data=surecol(data);
ordered=flipud(sort(data));
ordered=ordered(ordered>0);
n=length(ordered);
loggs=log(ordered);
avesumlog=cumsum(loggs)./(1:n)';
diffi=avesumlog-loggs;
diffs=[NaN; diffi(2:n)];
if strcmp(param,'alpha'),
    diffs=1./diffs;
elseif strcmp(param,'xi'),
    diffs=diffs;
else disp('Parameter should be one of ''xi'' or ''alpha'' ');
    return
end
 
ses=diffs./sqrt(1:n)';
x=1:n;
y=diffs;
last=n;
if nargin==5,
   if isstr(plotlen)==1,
		perc=str2num(plotlen(2:end))/100;
		last=floor(perc*n);
   else
      last=plotlen;
   end
end
if strcmp(opt,'n') 
plot(x(1:last),y(1:last))
c=y;
qq=norminv(1-(1-ci)/2);
uband=y+ses(x)*qq;
lband=y-ses(x)*qq;
hold on
plot(x(1:last),uband(1:last),'r:');
plot(x(1:last),lband(1:last),'r:');
xlabel('Order Statistics')
ul=[uband;lband];
scx=range(x(1:last))/20;
scy=range(ul(1:last,:))/20;
axis([min(x(1:last))-scx max(x(1:last))+scx min(ul)-scy max(ul(1:last,:))+scy]);


elseif strcmp(opt,'t'),
   
   threshs=findthresh(data,1:n);
   plot(threshs(1:last),y(1:last));
	c=y;
	qq=norminv(1-(1-ci)/2);
	uband=y+ses(x)*qq;
	lband=y-ses(x)*qq;
	hold on
	plot(threshs(1:last),uband(1:last),'r:');
	plot(threshs(1:last),lband(1:last),'r:');

	xlabel('Thresholds')
   ul=[uband;lband];
   scx=range(threshs(1:last))/20;
   scy=range(ul(1:last,:))/20;
	axis([min(threshs(1:last))-scx max(threshs(1:last))+scx min(ul)-scy max(ul(1:last,:))+scy]);
   set(gca,'XDir','Reverse');
else 
   error('Last argument should be ''n'' or ''t''')
end

if strcmp(param,'alpha'),
    ylabel('alpha');
end
if strcmp(param,'xi'),
    ylabel('xi');
end
title('Hillplot')
hold off