den.m

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

function res=gpd(data,threshold,nextremes,information),
%Fits a generalized Pareto model to excesses over a high threshold
%
%USAGE: res=gpd(data,threshold,nextremes,information)
%
%!Either threshold or nextremes should be defined. The undefined one should be entered as '[]'
%
%data: Data vector
%threshold: Excesses over this value will be fitted a model
%nextremes: Implies a threshold value that number of observations remaining above is nextremes
%information: Default is 'observed'. Can be entered as 'expected' also. Determines whether standard
%             errors will be calculated with observed or expected information
%
%res: Fitted distribution
%   res.par_ests: Estimated parameters. 1X2 vector: 
%                                                   1st element: xi
%                                                   2nd element: beta
%   res.funval: Value of the negative log likelihood
%   res.terminated: Termination condition. 1 if successfully terminated
%   res.details: Details of the nonlinear minimization process of the negative likelihood
%   res.varcov: Variance-covariance matrix of the parameters
%   res.par_ses: Standard deviations of the parameters of the distribution
%   res.data: Excesses over the threshold


warning off
n=length(data);
if (isempty(threshold)&isempty(nextremes))
    disp('Enter either a threshold or the number of upper extremes')
    return
end
if (~isempty(threshold)&~isempty(nextremes))
    disp('Enter Either a threshold or the number of upper extremes')
    return
end

if (~isempty(nextremes))
    threshold=findthresh(data,nextremes);
end
if nargin<4,
    information='observed';
end
    

exceedances=data(data>threshold);
excess=exceedances-threshold;
xbar=mean(excess);
s2=var(excess);
xi0=-0.5*xbar*(((xbar^2)/s2)-1);
beta0=0.5*xbar*(((xbar^2)/s2)+1);

theta=[xi0,beta0];

opts=optimset('MaxFunEvals',1000,'MaxIter',200,'TolX',1e-8,'TolFun',1e-8);
%xi=theta(1);
%beta=theta(2);
%cond1 = beta <= 0;
%cond2 = ((xi <= 0) & (max(excess) > ( - beta/xi)));
%if (cond1 | cond2),
%    theta(1)=1;
%    theta(2)=1;
%end
[res.par_ests,res.funval,res.terminated,res.details] = fmincon('negloglikgpd',theta,[0 -1],0,[],[10,inf],[],[0,0],'nlcongpd',opts,excess);
%[res.par_ests,res.funval,res.terminated,res.details] = fminunc('negloglikgpd',res.par_ests,opts,excess);
if strcmp(information,'observed'),
    res.varcov=hessigpd('negloglikgpd',res.par_ests,excess);
    res.par_ses=sqrt(diag(res.varcov))';
elseif strcmp(information,'expected'),
    one = (1 + res.par_ests(1))^2/length(excess);
	two = (2 * (1 + res.par_ests(1)) * res.par_ests(2)^2)/length(excess);
	cov =  - ((1 + res.par_ests(1)) * res.par_ests(2))/length(excess);
	res.varcov = [one,cov;cov,two];
    res.par_ses=sqrt(diag(res.varcov))';
else
    
    
    disp('WARNING 4th input should be either observed or expected');
    return
end

res.threshold=threshold;
res.data=exceedances;
res.p_less_thresh=1-length(excess)/n;


function f=negloglikgpd(theta,excess),
    xi = theta(1);
    
    beta = theta(2);
	
    cond1 = beta <= 0;
	cond2 = (xi <= 0) & (max(excess) > ( - beta/xi));
    if(cond1 | cond2)
		f=inf;
    else
        y = log(1 + (xi * excess)./beta);
		y = y./xi;
		f = length(excess) * log(beta) + (1 + xi) * sum(y);
    end





warning on