pot.m

来自「极值理论中各种函数及图像的程序。matlab实现。」· M 代码 · 共 90 行

function [res,data]=pot(data,threshold,nextremes,run),
%Fits a Poisson point process to the data, an approach sometimes known as peaks over
%thresholds (POT)
%
%USAGE: [res,data2]=pot(data,threshold,nextremes,run)
%
%Either threshold or nextremes should be supplied. The undefined one should be entered
%as [] !!!
%
%
%     data: Data vector.
%threshold: Threshold value. Either this or nextremes should be supplied. The undefined 
%           one should be entered as [].
%nextremes: Implies a threshold value that number of observations remaining above is
%           nextremes.
%      run: If the data is to be declustered enter the separation here else enter [] .
%
%      res: Fitted process.
%    data2: Data vector assigned time index from 0 to 1.
warning off
times=(1:length(data))/length(data);
data=surecol(data);
data=[times',data];

span=max(cumsum(diff(times)));
n=length(data);
first=times(1);
last=times(n);
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(:,2),nextremes);
end
if threshold > 10,
    factor=10^(floor(log10(threshold)));
end

exceedances=data(data(:,2)>threshold,:);
n_exceed=length(exceedances);
p_less_thresh=1-n_exceed/n;

if ~isempty(run),
    exceedances=data(data(:,2)>threshold,:); 
    exceedances=decluster(exceedances,run,1);
    %exceedances=exceedances(:,2);
    n_exceed=length(exceedances);
end

intensity=n_exceed/span;
xbar=mean(exceedances(:,2))-threshold;
s2=var(exceedances(:,2));
shape0=-0.5*(((xbar*xbar)/s2)-1);
extra=((length(exceedances)/span)^(-shape0)-1)/shape0;
betahat=0.5*xbar*((xbar*xbar/s2)+1);
scale0=betahat/(1+shape0*extra);
loc0=0;
theta=[shape0,scale0,loc0];
cond1=theta(2)<=0 ;
cond2=min(1+(theta(1)*(exceedances(:,2)-theta(3)))/theta(2)<=0);
cond3=imag((1+(theta(1)*(threshold-theta(3)))/theta(2))^(-1/theta(1)))~=0;
if cond1|cond2|cond3;
    theta(1)=0;
    theta(2)=0;
    theta(3)=0;
end

res.data=exceedances;
opts=optimset('MaxFunEvals',10000,'MaxIter',10000,'TolX',1e-7,'TolFun',1e-7,'Display','off');
%[res.par_ests,res.funval,res.terminated,res.details] = fminsearch('negloglikpot',theta,opts,exceedances(:,2),span,threshold);
[res.par_ests,res.funval,res.terminated,res.details] = fminunc('negloglikpot',theta,opts,exceedances(:,2),span,threshold);
%[res.par_ests,res.funval,res.terminated,res.details] = fminsearch('negloglikpot',res.par_ests,opts,exceedances(:,2),span,threshold);
res.varcov=hessipot('negloglikpot',res.par_ests,exceedances(:,2),span,threshold);
res.par_ses=sqrt(diag(res.varcov))';
beta = res.par_ests(2) + res.par_ests(1) * (threshold - res.par_ests(3));
res.par_ests=[res.par_ests beta];
res.intensity=intensity;
res.threshold=threshold;
res.span=span;
res.p_less_thresh=p_less_thresh;
res.n_exceed=n_exceed;
warning on