gpd_q.m

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

function [out,xp,parmax]=gpd_q(resgpd,firstdata,cprob,ci_type,ci_p,like_num),
%Calculates quantile estimates and confidence intervals for high quantiles 
%above the threshold in a GPD analysis
%
%    USAGE: [out,xp,parmax]=gpd_q(resgpd,firstdata,cprob,ci_type,ci_p,like_num),
%
%   resgpd: gpd fit structure
%firstdata: Data obtained from gpd_plot as: firstdata=gpd_plot(...), choose option 2 then 
%           choose 0 and exit
%    cprob: Cumulative tail probability.
%  ci_type: Confidence interval calculation method, 'wald' or 'likelihood'
%     ci_p: Confidence interval probability
% like_num: number of likelihoods to be calculated
%
%      out: Structure for confidence intervals calculated
%       xp: index for likelihoods
%   parmax: likelihoods


threshold=resgpd.threshold;
par_ests=resgpd.par_ests;
xihat=par_ests(1);
betahat=par_ests(2);
varcov=resgpd.varcov;
p_less_thresh=resgpd.p_less_thresh;    
lambda=1/(1-p_less_thresh);

a=lambda*(1-cprob);
q = threshold + betahat * gfunc(a, xihat);
if strcmp(ci_type,'wald'),
   scaling = threshold;
	betahat = betahat/scaling;
	xivar = varcov(1, 1);
	betavar = varcov(2, 2)/(scaling^2);
	covar = varcov(1, 2)/scaling;
	term1 = betavar * (gfunc(a, xihat))^2;
	term2 = xivar * (betahat^2) * (gfunc_deriv(a, xihat))^2;
	term3 = 2 * covar * betavar * gfunc(a, xihat) * gfunc_deriv(a, xihat);
	qvar = term1 + term2 + term3;
   if(qvar < 0)
		disp('Negative estimate of quantile variance');
      return
   end
    
   qse = scaling * sqrt(qvar);
	qq = norminv(1 - (1 - ci_p)/2);
	upper = q + qse * qq;
	lower = q - qse * qq;
	out.Upper=upper;
   out.Lower=lower;
   out.Estimate=q;
   out.Std_Err=qse;
    
   loglog(firstdata(:,1),firstdata(:,2),'.');
   set(gca,'PlotBoxAspectRatio',[1 0.7 1]);
   set(gca,'yticklabel',get(gca,'ytick'));
   set(gca,'ytickmode','manual');
   xlabel('x on log scale')
   ylabel('1-F(x) on log scale')
    
   hold on
   
   x=qgpd(linspace(0,1,1000),xihat,threshold,betahat*scaling);
   y=pgpd(x,xihat,threshold,betahat*scaling);
   prob = resgpd.p_less_thresh;
   y = (1 - prob) * (1 - y);
   plot(x(y>0),y(y>0),'k');
   set(gca,'FontSize',8)
   yy=get(gca,'ylim');yl=yy(end);
   
   y2=pgpd(out.Estimate,xihat,threshold,betahat*scaling);
   y2 = (1 - prob) * (1 - y2);

   h=plot([out.Lower,out.Upper],[y2,y2],'r');
   set(h,'linewidth',4);
   plot(out.Estimate,y2,'y.');

   plot([out.Lower,out.Lower],[yy(1),yy(end)],':');
   h=text(out.Lower,yl/0.95,num2str(out.Lower));
   set(h,'fontsize',8);set(h,'rotation',90);
        
   plot([out.Upper,out.Upper],[yy(1),yy(end)],':');
   h=text(out.Upper,yl/0.95,num2str(out.Upper));
   set(h,'fontsize',8);set(h,'rotation',90);
        
   plot([out.Estimate,out.Estimate],[yy(1),yy(end)],'r:');
   h=text(out.Estimate,yl/0.95,num2str(out.Estimate));
   set(h,'fontsize',8);set(h,'rotation',90);
    
   hold off
   
