Home > wafo > cycles > wgpdfit_mld.m

wgpdfit_mld

PURPOSE ^

Parameter estimates for GPD data

SYNOPSIS ^

[f,k,s] = wgpdfit_mld(x,data)

DESCRIPTION ^

 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

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

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

Mathematical Statistics
Centre for Mathematical Sciences
Lund University with Lund Institute of Technology

Comments or corrections to the WAFO group


Generated on Thu 06-Oct-2005 02:21:16 for WAFO by m2html © 2003