Home > wafo > wstats > wgpdpdf.m

# wgpdpdf

## PURPOSE

Generalized Pareto probability density function

## SYNOPSIS

f = wgpdpdf(x,k,s,m);

## DESCRIPTION

``` WGPDPDF Generalized Pareto probability density function

CALL:  f = wgpdpdf(x,k,s,m);

f = density 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) =
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 = wgpdpdf(x,1.25,1); p2 = wgpdpdf(x,1,1);
p3 = wgpdpdf(x,0.75,1); p4 = wgpdpdf(x,0.5,1);
plot(x,p1,x,p2,x,p3,x,p4)```

## CROSS-REFERENCE INFORMATION

This function calls:
 comnsize Check if all input arguments are either scalar or of common size. error Display message and abort function. nan Not-a-Number.
This function is called by:

## SOURCE CODE

```001 function f = wgpdpdf(x,k,s,m);
002 %WGPDPDF Generalized Pareto probability density function
003 %
004 % CALL:  f = wgpdpdf(x,k,s,m);
005 %
006 %        f = density 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) =
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 = wgpdpdf(x,1.25,1); p2 = wgpdpdf(x,1,1);
022 %   p3 = wgpdpdf(x,0.75,1); p4 = wgpdpdf(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 % -fixed a bug for k<0
034 % revised pab June 2005
035 % fixed bug for k==0 now fixed
036 % revised jr 14.08.2001
037 %  - a bug in the last if-statement condition fixed
038 %    (thanks to D Eddelbuettel <edd@debian.org>)
039 % revised pab 24.10.2000
040 %  - added comnsize, nargchk
041 %  - added extra check on the scale parameter s
042 %  - added m + default values on m and s
044
045 error(nargchk(2,4,nargin))
046
047 if nargin<4|isempty(m), m=0;end
048 if nargin<3|isempty(s), s=1;end
049
050 [errorcode x k s m] = comnsize(x,k,s,m);
051 if errorcode > 0
052     error('x, k s and m must be of common size or scalar.');
053 end
054 epsilon=1e-4; % treshold defining k to zero
055
056 xm = x-m;
057 f = zeros(size(x));
058
059 k0 = find((xm>=0) & (abs(k)<=epsilon) & (s>0));
060 if any(k0),
061    f(k0) = exp(-xm(k0)./s(k0))./s(k0);
062 end
063
064 k1=find((xm>=0) & (k.*xm < s) & (abs(k)>epsilon) & (s>0));
065 if any(k1),
066   f(k1) = (1-k(k1).*xm(k1)./s(k1)).^(1./k(k1)-1)./s(k1);
067 end
068
069
070 k3 = find(s<=0);
071 if any(k3),
072    tmp   = NaN;
073    f(k3) = tmp(ones(size(k3)));
074 end
075
076
077
078
079
080
081
082
083
084
085
086
087```

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