elseif strcmp(ci_type,'likelihood'),
   
   tmp=resgpd.data-threshold;
   theta = xihat;
	parmax = [];
   xp=exp(linspace(log(threshold),log(max(resgpd.data)*1.5),like_num));
   u=threshold;    
   opts=optimset('MaxFunEvals',1000,'MaxIter',200,'TolX',1e-8,'TolFun',1e-8,'Display','off');
 
   for i=1:length(xp),
      xpi=xp(i);
      [e,f,g,h]=fminsearch('parloglik',theta,opts,tmp,a,u,xpi);
      %[e,f,g,h]=fminunc('parloglik',theta,opts,tmp,a,u,xpi);                       
      parmax = [parmax, parloglik(e,tmp,a,u,xpi)];
      
   end
        
   parmax =  - parmax;
   xpi=q;
	overallmax =  - parloglik(xihat,tmp,a,u,xpi);
	crit = overallmax - chi2inv(0.999, 1)/2;
	cond = (parmax > crit);
	xp = xp(cond);
   parmax = parmax(cond);
	    
   loglog(firstdata(:,1),firstdata(:,2),'.');
   set(gca,'PlotBoxAspectRatio',[1 0.7 1]);
   x=qgpd(linspace(0,1,1000),xihat,threshold,betahat);
   y=pgpd(x,xihat,threshold,betahat);
   prob = resgpd.p_less_thresh;
   y = (1 - prob) * (1 - y);
   hold on
   plot(x(y>0),y(y>0),'k');
   set(gca,'FontSize',8)


	%parmax2=parmax/((max(parmax)-min(parmax))/max(firstdata(:,2))-min(firstdata(:,2)));
   %parmax2=parmax2-max(parmax2)+max(firstdata(:,2));
        
   xpp=linspace(xp(1),xp(end),4*like_num);
   parmaxx=spline(xp,parmax,xpp);
   xlabel('x on log scale')
   ylabel('1-F(x) on log scale')

        %subplot(212)
        %semilogx(xpp,parmaxx,'r');
        %set(gca,'XGrid','on')

   con = overallmax-chi2inv(ci_p, 1)/2;
   [y,j]=max(parmax);
        
   
   %subplot(212)
   %hold on
   %xs=get(gca,'Xtick');
   %plot(xs,repmat(con,1,length(xs)),':')
   %con2 = overallmax-chi2inv(0.99, 1)/2;
	%plot(xs,repmat(con2,1,length(xs)),'g:')

   %set(gca,'FontSize',8)
   %xlabel(['                                       x on log scale                                  ';'(Dotted blue is for the specific confidence entered, 0.99 confidence band dotted green)']);
   %hold off
   
   a=xpp(find(parmaxx>=con));
   out.Lower=min(a);
   out.Estimate=q;
   out.Upper=max(a);
   y2=pgpd(out.Estimate,xihat,threshold,betahat);
   y2 = (1 - prob) * (1 - y2);
   
   %hold on
   %subplot(211)    
   
   h=plot([out.Lower,out.Upper],[y2,y2],'r');
   set(h,'linewidth',4);
   plot(out.Estimate,y2,'y.');
   set(gca,'yticklabel',get(gca,'ytick'));
   set(gca,'ytickmode','manual');
        
   
   yy=get(gca,'ylim');yl=yy(end);
        
   plot([out.Lower,out.Lower],[yy(1),yy(end)],':');
   h=text(out.Lower,yl/0.95,num2str(out.Lower));
   set(h,'fontsize',8);set(h,'rotation',90);
   
   plot([out.Upper,out.Upper],[yy(1),yy(end)],':');
   h=text(out.Upper,yl/0.95,num2str(out.Upper));
   set(h,'fontsize',8);set(h,'rotation',90);
        
   plot([out.Estimate,out.Estimate],[yy(1),yy(end)],'r:');
   h=text(out.Estimate,yl/0.95,num2str(out.Estimate));
   set(h,'fontsize',8);set(h,'rotation',90);
   
   hold off
   
else 
   disp('The 4th input argument should be ''likelihood'' or ''wald'' ')
   return
   
end
     
    
function c=gfunc(a,xihat),
c=(a^(-xihat)-1)/xihat;
    
function c=gfunc_deriv(a,xihat),
c=(-(a^(-xihat)-1)/xihat-a^(-xihat)*log(a))/xihat;