Mercurial > forge
changeset 11336:a7de0aafd809 octave-forge
new gev* functions added to statistics
| author | nir-krakauer |
|---|---|
| date | Fri, 28 Dec 2012 00:04:09 +0000 |
| parents | a9cdb3c433f4 |
| children | 3f7da417c367 |
| files | main/statistics/INDEX main/statistics/NEWS main/statistics/inst/gevcdf.m main/statistics/inst/gevfit.m main/statistics/inst/gevinv.m main/statistics/inst/gevlike.m main/statistics/inst/gevpdf.m main/statistics/inst/gevrnd.m main/statistics/inst/gevstat.m |
| diffstat | 9 files changed, 943 insertions(+), 1 deletions(-) [+] |
line wrap: on
line diff
--- a/main/statistics/INDEX Thu Dec 27 22:11:24 2012 +0000 +++ b/main/statistics/INDEX Fri Dec 28 00:04:09 2012 +0000 @@ -11,6 +11,7 @@ gamlike gamstat geostat + gevcdf gevfit gevinv gevlike gevpdf gevrnd gevstat hygestat jsucdf jsupdf
--- a/main/statistics/NEWS Thu Dec 27 22:11:24 2012 +0000 +++ b/main/statistics/NEWS Fri Dec 28 00:04:09 2012 +0000 @@ -4,7 +4,11 @@ ** The following functions are new: regress_gp dendogram plsregress - + + ** New functions for the generalized extreme value (GEV) distribution: + + gevcdf gevfit gevinv gevlike gevpdf gevrnd gevstat + ** The interface of the following functions has been modified: mvnrnd
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/main/statistics/inst/gevcdf.m Fri Dec 28 00:04:09 2012 +0000 @@ -0,0 +1,131 @@ +## Copyright (C) 2012 Nir Krakauer <nkrakauer@ccny.cuny.edu> +## +## This program is free software; you can redistribute it and/or modify it under +## the terms of the GNU General Public License as published by the Free Software +## Foundation; either version 3 of the License, or (at your option) any later +## version. +## +## This program is distributed in the hope that it will be useful, but WITHOUT +## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more +## details. +## +## You should have received a copy of the GNU General Public License along with +## this program; if not, see <http://www.gnu.org/licenses/>. + +## -*- texinfo -*- +## @deftypefn {Function File} {@var{p} =} gevcdf (@var{x}, @var{k}, @var{sigma}, @var{mu}) +## Compute the cumulative distribution function of the generalized extreme value (GEV) distribution. +## +## @subheading Arguments +## +## @itemize @bullet +## @item +## @var{x} is the support. +## +## @item +## @var{k} is the shape parameter of the GEV distribution. (Also denoted gamma or xi.) +## @item +## @var{sigma} is the scale parameter of the GEV distribution. The elements +## of @var{sigma} must be positive. +## @item +## @var{mu} is the location parameter of the GEV distribution. +## @end itemize +## The inputs must be of common size, or some of them must be scalar. +## +## @subheading Return values +## +## @itemize @bullet +## @item +## @var{p} is the cumulative distribution of the GEV distribution at each +## element of @var{x} and corresponding parameter values. +## @end itemize +## +## @subheading Examples +## +## @example +## @group +## x = 0:0.5:2.5; +## sigma = 1:6; +## k = 1; +## mu = 0; +## y = gevcdf (x, k, sigma, mu) +## @end group +## +## @group +## y = gevcdf (x, k, 0.5, mu) +## @end group +## @end example +## +## @subheading References +## +## @enumerate +## @item +## Rolf-Dieter Reiss and Michael Thomas. @cite{Statistical Analysis of Extreme Values with Applications to Insurance, Finance, Hydrology and Other Fields}. Chapter 1, pages 16-17, Springer, 2007. +## +## @end enumerate +## @seealso{gevfit, gevinv, gevlike, gevpdf, gevrnd, gevstat} +## @end deftypefn + +## Author: Nir Krakauer <nkrakauer@ccny.cuny.edu> +## Description: CDF of the generalized extreme value distribution + +function p = gevcdf (x, k, sigma, mu) + + # Check arguments + if (nargin != 4) + print_usage (); + endif + + if (isempty (x) || isempty (k) || isempty (sigma) || isempty (mu) || ~ismatrix (x) || ~ismatrix (k) || ~ismatrix (sigma) || ~ismatrix (mu)) + error ("gevcdf: inputs must be a numeric matrices"); + endif + + [retval, x, k, sigma, mu] = common_size (x, k, sigma, mu); + if (retval > 0) + error ("gevcdf: inputs must be of common size or scalars"); + endif + + z = 1 + k .* (x - mu) ./ sigma; + + # Calculate pdf + p = exp(-(z .^ (-1 ./ k))); + + p(z <= 0 & x < mu) = 0; + p(z <= 0 & x > mu) = 1; + + inds = (k == 0); %use a different formula + if any(inds) + z = (mu(inds) - x(inds)) ./ sigma(inds); + p(inds) = exp(-exp(z)); + endif + + +endfunction + +%!test +%! x = 0:0.5:2.5; +%! sigma = 1:6; +%! k = 1; +%! mu = 0; +%! p = gevcdf (x, k, sigma, mu); +%! expected_p = [0.36788 0.44933 0.47237 0.48323 0.48954 0.49367]; +%! assert (p, expected_p, 0.001); + +%!test +%! x = -0.5:0.5:2.5; +%! sigma = 0.5; +%! k = 1; +%! mu = 0; +%! p = gevcdf (x, k, sigma, mu); +%! expected_p = [0 0.36788 0.60653 0.71653 0.77880 0.81873 0.84648]; +%! assert (p, expected_p, 0.001); + +%!test #check for continuity for k near 0 +%! x = 1; +%! sigma = 0.5; +%! k = -0.03:0.01:0.03; +%! mu = 0; +%! p = gevcdf (x, k, sigma, mu); +%! expected_p = [0.88062 0.87820 0.87580 0.87342 0.87107 0.86874 0.86643]; +%! assert (p, expected_p, 0.001);
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/main/statistics/inst/gevfit.m Fri Dec 28 00:04:09 2012 +0000 @@ -0,0 +1,85 @@ +## Copyright (C) 2012 Nir Krakauer <nkrakauer@ccny.cuny.edu> +## +## This program is free software; you can redistribute it and/or modify it under +## the terms of the GNU General Public License as published by the Free Software +## Foundation; either version 3 of the License, or (at your option) any later +## version. +## +## This program is distributed in the hope that it will be useful, but WITHOUT +## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more +## details. +## +## You should have received a copy of the GNU General Public License along with +## this program; if not, see <http://www.gnu.org/licenses/>. + +## -*- texinfo -*- +## @deftypefn {Function File} {@var{paramhat}, @var{paramci} =} gevfit (@var{data}) +## Find the maximum likelihood estimator (@var{paramhat}) of the generalized extreme value (GEV) distribution to fit @var{data}. +## +## @subheading Arguments +## +## @itemize @bullet +## @item +## @var{data} is the vector of given values. +## +## @end itemize +## +## @subheading Return values +## +## @itemize @bullet +## @item +## @var{parmhat} is the 3-parameter maximum-likelihood parameter vector [@var{k}, @var{sigma}, @var{mu}], where @var{k} is the shape parameter of the GEV distribution, @var{sigma} is the scale parameter of the GEV distribution, and @var{mu} is the location parameter of the GEV distribution. +## @item +## @var{paramci} has the approximate 95% confidence intervals of the parameter values based on the Fisher information matrix at the maximum-likelihood position. +## +## @end itemize +## +## @subheading Examples +## +## @example +## @group +## data = 1:50; +## [pfit, pci] = gevfit (data); +## p1 = gevcdf(data,pfit(1),pfit(2),pfit(3)); +## plot(data, p1) +## @end group +## @end example +## @seealso{gevcdf, gevinv, gevlike, gevpdf, gevrnd, gevstat} +## @end deftypefn + +## Author: Nir Krakauer <nkrakauer@ccny.cuny.edu> +## Description: Maximum likelihood parameter estimation for the generalized extreme value distribution + +function [paramhat, paramci] = gevfit (data) + + # Check arguments + if (nargin < 1) + print_usage; + endif + + #initial guess for parameter values + p0 = [0 1 0]'; + + #cost function to minimize + f = @(p) gevlike(p, data); + + paramhat = fminunc(f, p0); + + if nargout > 1 + [nlogL, ACOV] = gevlike (paramhat, data); + param_se = sqrt(diag(-ACOV)); + paramci(:, 1) = paramhat - 1.96*param_se; + paramci(:, 2) = paramhat + 1.96*param_se; + endif + +endfunction + + + +%!test +%! data = 1:50; +%! [pfit, pci] = gevfit (data); +%! sigma = 0.3; +%! expected_p = [-0.44 15.19 21.53]'; +%! assert (pfit, expected_p, 0.1);
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/main/statistics/inst/gevinv.m Fri Dec 28 00:04:09 2012 +0000 @@ -0,0 +1,88 @@ +## Copyright (C) 2012 Nir Krakauer <nkrakauer@ccny.cuny.edu> +## +## This program is free software; you can redistribute it and/or modify it under +## the terms of the GNU General Public License as published by the Free Software +## Foundation; either version 3 of the License, or (at your option) any later +## version. +## +## This program is distributed in the hope that it will be useful, but WITHOUT +## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more +## details. +## +## You should have received a copy of the GNU General Public License along with +## this program; if not, see <http://www.gnu.org/licenses/>. + +## -*- texinfo -*- +## @deftypefn {Function File} {@var{X} =} gevinv (@var{P}, @var{k}, @var{sigma}, @var{mu}) +## Compute a desired quantile (inverse CDF) of the generalized extreme value (GEV) distribution. +## +## @subheading Arguments +## +## @itemize @bullet +## @item +## @var{P} is the desired quantile of the GEV distribution. (Between 0 and 1.) +## @item +## @var{k} is the shape parameter of the GEV distribution. (Also denoted gamma or xi.) +## @item +## @var{sigma} is the scale parameter of the GEV distribution. The elements +## of @var{sigma} must be positive. +## @item +## @var{mu} is the location parameter of the GEV distribution. +## @end itemize +## The inputs must be of common size, or some of them must be scalar. +## +## @subheading Return values +## +## @itemize @bullet +## @item +## @var{X} is the value corresponding to each quantile of the GEV distribution +## @end itemize +## @subheading References +## +## @enumerate +## @item +## Rolf-Dieter Reiss and Michael Thomas. @cite{Statistical Analysis of Extreme Values with Applications to Insurance, Finance, Hydrology and Other Fields}. Chapter 1, pages 16-17, Springer, 2007. +## @item +## J. R. M. Hosking (2012). @cite{L-moments}. R package, version 1.6. URL: http://CRAN.R-project.org/package=lmom. +## +## @end enumerate +## @seealso{gevcdf, gevfit, gevlike, gevpdf, gevrnd, gevstat} +## @end deftypefn + +## Author: Nir Krakauer <nkrakauer@ccny.cuny.edu> +## Description: Inverse CDF of the generalized extreme value distribution + +function [X] = gevinv (P, k = 0, sigma = 1, mu = 0) + + [retval, P, k, sigma, mu] = common_size (P, k, sigma, mu); + if (retval > 0) + error ("gevinv: inputs must be of common size or scalars"); + endif + + X = P; + + ii = (k == 0); + X(ii) = mu(ii) - sigma(ii) .* log(-log(P(ii))); + X(~ii) = mu(~ii) + sigma(~ii) .* (-1 - (log(P(~ii)).^ -k(~ii))) ./ k(~ii); + +endfunction + +%!test +%! p = 0.1:0.1:0.9; +%! k = 0; +%! sigma = 1; +%! mu = 0; +%! x = gevinv (p, k, sigma, mu); +%! c = gevcdf(x, k, sigma, mu); +%! assert (c, p, 0.001); + +%!test +%! p = 0.1:0.1:0.9; +%! k = 1; +%! sigma = 1; +%! mu = 0; +%! x = gevinv (p, k, sigma, mu); +%! c = gevcdf(x, k, sigma, mu); +%! assert (c, p, 0.001); +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/main/statistics/inst/gevlike.m Fri Dec 28 00:04:09 2012 +0000 @@ -0,0 +1,292 @@ +## Copyright (C) 2012 Nir Krakauer <nkrakauer@ccny.cuny.edu> +## +## This program is free software; you can redistribute it and/or modify it under +## the terms of the GNU General Public License as published by the Free Software +## Foundation; either version 3 of the License, or (at your option) any later +## version. +## +## This program is distributed in the hope that it will be useful, but WITHOUT +## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more +## details. +## +## You should have received a copy of the GNU General Public License along with +## this program; if not, see <http://www.gnu.org/licenses/>. + +## -*- texinfo -*- +## @deftypefn {Function File} {@var{nlogL}, @var{ACOV}, @var{Grad} =} gevlike (@var{params}, @var{data}) +## Compute the negative log-likelihood of data under the generalized extreme value (GEV) distribution with given parameter values. +## +## @subheading Arguments +## +## @itemize @bullet +## @item +## @var{params} is the 3-parameter vector [@var{k}, @var{sigma}, @var{mu}], where @var{k} is the shape parameter of the GEV distribution, @var{sigma} is the scale parameter of the GEV distribution, and @var{mu} is the location parameter of the GEV distribution. +## @item +## @var{data} is the vector of given values. +## +## @end itemize +## +## @subheading Return values +## +## @itemize @bullet +## @item +## @var{nlogL} is the negative log-likelihood. +## @item +## @var{ACOV} is the 3 by 3 Fisher information matrix (negative second derivative of the log likelihood with respect to the parameter values) +## @var{Grad} is the 3 by 1 gradient vector (first derivative of the log likelihood with respect to the parameter values) +## +## @end itemize +## +## @subheading Examples +## +## @example +## @group +## x = -5:-1; +## k = -0.2; +## sigma = 0.3; +## mu = 0.5; +## [L, C] = gevlike ([k sigma mu], x); +## @end group +## @end example +## +## @subheading References +## +## @enumerate +## @item +## Rolf-Dieter Reiss and Michael Thomas. @cite{Statistical Analysis of Extreme Values with Applications to Insurance, Finance, Hydrology and Other Fields}. Chapter 1, pages 16-17, Springer, 2007. +## +## @end enumerate +## @seealso{gevcdf, gevfit, gevinv, gevpdf, gevrnd, gevstat} +## @end deftypefn + +## Author: Nir Krakauer <nkrakauer@ccny.cuny.edu> +## Description: Negative log-likelihood for the generalized extreme value distribution + +function [nlogL, ACOV, Grad] = gevlike (params, data) + + # Check arguments + if (nargin != 2) + print_usage; + endif + + k = params(1); + sigma = params(2); + mu = params(3); + + #calculate negative log likelihood + nlogL = sum(gevnll (data, k, sigma, mu) (:)); + + if nargout > 1 + ACOV = -gevfim (data, k, sigma, mu); + if nargout > 2 + Grad = gevgrad (data, k, sigma, mu); + endif + endif + +endfunction + + +function nlogL = gevnll (x, k, sigma, mu) +#internal function to calculate negative log likelihood for gevlike +#no input checking done + + if k == 0 + z = (mu - x) ./ sigma; + nlogL = exp(z) - z + log(sigma); + else + z = 1 + k .* (x - mu) ./ sigma; + nlogL = z .^ (-1 ./ k) + (1 + 1 ./ k) .* log(z) + log(sigma); + nlogL(z <= 0) = Inf; + endif + +endfunction + +function G = gevgrad (x, k, sigma, mu) +#calculate the gradient of the negative log likelihood of data x with respect to the parameters of the generalized extreme value distribution for gevlike +#no input checking done + +G = ones(3, 1); + +if k == 0 ##use the expressions for first derivatives that are the limits as k --> 0 + a = (x - mu) ./ sigma; + f = exp(a); + #k + g = -(2 * x .* (mu .* (1 - f) - sigma .* f) + 2 .* sigma .* mu .* f + (x.^2 + mu.^2).*(f - 1)) ./ (2 * f .* sigma .^ 2); + G(1) = sum(g(:)); + + #sigma + g = ((a ./ f) - a + 1) ./ sigma; + G(2) = sum(g(:)); + + #mu + g = (1 ./ f - 1) ./ sigma; + G(3) = sum(g(:)); + + return +endif + +z = 1 + k .* (x - mu) ./ sigma; +if any (z <= 0) + G(:) = 0; #negative log likelihood is locally infinite + return +end + +#k +a = (x - mu) ./ sigma; +b = k .* a + 1; +c = log(b); +d = 1 ./ k + 1; +g = (c ./ k - a ./ b) ./ (k .* b .^ (1/k)) - c ./ (k .^ 2) + a .* d ./ b; +G(1) = sum(g(:)); + +#sigma +g = (a .* b .^ (-d) - d .* k .* a ./ b + 1) ./ sigma; +G(2) = sum(g(:)); + +#mu +g = (b .^ (-d) - d .* k ./ b) ./ sigma; +G(3) = sum(g(:)); + + +endfunction + +function ACOV = gevfim (x, k, sigma, mu) +#internal function to calculate the (negative) Fisher information matrix for gevlike +#no input checking done + +#calculate the various second derivatives (used Maxima to find the expressions) + +ACOV = ones(3); + +if k == 0 ##use the expressions for second derivatives that are the limits as k --> 0 + #k, k + a = (x - mu) ./ sigma; + f = exp(a); + der = (x .* (24 * mu .^ 2 .* sigma .* (f - 1) + 24 * mu .* sigma .^ 2 .* f - 12 * mu .^ 3) + x .^ 3 .* (8 * sigma .* (f - 1) - 12*mu) + x .^ 2 .* (-12 * sigma .^ 2 .* f + 24 * mu .* sigma .* (1 - f) + 18 * mu .^ 2) - 12 * mu .^ 2 .* sigma .^ 2 .* f + 8 * mu .^ 3 .* sigma .* (1 - f) + 3 * (x .^ 4 + mu .^ 4)) ./ (12 .* f .* sigma .^ 4); + ACOV(1, 1) = sum(der(:)); + + #sigma, sigma + der = (sigma .^ -2) .* ... + ((-2*a ./ f) + ... + a .^ 2 ./ f + ... + 2*a - 1); + ACOV(2, 2) = sum(der(:)); + + #mu, mu + der = (sigma .^ -2) ./ f; + ACOV(3, 3) = sum(der(:)); + + #k, sigma + der = (x .^2 .* (2*sigma .* (f - 1) - 3*mu) + ... + x .* (-2 * sigma .^ 2 .* f + 4 * mu .* sigma .* (1 - f) + 3 .* mu .^ 2) + ... + 2 * mu .^ 2 .* sigma .* (f - 1) + 2 * mu * sigma .^ 2 * f + x .^ 3 - mu .^ 3) ... + ./ (2 .* f .* sigma .^ 4); + ACOV(1, 2) = ACOV(2, 1) = sum(der(:)); + + #k, mu + der = (x .* (2*sigma .* (f - 1) - 2*mu) - 2 * f .* sigma .^ 2 + ... + 2 .* mu .* sigma .* (1 - f) + x .^ 2 + mu .^ 2)... + ./ (2 .* f .* sigma .^ 3); + ACOV(1, 3) = ACOV(3, 1) = sum(der(:)); + + #sigma, mu + der = (sigma + (x - sigma - mu) ./ f) ./ (sigma .^ 3); + ACOV(2, 3) = ACOV(3, 2) = sum(der(:)); + + return +endif + + +#general case + +z = 1 + k .* (x - mu) ./ sigma; + +if any (z <= 0) + ACOV(:) = 0; #negative log likelihood is locally infinite + return +end + +#k, k +a = (x - mu) ./ sigma; +b = k .* a + 1; +c = log(b); +d = 1 ./ k + 1; +der = ((((c ./ k.^2) - (a ./ (k .* b))) .^ 2) ./ (b .^ (1 ./ k))) + ... + ((-2*c ./ k.^3) + (2*a ./ (k.^2 .* b)) + ((a ./ b) .^ 2 ./ k)) ./ (b .^ (1 ./ k)) + ... + 2*c ./ k.^3 - ... + (2*a ./ (k.^2 .* b)) - (d .* (a ./ b) .^ 2); +ACOV(1, 1) = sum(der(:)); + +#sigma, sigma +der = (sigma .^ -2) .* (... + -2*a .* b .^ (-d) + ... + d .* k .* a .^ 2 .* (b .^ (-d-1)) + ... + 2 .* d .* k .* a ./ b - ... + d .* (k .* a ./ b) .^ 2 - 1); +ACOV(2, 2) = sum(der(:)); + +#mu, mu +der = (d .* (sigma .^ -2)) .* (... + k .* (b .^ (-d-1)) - ... + (k ./ b) .^ 2); +ACOV(3, 3) = sum(der(:)); + +#k, sigma +der = (a ./ sigma) .* (... +(b .^ (-d)) .* (c ./ k - a ./ b) ./ k + ... +-a .* (b .^ (-d-1)) + ... +((1 ./ k) - d) ./ b + ... +a .* k .* d ./ (b .^ 2)); +ACOV(1, 2) = ACOV(2, 1) = sum(der(:)); + +#k, mu +der = ( (b .^ (-d)) .* (c ./ k - a ./ b) ./ k - ... +a .* (b .^ (-d-1)) + ... +((1 ./ k) - d) ./ b + +a .* k .* d ./ (b .^ 2)) ./ sigma; +ACOV(1, 3) = ACOV(3, 1) = sum(der(:)); + +#sigma, mu +der = ( -(b .^ (-d)) + ... +a .* k .* d .* (b .^ (-d-1)) + ... +(d .* k ./ b) - a .* (k./b).^2 .* d) ./ (sigma .^ 2); +ACOV(2, 3) = ACOV(3, 2) = sum(der(:)); + +endfunction + + + + +%!test +%! x = 1; +%! k = 0.2; +%! sigma = 0.3; +%! mu = 0.5; +%! [L, C] = gevlike ([k sigma mu], x); +%! expected_L = 0.75942; +%! expected_C = -[-0.12547 1.77884 1.06731; 1.77884 16.40761 8.48877; 1.06731 8.48877 0.27979]; +%! assert (L, expected_L, 0.001); +%! assert (C, expected_C, 0.001); + +%!test +%! x = 1; +%! k = 0; +%! sigma = 0.3; +%! mu = 0.5; +%! [L, C] = gevlike ([k sigma mu], x); +%! expected_L = 0.65157; +%! expected_C = -[0.090036 3.41229 2.047337; 3.412229 24.760027 12.510190; 2.047337 12.510190 2.098618]; +%! assert (L, expected_L, 0.001); +%! assert (C, expected_C, 0.001); + +%!test +%! x = -5:-1; +%! k = -0.2; +%! sigma = 0.3; +%! mu = 0.5; +%! [L, C] = gevlike ([k sigma mu], x); +%! expected_L = 3786.4; +%! expected_C = [-1.4937e+06 1.0083e+06 -6.1837e+04; 1.0083e+06 -8.1138e+05 4.0917e+04; -6.1837e+04 4.0917e+04 -2.0422e+03]; +%! assert (L, expected_L, -0.001); +%! assert (C, expected_C, -0.001);
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/main/statistics/inst/gevpdf.m Fri Dec 28 00:04:09 2012 +0000 @@ -0,0 +1,130 @@ +## Copyright (C) 2012 Nir Krakauer <nkrakauer@ccny.cuny.edu> +## +## This program is free software; you can redistribute it and/or modify it under +## the terms of the GNU General Public License as published by the Free Software +## Foundation; either version 3 of the License, or (at your option) any later +## version. +## +## This program is distributed in the hope that it will be useful, but WITHOUT +## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more +## details. +## +## You should have received a copy of the GNU General Public License along with +## this program; if not, see <http://www.gnu.org/licenses/>. + +## -*- texinfo -*- +## @deftypefn {Function File} {@var{y} =} gevpdf (@var{x}, @var{k}, @var{sigma}, @var{mu}) +## Compute the probability density function of the generalized extreme value (GEV) distribution. +## +## @subheading Arguments +## +## @itemize @bullet +## @item +## @var{x} is the support. +## +## @item +## @var{k} is the shape parameter of the GEV distribution. (Also denoted gamma or xi.) +## @item +## @var{sigma} is the scale parameter of the GEV distribution. The elements +## of @var{sigma} must be positive. +## @item +## @var{mu} is the location parameter of the GEV distribution. +## @end itemize +## The inputs must be of common size, or some of them must be scalar. +## +## @subheading Return values +## +## @itemize @bullet +## @item +## @var{y} is the probability density of the GEV distribution at each +## element of @var{x} and corresponding parameter values. +## @end itemize +## +## @subheading Examples +## +## @example +## @group +## x = 0:0.5:2.5; +## sigma = 1:6; +## k = 1; +## mu = 0; +## y = gevpdf (x, k, sigma, mu) +## @end group +## +## @group +## y = gevpdf (x, k, 0.5, mu) +## @end group +## @end example +## +## @subheading References +## +## @enumerate +## @item +## Rolf-Dieter Reiss and Michael Thomas. @cite{Statistical Analysis of Extreme Values with Applications to Insurance, Finance, Hydrology and Other Fields}. Chapter 1, pages 16-17, Springer, 2007. +## +## @end enumerate +## @seealso{gevcdf, gevfit, gevinv, gevlike, gevrnd, gevstat} +## @end deftypefn + +## Author: Nir Krakauer <nkrakauer@ccny.cuny.edu> +## Description: PDF of the generalized extreme value distribution + +function y = gevpdf (x, k, sigma, mu) + + # Check arguments + if (nargin != 4) + print_usage (); + endif + + if (isempty (x) || isempty (k) || isempty (sigma) || isempty (mu) || ~ismatrix (x) || ~ismatrix (k) || ~ismatrix (sigma) || ~ismatrix (mu)) + error ("gevpdf: inputs must be a numeric matrices"); + endif + + [retval, x, k, sigma, mu] = common_size (x, k, sigma, mu); + if (retval > 0) + error ("gevpdf: inputs must be of common size or scalars"); + endif + + z = 1 + k .* (x - mu) ./ sigma; + + # Calculate pdf + y = exp(-(z .^ (-1 ./ k))) .* (z .^ (-1 - 1 ./ k)) ./ sigma; + + y(z <= 0) = 0; + + inds = (k == 0); %use a different formula + if any(inds) + z = (mu(inds) - x(inds)) ./ sigma(inds); + y(inds) = exp(z-exp(z)) ./ sigma(inds); + endif + + +endfunction + +%!test +%! x = 0:0.5:2.5; +%! sigma = 1:6; +%! k = 1; +%! mu = 0; +%! y = gevpdf (x, k, sigma, mu); +%! expected_y = [0.367879 0.143785 0.088569 0.063898 0.049953 0.040997]; +%! assert (y, expected_y, 0.001); + +%!test +%! x = -0.5:0.5:2.5; +%! sigma = 0.5; +%! k = 1; +%! mu = 0; +%! y = gevpdf (x, k, sigma, mu); +%! expected_y = [0 0.735759 0.303265 0.159229 0.097350 0.065498 0.047027]; +%! assert (y, expected_y, 0.001); + +%!test #check for continuity for k near 0 +%! x = 1; +%! sigma = 0.5; +%! k = -0.03:0.01:0.03; +%! mu = 0; +%! y = gevpdf (x, k, sigma, mu); +%! expected_y = [0.23820 0.23764 0.23704 0.23641 0.23576 0.23508 0.23438]; +%! assert (y, expected_y, 0.001);
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/main/statistics/inst/gevrnd.m Fri Dec 28 00:04:09 2012 +0000 @@ -0,0 +1,121 @@ +## Copyright (C) 2012 Nir Krakauer <nkrakauer@ccny.cuny.edu> +## +## This program is free software; you can redistribute it and/or modify it under +## the terms of the GNU General Public License as published by the Free Software +## Foundation; either version 3 of the License, or (at your option) any later +## version. +## +## Octave is distributed in the hope that it will be useful, but +## WITHOUT ANY WARRANTY; without even the implied warranty of +## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +## General Public License for more details. +## +## You should have received a copy of the GNU General Public License +## along with Octave; see the file COPYING. If not, see +## <http://www.gnu.org/licenses/>. + +## -*- texinfo -*- +## @deftypefn {Function File} {} gevrnd (@var{k}, @var{sigma}, @var{mu}) +## @deftypefnx {Function File} {} gevrnd (@var{k}, @var{sigma}, @var{mu}, @var{r}) +## @deftypefnx {Function File} {} gevrnd (@var{k}, @var{sigma}, @var{mu}, @var{r}, @var{c}, @dots{}) +## @deftypefnx {Function File} {} gevrnd (@var{k}, @var{sigma}, @var{mu}, [@var{sz}]) +## Return a matrix of random samples from the generalized extreme value (GEV) distribution with parameters +## @var{k}, @var{sigma}, @var{mu}. +## +## When called with a single size argument, returns a square matrix with +## the dimension specified. When called with more than one scalar argument the +## first two arguments are taken as the number of rows and columns and any +## further arguments specify additional matrix dimensions. The size may also +## be specified with a vector @var{sz} of dimensions. +## +## If no size arguments are given, then the result matrix is the common size of +## the input parameters. +## @seealso{gevcdf, gevfit, gevinv, gevlike, gevpdf, gevstat} +## @end deftypefn + +## Author: Nir Krakauer <nkrakauer@ccny.cuny.edu> +## Description: Random deviates from the generalized extreme value distribution + +function rnd = gevrnd (k, sigma, mu, varargin) + + if (nargin < 3) + print_usage (); + endif + + if any (sigma <= 0) + error ("gevrnd: sigma must be positive"); + endif + + if (!isscalar (k) || !isscalar (sigma) || !isscalar (mu)) + [retval, k, sigma, mu] = common_size (k, sigma, mu); + if (retval > 0) + error ("gevrnd: k, sigma, mu must be of common size or scalars"); + endif + endif + + if (iscomplex (k) || iscomplex (sigma) || iscomplex (mu)) + error ("gevrnd: k, sigma, mu must not be complex"); + endif + + if (nargin == 3) + sz = size (k); + elseif (nargin == 4) + if (isscalar (varargin{1}) && varargin{1} >= 0) + sz = [varargin{1}, varargin{1}]; + elseif (isrow (varargin{1}) && all (varargin{1} >= 0)) + sz = varargin{1}; + else + error ("gevrnd: dimension vector must be row vector of non-negative integers"); + endif + elseif (nargin > 4) + if (any (cellfun (@(x) (!isscalar (x) || x < 0), varargin))) + error ("gevrnd: dimensions must be non-negative integers"); + endif + sz = [varargin{:}]; + endif + + if (!isscalar (k) && !isequal (size (k), sz)) + error ("gevrnd: k, sigma, mu must be scalar or of size SZ"); + endif + + if (isa (k, "single") || isa (sigma, "single") || isa (mu, "single")) + cls = "single"; + else + cls = "double"; + endif + + rnd = gevinv (rand(sz), k, sigma, mu); + + if (strcmp (cls, "single")) + rnd = single (rnd); + endif + +endfunction + + +%!assert(size (gevrnd (1,2,1)), [1, 1]); +%!assert(size (gevrnd (ones(2,1), 2, 1)), [2, 1]); +%!assert(size (gevrnd (ones(2,2), 2, 1)), [2, 2]); +%!assert(size (gevrnd (1, 2*ones(2,1), 1)), [2, 1]); +%!assert(size (gevrnd (1, 2*ones(2,2), 1)), [2, 2]); +%!assert(size (gevrnd (1, 2, 1, 3)), [3, 3]); +%!assert(size (gevrnd (1, 2, 1, [4 1])), [4, 1]); +%!assert(size (gevrnd (1, 2, 1, 4, 1)), [4, 1]); + +%% Test input validation +%!error gevrnd () +%!error gevrnd (1, 2) +%!error gevrnd (ones(3),ones(2),1) +%!error gevrnd (ones(2),ones(3),1) +%!error gevrnd (i, 2, 1) +%!error gevrnd (2, i, 1) +%!error gevrnd (2, 0, 1) +%!error gevrnd (1,2, 1, -1) +%!error gevrnd (1,2, 1, ones(2)) +%!error gevrnd (1,2, 1, [2 -1 2]) +%!error gevrnd (1,2, 1, 1, ones(2)) +%!error gevrnd (1,2, 1, 1, -1) +%!error gevrnd (ones(2,2), 2, 1, 3) +%!error gevrnd (ones(2,2), 2, 1, [3, 2]) +%!error gevrnd (ones(2,2), 2, 1, 2, 3) +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/main/statistics/inst/gevstat.m Fri Dec 28 00:04:09 2012 +0000 @@ -0,0 +1,90 @@ +## Copyright (C) 2012 Nir Krakauer <nkrakauer@ccny.cuny.edu> +## +## This program is free software; you can redistribute it and/or modify it under +## the terms of the GNU General Public License as published by the Free Software +## Foundation; either version 3 of the License, or (at your option) any later +## version. +## +## This program is distributed in the hope that it will be useful, but WITHOUT +## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more +## details. +## +## You should have received a copy of the GNU General Public License along with +## this program; if not, see <http://www.gnu.org/licenses/>. + +## -*- texinfo -*- +## @deftypefn {Function File} {@var{y} =} gevstat (@var{k}, @var{sigma}, @var{mu}) +## Compute the mean and variance of the generalized extreme value (GEV) distribution. +## +## @subheading Arguments +## +## @itemize @bullet +## @item +## @var{k} is the shape parameter of the GEV distribution. (Also denoted gamma or xi.) +## @item +## @var{sigma} is the scale parameter of the GEV distribution. The elements +## of @var{sigma} must be positive. +## @item +## @var{mu} is the location parameter of the GEV distribution. +## @end itemize +## The inputs must be of common size, or some of them must be scalar. +## +## @subheading Return values +## +## @itemize @bullet +## @item +## @var{m} is the mean of the GEV distribution +## +## @item +## @var{v} is the variance of the GEV distribution +## @end itemize +## @seealso{gevcdf, gevfit, gevinv, gevlike, gevpdf, gevrnd} +## @end deftypefn + +## Author: Nir Krakauer <nkrakauer@ccny.cuny.edu> +## Description: Moments of the generalized extreme value distribution + +function [m, v] = gevstat (k, sigma, mu) + + # Check arguments + if (nargin < 3) + print_usage (); + endif + + if (isempty (k) || isempty (sigma) || isempty (mu) || ~ismatrix (k) || ~ismatrix (sigma) || ~ismatrix (mu)) + error ("gevstat: inputs must be numeric matrices"); + endif + + [retval, k, sigma, mu] = common_size (k, sigma, mu); + if (retval > 0) + error ("gevstat: inputs must be of common size or scalars"); + endif + + eg = 0.57721566490153286; %Euler-Mascheroni constant + + m = v = k; + + #find the mean + m(k >= 1) = Inf; + m(k == 0) = mu(k == 0) + eg*sigma(k == 0); + m(k < 1 & k ~= 0) = mu(k < 1 & k ~= 0) + sigma(k < 1 & k ~= 0) .* (gamma(1-k(k < 1 & k ~= 0)) - 1) ./ k(k < 1 & k ~= 0); + + #find the variance + v(k >= 0.5) = Inf; + v(k == 0) = (pi^2 / 6) * sigma(k == 0) .^ 2; + v(k < 0.5 & k ~= 0) = (gamma(1-2*k(k < 0.5 & k ~= 0)) - gamma(1-k(k < 0.5 & k ~= 0)).^2) .* (sigma(k < 0.5 & k ~= 0) ./ k(k < 0.5 & k ~= 0)) .^ 2; + + + +endfunction + +%!test +%! k = [-1 -0.5 0 0.2 0.4 0.5 1]; +%! sigma = 2; +%! mu = 1; +%! [m, v] = gevstat (k, sigma, mu); +%! expected_m = [1 1.4551 2.1544 2.6423 3.4460 4.0898 Inf]; +%! expected_v = [4 3.4336 6.5797 13.3761 59.3288 Inf Inf]; +%! assert (m, expected_m, -0.001); +%! assert (v, expected_v, -0.001);