WGPDFIT_MLD Parameter estimates for GPD data CALL: [f,k,s] = wgpdfit_mld(x,data) f = function values. k = shape parameter of GPD. s = scale parameter of GPD. This function is used by cmat2extralc for numerical solution of the ML estimate, i.e. solve f=0 for x. data0 = wgpdrnd(0.3,1,200,1); x_ML = fzero('wgpdfit_ml',0,[],data0); [f,k_ML,s_ML] = wgpdfit_ml(x_ML,data0) % Estimates k_ML and s_ML data1 = floor(data0*10)/10; x=(0:0.1:(max(data1)+0.1))'; N = histc(data1+0.05,x); x_MLD = fzero('wgpdfit_mld',0,[],[x N]); [f,k_MLD,s_MLD] = wgpdfit_mld(x_MLD,[x N]) % Estimates k_ML and s_ML See also wgpdfit, cmat2extralc
Extrapolate level crossing spectrum | |
Quick test of the routines in module 'cycles' |
001 function [f,k,s] = wgpdfit_mld(x,data) 002 %WGPDFIT_MLD Parameter estimates for GPD data 003 % 004 % CALL: [f,k,s] = wgpdfit_mld(x,data) 005 % 006 % f = function values. 007 % k = shape parameter of GPD. 008 % s = scale parameter of GPD. 009 % 010 % This function is used by cmat2extralc for numerical solution of 011 % the ML estimate, i.e. solve f=0 for x. 012 % data0 = wgpdrnd(0.3,1,200,1); 013 % x_ML = fzero('wgpdfit_ml',0,[],data0); 014 % [f,k_ML,s_ML] = wgpdfit_ml(x_ML,data0) % Estimates k_ML and s_ML 015 % data1 = floor(data0*10)/10; 016 % x=(0:0.1:(max(data1)+0.1))'; 017 % N = histc(data1+0.05,x); 018 % x_MLD = fzero('wgpdfit_mld',0,[],[x N]); 019 % [f,k_MLD,s_MLD] = wgpdfit_mld(x_MLD,[x N]) % Estimates k_ML and s_ML 020 % 021 % See also wgpdfit, cmat2extralc 022 023 % References 024 % 025 % Davidson & Smith (1990) 026 % Models for Exceedances over high Threholds. 027 % Journal of the Royal Statistical Society B,52, pp. 393-442. 028 029 % Tested on; Matlab 5.3 030 % History: 031 % Created by PJ 10-Oct-2000 032 % - created from wgpdfit_ml 033 034 % In order to avoid boundary problems in numerical solution we use a transformation 035 % Transformation: x = log(1/max_data - t), -Inf < t < 1/max_data 036 % Inverse Trans.: t = 1/max(data) - exp(x), -Inf < x < Inf 037 038 M = max(data(:,1)); % Max of data 039 t = 1/M - exp(x); % Inverse Transformation 040 041 N = sum(data(:,2)); 042 043 if t<1/M 044 k = -1/N*sum(data(:,2).*log(1-t*data(:,1))); % Shape parameter 045 s = k/t; % Scale parameter 046 % Evaluate function 047 f = (1/k-1)*sum(data(:,2).*data(:,1)./(1-t*data(:,1))) - N/t; 048 else 049 I = (data(:,1)==M); 050 k = -1/N*(sum(data(~I,2).*log(1-t*data(~I,1)))+sum(data(I,2).*(x+log(data(I,1))))); % Shape parameter 051 s = k/t; % Scale parameter 052 % Evaluate function 053 f = (1/k-1)*(sum(data(~I,2).*data(~I,1)./(1-t*data(~I,1)))+sum(data(~I,2)*exp(-x)))- N/t; 054 end 055 056 057 % Evaluate function 058 %f = (1/k-1)*sum(data(:,2).*data(:,1)./(1-t*data(:,1))) - N/t; 059 060 if isinf(f) 061 f=realmax*sign(f); 062 end 063 if isnan(f) %if x<0, f=realmax; else, f=-realmax; end 064 f=realmax; 065 end 066
Comments or corrections to the WAFO group