hessipot.m

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

function c=hessipot(func,x,exceedances,span,threshold),
ep=0.0001;
eps=ep*x;
n=length(x);
m=zeros(n,n);
for i=1:n,
    for j=1:n,
        x1=x;
        x1(i)=x1(i)+eps(i);
        x1(j)=x1(j)+eps(j);
        x2=x;
        x2(i)=x2(i)+eps(i);
        x2(j)=x2(j)-eps(j);
        x3=x;
        x3(i)=x3(i)-eps(i);
        x3(j)=x3(j)+eps(j);
        x4=x;
        x4(i)=x4(i)-eps(i);
        x4(j)=x4(j)-eps(j);
        
        
        m(i,j)=eval(['(' func '(x1,exceedances,span,threshold)-' func '(x2,exceedances,span,threshold)-' func '(x3,exceedances,span,threshold)+' func '(x4,exceedances,span,threshold)' ')' '/(4*eps(' num2str(i) ')*eps(' num2str(j) '))' ]);
    end
end

        c=inv(m);
    function f=negloglik(theta,exceedances,span,threshold),
    if theta(2)<=0 | min(1+(theta(1)*(exceedances-theta(3)))/theta(2)<=0),
        f=NaN;
    else
        y=log(1+(theta(1)*(exceedances-theta(3)))/theta(2));
        term3=(1/theta(1)+1)*sum(y);
        term1=span*(1+(theta(1)*(threshold-theta(3)))/theta(2))^(-1/theta(1));
        term2=length(y)*log(theta(2));
        f=term1+term2+term3;
    end
		

        
%         "hess" <- function(f, x)
% {
% #a function by Stuart Coles
% 	ep <- 0.0001
% 	eps <- ep * x
% 	n <- length(x)
% 	m <- matrix(0, ncol = n, nrow = n)
% 	for(i in 1:n) {
% 		for(j in 1:n) {
% 			x1 <- x
% 			x1[i] <- x1[i] + eps[i]
% 			x1[j] <- x1[j] + eps[j]
% 			x2 <- x
% 			x2[i] <- x2[i] + eps[i]
% 			x2[j] <- x2[j] - eps[j]
% 			x3 <- x
% 			x3[i] <- x3[i] - eps[i]
% 			x3[j] <- x3[j] + eps[j]
% 			x4 <- x
% 			x4[i] <- x4[i] - eps[i]
% 			x4[j] <- x4[j] - eps[j]
% 			m[i, j] <- (f(x1) - f(x2) - f(x3) + f(x4))/(4 * 
% 				eps[i] * eps[j])
% 		}
% 	}
% 	solve(m)
% }
%