WGPDCDF Generalized Pareto cumulative distribution function CALL: F = wgpdcdf(x,k,s,m); F = distribution function evaluated at x k = shape parameter in the GPD s = scale parameter in the GPD (default 1) m = location parameter in the GPD (default 0) The Generalized Pareto distribution is defined by its cdf 1 - (1-k(x-m)/s)^1/k, k~=0 F(x;k,s,m) = 1 - exp(-(x-m)/s), k==0 for x>=m (when k<=0) and 0 <= x-m < s/k (when k>0), s>0. Example: x = linspace(0,2,200); p1 = wgpdcdf(x,1.25,1); p2 = wgpdcdf(x,1,1); p3 = wgpdcdf(x,0.75,1); p4 = wgpdcdf(x,0.5,1); plot(x,p1,x,p2,x,p3,x,p4)
Check if all input arguments are either scalar or of common size. | |
Display message and abort function. | |
Not-a-Number. |
% CHAPTER5 contains the commands used in Chapter 5 of the tutorial | |
Extrapolate level crossing spectrum | |
Extrapolate level crossing spectrum | |
Parameter estimates for Generalized Pareto data |
001 function F = wgpdcdf(x,k,s,m); 002 %WGPDCDF Generalized Pareto cumulative distribution function 003 % 004 % CALL: F = wgpdcdf(x,k,s,m); 005 % 006 % F = distribution function evaluated at x 007 % k = shape parameter in the GPD 008 % s = scale parameter in the GPD (default 1) 009 % m = location parameter in the GPD (default 0) 010 % 011 % The Generalized Pareto distribution is defined by its cdf 012 % 013 % 1 - (1-k(x-m)/s)^1/k, k~=0 014 % F(x;k,s,m) = 015 % 1 - exp(-(x-m)/s), k==0 016 % 017 % for x>=m (when k<=0) and 0 <= x-m < s/k (when k>0), s>0. 018 % 019 % Example: 020 % x = linspace(0,2,200); 021 % p1 = wgpdcdf(x,1.25,1); p2 = wgpdcdf(x,1,1); 022 % p3 = wgpdcdf(x,0.75,1); p4 = wgpdcdf(x,0.5,1); 023 % plot(x,p1,x,p2,x,p3,x,p4) 024 025 % References 026 % Johnson N.L., Kotz S. and Balakrishnan, N. (1994) 027 % Continuous Univariate Distributions, Volume 1. Wiley. 028 029 % Tested on: Matlab 5.3 030 % History: 031 % Revised pab oct 2005 032 % - limits m<x<s/k replaced with 0 <= x-m < s/k in help header. 033 % revised pab June 2005 034 % - distribution for k==0 was wrong, now fixed 035 % Revised by jr 22.12.1999 036 % Modified by PJ 08-Mar-2000 037 % Hjälptext 038 % revised ms 14.06.2000 039 % - updated header info 040 % - changed name to wgpdcdf (from gpdcdf) 041 % revised pab 24.10.2000 042 % - added nargchk, comnsize and default values for m, s 043 % revised jr 14.08.2001 044 % - a bug in the last if-statement condition fixed 045 % (thanks to D Eddelbuettel <edd@debian.org>) 046 047 error(nargchk(2,4,nargin)) 048 049 if nargin<4|isempty(m), m=0;end 050 if nargin<3|isempty(s), s=1;end 051 052 [errorcode x k s,m] = comnsize(x,k,s,m); 053 if errorcode > 0 054 error('x, k, s and m must be of common size or scalar.'); 055 end 056 057 epsilon = 1e-4; % treshold defining k to zero 058 059 060 F = zeros(size(x)); 061 xm = x-m; 062 k0 = find((0 < xm) & (abs(k)<=epsilon) & (0<s)); 063 if any(k0), 064 F(k0) = 1-exp(-xm(k0)./s(k0)); 065 end 066 067 k1 = find((0<xm) & (k.*xm < s ) & (abs(k)>epsilon) & (0<s)); 068 if any(k1), 069 F(k1) = 1-(1-k(k1).*xm(k1)./s(k1)).^(1./k(k1)); 070 end 071 072 k2 = find((k.*xm>=s) & (k>epsilon)); 073 if any(k2), 074 F(k2)=ones(size(k2)); 075 end 076 077 k3 = find(s<=0 ); 078 if any(k3), 079 tmp = NaN; 080 F(k3) = tmp(ones(size(k3))); 081 end 082 083 return 084 % old call 085 if abs(k) < epsilon, 086 p = (x>0).*(1-exp(-x/s)); 087 elseif k>0, 088 p=(x>=s/k)+(x>0).*(x<s/k).*... 089 (1-(1-k*x/s).^(1/k)); 090 else 091 p=zeros(size(x)); 092 p=(x>0).*(1-(1-k*x/s).^(1/k)); 093 end; 094 095 096 097 098
Comments or corrections to the WAFO group