# HG changeset patch
# User Dag Lyberg <daglyberg80@gmail.com>
# Date 1463176164 -7200
#      Fri May 13 23:49:24 2016 +0200
# Node ID f047e12820ae3ca21e4cc139aba16d31e1df8f78
# Parent  1449e3b98941b6400fafb0e0cfc11f411d07be28
Add statistical distributions to package 'statistics'

* bbspdf.m, bbscdf.m, bbsinv.m, bbsrnd.m: Birnbaum-Saunders distribution
* burrpdf.m, burrcdf.m, burrinv.m, burrrnd.m: Burr distribution
* gevpdf.m, gevcdf.m, gevinv.m, gevrnd.m: generalized extreme-value distribution
* gppdf.m, gpcdf.m, gpinv.m, gprnd.m: generalized Pareto distribution
* nakapdf.m, nakacdf.m, nakainv.m, nakarnd.m: Nakagami distribution
* raylpdf.m, raylcdf.m, raylinv.m, raylrnd.m: Rayleigh distribution
* tripdf.m, tricdf.m, triinv.m, trirnd.m: triangular distribution

diff -r 1449e3b98941 -r f047e12820ae scripts/statistics/distributions/bbscdf.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/statistics/distributions/bbscdf.m	Fri May 13 23:49:24 2016 +0200
@@ -0,0 +1,106 @@
+## Copyright (C) 2016 Dag Lyberg
+## Copyright (C) 1995-2015 Kurt Hornik
+##
+## This file is part of Octave.
+##
+## Octave 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 {} {} bbscdf (@var{x}, @var{location}, @var{scale}, @var{shape})
+## For each element of @var{x}, compute the cumulative distribution function
+## (CDF) at @var{x} of the Birnbaum-Saunders distribution with parameters
+## @var{location}, @var{scale} and @var{shape}.
+## @end deftypefn
+
+## Author: Dag Lyberg <daglyberg80@gmail.com>
+## Description: CDF of the Birnbaum-Saunders distribution
+
+function cdf = bbscdf (x, location, scale, shape)
+
+  if (nargin != 4)
+    print_usage ();
+  endif
+
+  if (! isscalar (location) || ! isscalar (scale) || ! isscalar(shape))
+    [retval, x, location, scale, shape] = ...
+        common_size (x, location, scale, shape);
+    if (retval > 0)
+      error ("bbscdf: X, LOCATION, SCALE and SHAPE must be of common size or scalars");
+    endif
+  endif
+
+  if (iscomplex (x) || iscomplex (location) || iscomplex (scale) ...
+      || iscomplex(shape))
+    error ("bbscdf: X, LOCATION, SCALE and SHAPE must not be complex");
+  endif
+
+  if (isa (x, "single") || isa (location, "single") || isa (scale, "single") ...
+      || isa (shape, "single"))
+    cdf = zeros (size (x), "single");
+  else
+    cdf = zeros (size (x));
+  endif
+
+  k = isnan(x) | ! (-Inf < location) | ! (location < Inf) ...
+      | ! (scale > 0) | ! (scale < Inf) | ! (shape > 0) | ! (shape < Inf);
+  cdf(k) = NaN;
+
+  k = (x > location) & (x <= Inf) & (-Inf < location) & (location < Inf) ...
+      & (0 < scale) & (scale < Inf) & (0 < shape) & (shape < Inf);
+  if (isscalar (location) && isscalar(scale) && isscalar(shape))
+    a = x(k) - location;
+    b = sqrt(a ./ scale);
+    cdf(k) = normcdf ((b - b.^-1) / shape);
+  else
+    a = x(k) - location(k);
+    b = sqrt(a ./ scale(k));
+    cdf(k) = normcdf ((b - b.^-1) ./ shape(k));
+  endif
+endfunction
+
+
+%!shared x,y
+%! x = [-1, 0, 1, 2, Inf];
+%! y = [0, 0, 1/2, 0.76024993890652337, 1];
+%!assert (bbscdf (x, zeros (1,5), ones (1,5), ones (1,5)), y, eps)
+%!assert (bbscdf (x, zeros (1,5), 1, 1), y, eps)
+%!assert (bbscdf (x, 0, ones (1,5), 1), y, eps)
+%!assert (bbscdf (x, 0, 1, ones (1,5)), y, eps)
+%!assert (bbscdf (x, 0, 1, 1), y, eps)
+%!assert (bbscdf (x, [0, 0, NaN, 0, 0], 1, 1), [y(1:2), NaN, y(4:5)], eps)
+%!assert (bbscdf (x, 0, [1, 1, NaN, 1, 1], 1), [y(1:2), NaN, y(4:5)], eps)
+%!assert (bbscdf (x, 0, 1, [1, 1, NaN, 1, 1]), [y(1:2), NaN, y(4:5)], eps)
+%!assert (bbscdf ([x, NaN], 0, 1, 1), [y, NaN], eps)
+
+## Test class of input preserved
+%!assert (bbscdf (single ([x, NaN]), 0, 1, 1), single ([y, NaN]), eps('single'))
+%!assert (bbscdf ([x, NaN], single (0), 1, 1), single ([y, NaN]), eps('single'))
+%!assert (bbscdf ([x, NaN], 0, single (1), 1), single ([y, NaN]), eps('single'))
+%!assert (bbscdf ([x, NaN], 0, 1, single (1)), single ([y, NaN]), eps('single'))
+
+## Test input validation
+%!error bbscdf ()
+%!error bbscdf (1)
+%!error bbscdf (1,2,3)
+%!error bbscdf (1,2,3,4,5)
+%!error bbscdf (ones (3), ones (2), ones(2), ones(2))
+%!error bbscdf (ones (2), ones (3), ones(2), ones(2))
+%!error bbscdf (ones (2), ones (2), ones(3), ones(2))
+%!error bbscdf (ones (2), ones (2), ones(2), ones(3))
+%!error bbscdf (i, 2, 3, 4)
+%!error bbscdf (1, i, 3, 4)
+%!error bbscdf (1, 2, i, 4)
+%!error bbscdf (1, 2, 3, i)
+
diff -r 1449e3b98941 -r f047e12820ae scripts/statistics/distributions/bbsinv.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/statistics/distributions/bbsinv.m	Fri May 13 23:49:24 2016 +0200
@@ -0,0 +1,114 @@
+## Copyright (C) 2016 Dag Lyberg
+## Copyright (C) 1995-2015 Kurt Hornik
+##
+## This file is part of Octave.
+##
+## Octave 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 {} {} bbsinv (@var{x}, @var{location}, @var{scale}, @var{shape})
+## For each element of @var{x}, compute the quantile (the inverse of the CDF)
+## at @var{x} of the Birnbaum-Saunders distribution with parameters
+## @var{location}, @var{scale}, and @var{shape}.
+## @end deftypefn
+
+## Author: Dag Lyberg <daglyberg80@gmail.com>
+## Description: Quantile function of the Birnbaum-Saunders distribution
+
+function inv = bbsinv (x, location, scale, shape)
+
+  if (nargin != 4)
+    print_usage ();
+  endif
+
+  if (! isscalar (location) || ! isscalar (scale) || ! isscalar(shape))
+    [retval, x, location, scale, shape] = ...
+        common_size (x, location, scale, shape);
+    if (retval > 0)
+      error ("bbsinv: X, LOCATION, SCALE and SHAPE must be of common size or scalars");
+    endif
+  endif
+
+  if (iscomplex (x) || iscomplex (location) ...
+      || iscomplex (scale) || iscomplex(shape))
+    error ("bbsinv: X, LOCATION, SCALE and SHAPE must not be complex");
+  endif
+
+  if (isa (x, "single") || isa (location, "single") ...
+      || isa (scale, "single") || isa (shape, "single"))
+    inv = zeros (size (x), "single");
+  else
+    inv = zeros (size (x));
+  endif
+
+  k = isnan (x) | (x < 0) | (x > 1) | ! (-Inf < location) | ! (location < Inf) ...
+      | ! (scale > 0) | ! (scale < Inf) | ! (shape > 0) | ! (shape < Inf);
+  inv(k) = NaN;
+
+  k = (x <= 0) & (-Inf < location) & (location < Inf) ...
+      & (scale > 0) & (scale < Inf) & (shape > 0) & (shape < Inf);
+  inv(k) = 0;
+
+  k = (x == 1) & (-Inf < location) & (location < Inf) ...
+      & (scale > 0) & (scale < Inf) & (shape > 0) & (shape < Inf);
+  inv(k) = Inf;
+
+  k = (0 < x) & (x < 1) & (location < Inf) & (0 < scale) & (scale < Inf) ...
+    & (0 < shape) & (shape < Inf);
+  if (isscalar (location) && isscalar(scale) && isscalar(shape))
+    y = shape * norminv (x(k));
+    inv(k) = location + scale * (y + sqrt (4 + y.^2)).^2 / 4;
+  else
+    y = shape(k) .* norminv (x(k));
+    inv(k) = location(k) + scale(k) .* (y + sqrt (4 + y.^2)).^2 ./ 4;
+  endif
+
+endfunction
+
+
+%!shared x,y
+#%! f = @(x,a,b,c) (a + b * (c * norminv (x) + sqrt (4 + (c * norminv(x))^2))^2) / 4;
+%! x = [-1, 0, 1/4, 1/2, 1, 2];
+%! y = [0, 0, 0.51565488603842446, 1, Inf, NaN];
+%!assert (bbsinv (x, zeros (1,6), ones (1,6), ones (1,6)), y, eps)
+%!assert (bbsinv (x, zeros (1,6), 1, 1), y, eps)
+%!assert (bbsinv (x, 0, ones (1,6), 1), y, eps)
+%!assert (bbsinv (x, 0, 1, ones (1,6)), y, eps)
+%!assert (bbsinv (x, 0, 1, 1), y, eps)
+%!assert (bbsinv (x, [0, 0, 0, NaN, 0, 0], 1, 1), [y(1:3), NaN, y(5:6)], eps)
+%!assert (bbsinv (x, 0, [1, 1, 1, NaN, 1, 1], 1), [y(1:3), NaN, y(5:6)], eps)
+%!assert (bbsinv (x, 0, 1, [1, 1, 1, NaN, 1, 1]), [y(1:3), NaN, y(5:6)], eps)
+%!assert (bbsinv ([x, NaN], 0, 1, 1), [y, NaN], eps)
+
+## Test class of input preserved
+%!assert (bbsinv (single ([x, NaN]), 0, 1, 1), single ([y, NaN]))
+%!assert (bbsinv ([x, NaN], single (0), 1, 1), single ([y, NaN]))
+%!assert (bbsinv ([x, NaN], 0, single (1), 1), single ([y, NaN]))
+%!assert (bbsinv ([x, NaN], 0, 1, single (1)), single ([y, NaN]))
+
+## Test input validation
+%!error bbsinv ()
+%!error bbsinv (1)
+%!error bbsinv (1,2,3)
+%!error bbsinv (1,2,3,4,5)
+%!error bbsinv (ones (3), ones (2), ones(2), ones(2))
+%!error bbsinv (ones (2), ones (3), ones(2), ones(2))
+%!error bbsinv (ones (2), ones (2), ones(3), ones(2))
+%!error bbsinv (ones (2), ones (2), ones(2), ones(3))
+%!error bbsinv (i, 2, 3, 4)
+%!error bbsinv (1, i, 3, 4)
+%!error bbsinv (1, 2, i, 4)
+%!error bbsinv (1, 2, 3, i)
+
diff -r 1449e3b98941 -r f047e12820ae scripts/statistics/distributions/bbspdf.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/statistics/distributions/bbspdf.m	Fri May 13 23:49:24 2016 +0200
@@ -0,0 +1,112 @@
+## Copyright (C) 2016 Dag Lyberg
+## Copyright (C) 1995-2015 Kurt Hornik
+##
+## This file is part of Octave.
+##
+## Octave 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 {} {} bbspdf (@var{x}, @var{location}, @var{scale}, @var{shape})
+## For each element of @var{x}, compute the probability density function (PDF)
+## at @var{x} of the Birnbaum-Saunders distribution with parameters
+## @var{location}, @var{scale} and @var{shape}.
+## @end deftypefn
+
+## Author: Dag Lyberg <daglyberg80@gmail.com>
+## Description: PDF of the Birnbaum-Saunders distribution
+
+function pdf = bbspdf (x, location, scale, shape)
+
+  if (nargin != 4)
+    print_usage ();
+  endif
+
+  if (! isscalar (location) || ! isscalar (scale) || ! isscalar(shape))
+    [retval, x, location, scale, shape] = ...
+        common_size (x, location, scale, shape);
+    if (retval > 0)
+      error ("bbspdf: X, LOCATION, SCALE and SHAPE must be of common size or scalars");
+    endif
+  endif
+
+  if (iscomplex (x) || iscomplex (location) ...
+      || iscomplex (scale) || iscomplex(shape))
+    error ("bbspdf: X, LOCATION, SCALE and SHAPE must not be complex");
+  endif
+
+  if (isa (x, "single") || isa (location, "single") || isa (scale, "single") ...
+      || isa (shape, "single"))
+    pdf = zeros (size (x), "single");
+  else
+    pdf = zeros (size (x));
+  endif
+
+  k = isnan (x) | ! (-Inf < location) | ! (location < Inf) ...
+      | ! (scale > 0) | ! (scale < Inf) ...
+      | ! (shape > 0) | ! (shape < Inf);
+
+  pdf(k) = NaN;
+
+  k = (x > location) & (x < Inf) & (-Inf < location) ...
+      & (location < Inf) & (0 < scale) & (scale < Inf) ...
+      & (0 < shape) & (shape < Inf);
+
+  if (isscalar (location) && isscalar(scale) && isscalar(shape))
+    a = x(k) - location;
+    b = sqrt(a ./ scale);
+    pdf(k) = ((b + b.^-1) ./ (2 * shape * a)) ...
+        .* normpdf ((b - b.^-1) / shape);
+  else
+    a = x(k) - location(k);
+    b = sqrt(a ./ scale(k));
+    pdf(k) = ((b + b.^-1) ./ (2 * shape(k).* a)) ...
+        .* normpdf ((b - b.^-1) ./ shape(k));
+  endif
+endfunction
+
+
+%!shared x,y
+%! x = [-1, 0, 1, 2, Inf];
+%! y = [0, 0, 0.3989422804014327, 0.1647717335503959, 0];
+%!assert (bbspdf (x, zeros (1,5), ones (1,5), ones (1,5)), y, eps)
+%!assert (bbspdf (x, zeros (1,5), 1, 1), y, eps)
+%!assert (bbspdf (x, 0, ones (1,5), 1), y, eps)
+%!assert (bbspdf (x, 0, 1, ones (1,5)), y, eps)
+%!assert (bbspdf (x, 0, 1, 1), y, eps)
+%!assert (bbspdf (x, [0, 0, NaN, 0, 0], 1, 1), [y(1:2), NaN, y(4:5)], eps)
+%!assert (bbspdf (x, 0, [1, 1, NaN, 1, 1], 1), [y(1:2), NaN, y(4:5)], eps)
+%!assert (bbspdf (x, 0, 1, [1, 1, NaN, 1, 1]), [y(1:2), NaN, y(4:5)], eps)
+%!assert (bbspdf ([x, NaN], 0, 1, 1), [y, NaN], eps)
+
+## Test class of input preserved
+%!assert (bbspdf (single ([x, NaN]), 0, 1, 1), single ([y, NaN]), eps('single'))
+%!assert (bbspdf ([x, NaN], single (0), 1, 1), single ([y, NaN]), eps('single'))
+%!assert (bbspdf ([x, NaN], 0, single (1), 1), single ([y, NaN]), eps('single'))
+%!assert (bbspdf ([x, NaN], 0, 1, single (1)), single ([y, NaN]), eps('single'))
+
+## Test input validation
+%!error bbspdf ()
+%!error bbspdf (1)
+%!error bbspdf (1,2,3)
+%!error bbspdf (1,2,3,4,5)
+%!error bbspdf (ones (3), ones (2), ones(2), ones(2))
+%!error bbspdf (ones (2), ones (3), ones(2), ones(2))
+%!error bbspdf (ones (2), ones (2), ones(3), ones(2))
+%!error bbspdf (ones (2), ones (2), ones(2), ones(3))
+%!error bbspdf (i, 2, 3, 4)
+%!error bbspdf (1, i, 3, 4)
+%!error bbspdf (1, 2, i, 4)
+%!error bbspdf (1, 2, 3, i)
+
diff -r 1449e3b98941 -r f047e12820ae scripts/statistics/distributions/bbsrnd.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/statistics/distributions/bbsrnd.m	Fri May 13 23:49:24 2016 +0200
@@ -0,0 +1,144 @@
+## Copyright (C) 2016 Dag Lyberg
+## Copyright (C) 1995-2015 Kurt Hornik
+##
+## This file is part of Octave.
+##
+## Octave 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  {} {} bbsrnd (@var{location}, @var{scale}, @var{shape})
+## @deftypefnx {} {} bbsrnd (@var{location}, @var{scale}, @var{shape}, @var{r})
+## @deftypefnx {} {} bbsrnd (@var{location}, @var{scale}, @var{shape}, @var{r}, @var{c}, @dots{})
+## @deftypefnx {} {} bbsrnd (@var{location}, @var{scale}, @var{shape}, [@var{sz}])
+## Return a matrix of random samples from the generalized Pareto distribution
+## with parameters @var{location}, @var{scale} and @var{shape}.
+##
+## When called with a single size argument, return 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 of dimensions @var{sz}.
+##
+## If no size arguments are given then the result matrix is the common size of
+## @var{location}, @var{scale} and @var{shape}.
+## @end deftypefn
+
+## Author: Dag Lyberg <daglyberg80@gmail.com>
+## Description: Random deviates from the Birnbaum-Saunders distribution
+
+function rnd = bbsrnd (location, scale, shape, varargin)
+
+  if (nargin < 3)
+    print_usage ();
+  endif
+
+  if (! isscalar (location) || ! isscalar (scale) || ! isscalar (shape))
+    [retval, location, scale, shape] = common_size (location, scale, shape);
+    if (retval > 0)
+      error ("bbsrnd: LOCATION, SCALE and SHAPE must be of common size or scalars");
+    endif
+  endif
+
+  if (iscomplex (location) || iscomplex (scale) || iscomplex (shape))
+    error ("bbsrnd: LOCATION, SCALE and SHAPE must not be complex");
+  endif
+
+  if (nargin == 3)
+    sz = size (location);
+  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 ("bbsrnd: dimension vector must be row vector of non-negative integers");
+    endif
+  elseif (nargin > 3)
+    if (any (cellfun (@(x) (! isscalar (x) || x < 0), varargin)))
+      error ("bbsrnd: dimensions must be non-negative integers");
+    endif
+    sz = [varargin{:}];
+  endif
+
+  if (! isscalar (location) && ! isequal (size (location), sz))
+    error ("bbsrnd: LOCATION, SCALE and SHAPE must be scalar or of size SZ");
+  endif
+
+  if (isa (location, "single") || isa (scale, "single") || isa (shape, "single"))
+    cls = "single";
+  else
+    cls = "double";
+  endif
+
+  if (isscalar (location) && isscalar (scale) && isscalar (shape))
+    if ((-Inf < location) && (location < Inf) ...
+        && (0 < scale) && (scale < Inf) ...
+        && (0 < shape) && (shape < Inf))
+      rnd = rand(sz,cls);
+      y = shape * norminv (rnd);
+      rnd = location + scale * (y + sqrt (4 + y.^2)).^2 / 4;
+    else
+      rnd = NaN (sz, cls);
+    endif
+  else
+    rnd = NaN (sz, cls);
+
+    k = (-Inf < location) & (location < Inf) ...
+        & (0 < scale) & (scale < Inf) ...
+        & (0 < shape) & (shape < Inf);
+    rnd(k) = rand(sum(k(:)),1);
+    y = shape(k) .* norminv (rnd(k));
+    rnd(k) = location(k) + scale(k) .* (y + sqrt (4 + y.^2)).^2 / 4;
+  endif
+endfunction
+
+
+%!assert (size (bbsrnd (0, 1, 1)), [1 1])
+%!assert (size (bbsrnd (zeros (2,1), 1, 1)), [2, 1])
+%!assert (size (bbsrnd (zeros (2,2), 1, 1)), [2, 2])
+%!assert (size (bbsrnd (0, ones (2,1), 1)), [2, 1])
+%!assert (size (bbsrnd (0, ones (2,2), 1)), [2, 2])
+%!assert (size (bbsrnd (0, 1, ones (2,1))), [2, 1])
+%!assert (size (bbsrnd (0, 1, ones (2,2))), [2, 2])
+%!assert (size (bbsrnd (0, 1, 1, 3)), [3, 3])
+%!assert (size (bbsrnd (0, 1, 1, [4 1])), [4, 1])
+%!assert (size (bbsrnd (0, 1, 1, 4, 1)), [4, 1])
+
+## Test class of input preserved
+%!assert (class (bbsrnd (0,1,1)), "double")
+%!assert (class (bbsrnd (single (0),1,1)), "single")
+%!assert (class (bbsrnd (single ([0 0]),1,1)), "single")
+%!assert (class (bbsrnd (0,single (1),1)), "single")
+%!assert (class (bbsrnd (0,single ([1 1]),1)), "single")
+%!assert (class (bbsrnd (0,1,single (1))), "single")
+%!assert (class (bbsrnd (0,1,single ([1 1]))), "single")
+
+## Test input validation
+%!error bbsrnd ()
+%!error bbsrnd (1)
+%!error bbsrnd (1,2)
+%!error bbsrnd (ones (3), ones (2), ones (2), 2)
+%!error bbsrnd (ones (2), ones (3), ones (2), 2)
+%!error bbsrnd (ones (2), ones (2), ones (3), 2)
+%!error bbsrnd (i, 2, 3)
+%!error bbsrnd (1, i, 3)
+%!error bbsrnd (1, 2, i)
+%!error bbsrnd (1,2,3, -1)
+%!error bbsrnd (1,2,3, ones (2))
+%!error bbsrnd (1,2,3, [2 -1 2])
+%!error bbsrnd (ones (2),1,2, 3)
+%!error bbsrnd (ones (2),1,2, [3, 2])
+%!error bbsrnd (ones (2),1,2, 3, 2)
+
diff -r 1449e3b98941 -r f047e12820ae scripts/statistics/distributions/burrcdf.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/statistics/distributions/burrcdf.m	Fri May 13 23:49:24 2016 +0200
@@ -0,0 +1,98 @@
+## Copyright (C) 2016 Dag Lyberg
+## Copyright (C) 1995-2015 Kurt Hornik
+##
+## This file is part of Octave.
+##
+## Octave 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 {} {} burrcdf (@var{x}, @var{c}, @var{k})
+## For each element of @var{x}, compute the cumulative distribution function
+## (CDF) at @var{x} of the Burr distribution with scale parameter @var{alpha}
+## and shape parameters @var{c} and @var{k}.
+## @end deftypefn
+
+## Author: Dag Lyberg <daglyberg80@gmail.com>
+## Description: CDF of the Burr distribution
+
+function cdf = burrcdf (x, alpha, c, k)
+
+  if (nargin != 4)
+    print_usage ();
+  endif
+
+  if (! isscalar (alpha) || ! isscalar (c) || ! isscalar (k) )
+    [retval, x, alpha, c, k] = common_size (x, alpha, c, k);
+    if (retval > 0)
+      error ("burrcdf: X, ALPHA, C AND K must be of common size or scalars");
+    endif
+  endif
+
+  if (iscomplex (x) || iscomplex(alpha) || iscomplex (c) || iscomplex (k))
+    error ("burrcdf: X, ALPHA, C AND K must not be complex");
+  endif
+
+  if (isa (x, "single") || isa (alpha, "single") || isa (c, "single") ...
+      || isa (k, "single"))
+    cdf = zeros (size (x), "single");
+  else
+    cdf = zeros (size (x));
+  endif
+
+  j = isnan (x) | ! (alpha > 0) | ! (c > 0) | ! (k > 0);
+  cdf(j) = NaN;
+
+  j = (x > 0) & (0 < alpha) & (alpha < Inf) & (0 < c) & (c < Inf) ...
+    & (0 < k) & (k < Inf);
+  if (isscalar (alpha) && isscalar(c) && isscalar(k))
+    cdf(j) = 1 - (1 + (x(j) / alpha).^c).^(-k);
+  else
+    cdf(j) = 1 - (1 + (x(j) ./ alpha(j)).^c(j)).^(-k(j));
+  endif
+
+endfunction
+
+
+%!shared x,y
+%! x = [-1, 0, 1, 2, Inf];
+%! y = [0, 0, 1/2, 2/3, 1];
+%!assert (burrcdf (x, ones(1,5), ones (1,5), ones (1,5)), y, eps)
+%!assert (burrcdf (x, 1, 1, 1), y, eps)
+%!assert (burrcdf (x, [1, 1, NaN, 1, 1], 1, 1), [y(1:2), NaN, y(4:5)], eps)
+%!assert (burrcdf (x, 1, [1, 1, NaN, 1, 1], 1), [y(1:2), NaN, y(4:5)], eps)
+%!assert (burrcdf (x, 1, 1, [1, 1, NaN, 1, 1]), [y(1:2), NaN, y(4:5)], eps)
+%!assert (burrcdf ([x, NaN], 1, 1, 1), [y, NaN], eps)
+
+## Test class of input preserved
+%!assert (burrcdf (single ([x, NaN]), 1, 1, 1), single ([y, NaN]), eps('single'))
+%!assert (burrcdf ([x, NaN], single (1), 1, 1), single ([y, NaN]), eps('single'))
+%!assert (burrcdf ([x, NaN], 1, single (1), 1), single ([y, NaN]), eps('single'))
+%!assert (burrcdf ([x, NaN], 1, 1, single (1)), single ([y, NaN]), eps('single'))
+
+## Test input validation
+%!error burrcdf ()
+%!error burrcdf (1)
+%!error burrcdf (1,2)
+%!error burrcdf (1,2,3)
+%!error burrcdf (1,2,3,4,5)
+%!error burrcdf (ones (3), ones (2), ones(2), ones(2))
+%!error burrcdf (ones (2), ones (3), ones(2), ones(2))
+%!error burrcdf (ones (2), ones (2), ones(3), ones(2))
+%!error burrcdf (ones (2), ones (2), ones(2), ones(3))
+%!error burrcdf (i, 2, 2, 2)
+%!error burrcdf (2, i, 2, 2)
+%!error burrcdf (2, 2, i, 2)
+%!error burrcdf (2, 2, 2, i)
+
diff -r 1449e3b98941 -r f047e12820ae scripts/statistics/distributions/burrinv.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/statistics/distributions/burrinv.m	Fri May 13 23:49:24 2016 +0200
@@ -0,0 +1,102 @@
+## Copyright (C) 2016 Dag Lyberg
+## Copyright (C) 1995-2015 Kurt Hornik
+##
+## This file is part of Octave.
+##
+## Octave 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 {} {} burrinv (@var{x}, @var{alpha}, @var{c}, @var{k})
+## For each element of @var{x}, compute the quantile (the inverse of the CDF)
+## at @var{x} of the Burr distribution with scale parameter @var{alpha} and
+## shape parameters @var{c} and @var{k}.
+## @end deftypefn
+
+## Author: Dag Lyberg <daglyberg80@gmail.com>
+## Description: Quantile function of the Burr distribution
+
+function inv = burrinv (x, alpha, c, k)
+
+  if (nargin != 4)
+    print_usage ();
+  endif
+
+  if (! isscalar (alpha) || ! isscalar (c) || ! isscalar (k) )
+    [retval, x, alpha, c, k] = common_size (x, alpha, c, k);
+    if (retval > 0)
+      error ("burrinv: X, ALPHA, C AND K must be of common size or scalars");
+    endif
+  endif
+
+  if (iscomplex (x) || iscomplex(alpha) || iscomplex (c) || iscomplex (k))
+    error ("burrinv: X, ALPHA, C AND K must not be complex");
+  endif
+
+  if (isa (x, "single") || isa (alpha, "single") || isa (c, "single") ...
+      || isa (k, "single"))
+    inv = zeros (size (x), "single");
+  else
+    inv = zeros (size (x));
+  endif
+
+  j = isnan (x) | (x < 0) | (x > 1) | ! (alpha > 0) | ! (c > 0) | ! (k > 0);
+  inv(j) = NaN;
+
+  j = (x == 1) & (0 < alpha) & (alpha < Inf) & (0 < c) & (c < Inf) ...
+      & (0 < k) & (k < Inf);
+  inv(j) = Inf;
+
+  j = (0 < x) & (x < 1) & (0 < alpha) & (alpha < Inf) & (0 < c) & (c < Inf) ...
+      & (0 < k) & (k < Inf);
+  if (isscalar (alpha) && isscalar(c) && isscalar(k))
+    inv(j) = ((1 - x(j) / alpha).^(-1 / k) - 1).^(1 / c) ;
+  else
+    inv(j) = ((1 - x(j) ./ alpha(j)).^(-1 ./ k(j)) - 1).^(1 ./ c(j)) ;
+  endif
+
+endfunction
+
+
+%!shared x,y
+%! x = [-Inf, -1, 0, 1/2, 1, 2, Inf];
+%! y = [NaN, NaN, 0, 1 , Inf, NaN, NaN];
+%!assert (burrinv (x, ones (1,7), ones (1,7), ones(1,7)), y, eps)
+%!assert (burrinv (x, 1, 1, 1), y, eps)
+%!assert (burrinv (x, [1, 1, 1, NaN, 1, 1, 1], 1, 1), [y(1:3), NaN, y(5:7)], eps)
+%!assert (burrinv (x, 1, [1, 1, 1, NaN, 1, 1, 1], 1), [y(1:3), NaN, y(5:7)], eps)
+%!assert (burrinv (x, 1, 1, [1, 1, 1, NaN, 1, 1, 1]), [y(1:3), NaN, y(5:7)], eps)
+%!assert (burrinv ([x, NaN], 1, 1, 1), [y, NaN], eps)
+
+## Test class of input preserved
+%!assert (burrinv (single ([x, NaN]), 1, 1, 1), single ([y, NaN]), eps('single'))
+%!assert (burrinv ([x, NaN], single (1), 1, 1), single ([y, NaN]), eps('single'))
+%!assert (burrinv ([x, NaN], 1, single (1), 1), single ([y, NaN]), eps('single'))
+%!assert (burrinv ([x, NaN], 1, 1, single (1)), single ([y, NaN]), eps('single'))
+
+## Test input validation
+%!error burrinv ()
+%!error burrinv (1)
+%!error burrinv (1,2)
+%!error burrinv (1,2,3)
+%!error burrinv (1,2,3,4,5)
+%!error burrinv (ones (3), ones (2), ones(2), ones(2))
+%!error burrinv (ones (2), ones (3), ones(2), ones(2))
+%!error burrinv (ones (2), ones (2), ones(3), ones(2))
+%!error burrinv (ones (2), ones (2), ones(2), ones(3))
+%!error burrinv (i, 2, 2, 2)
+%!error burrinv (2, i, 2, 2)
+%!error burrinv (2, 2, i, 2)
+%!error burrinv (2, 2, 2, i)
+
diff -r 1449e3b98941 -r f047e12820ae scripts/statistics/distributions/burrpdf.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/statistics/distributions/burrpdf.m	Fri May 13 23:49:24 2016 +0200
@@ -0,0 +1,100 @@
+## Copyright (C) 2016 Dag Lyberg
+## Copyright (C) 1995-2015 Kurt Hornik
+##
+## This file is part of Octave.
+##
+## Octave 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 {} {} burrpdf (@var{x}, @var{alpha}, @var{c}, @var{k})
+## For each element of @var{x}, compute the probability density function (PDF)
+## at @var{x} of the Burr distribution with scale parameter @var{alpha} and 
+## shape parameters @var{c} and @var{k}.
+## @end deftypefn
+
+## Author: Dag Lyberg <daglyberg80@gmail.com>
+## Description: PDF of the Burr distribution
+
+function pdf = burrpdf (x, alpha, c, k)
+
+  if (nargin != 4)
+    print_usage ();
+  endif
+
+  if (! isscalar (alpha) || ! isscalar (c) || ! isscalar (k) )
+    [retval, x, alpha, c, k] = common_size (x, alpha, c, k);
+    if (retval > 0)
+      error ("burrpdf: X, ALPHA, C AND K must be of common size or scalars");
+    endif
+  endif
+
+  if (iscomplex (x) || iscomplex(alpha) || iscomplex (c) || iscomplex (k))
+    error ("burrpdf: X, ALPHA, C AND K must not be complex");
+  endif
+
+  if (isa (x, "single") || isa (alpha, "single") ...
+      || isa (c, "single") || isa (k, "single"))
+    pdf = zeros (size (x), "single");
+  else
+    pdf = zeros (size (x));
+  endif
+
+  j = isnan (x) | ! (alpha > 0) | ! (c > 0) | ! (k > 0);
+  pdf(j) = NaN;
+
+  j = (x > 0) & (0 < alpha) & (alpha < Inf) & (0 < c) & (c < Inf) ...
+      & (0 < k) & (k < Inf);
+  if (isscalar (alpha) && isscalar (c) && isscalar(k))
+    pdf(j) = (c * k / alpha) .* (x(j) / alpha).^(c-1) ./ ...
+        (1 + (x(j) / alpha).^c).^(k + 1);
+  else
+    pdf(j) = (c(j) .* k(j) ./ alpha(j) ).* x(j).^(c(j)-1) ./ ...
+        (1 + (x(j) ./ alpha(j) ).^c(j) ).^(k(j) + 1);
+  endif
+
+endfunction
+
+
+%!shared x,y
+%! x = [-1, 0, 1, 2, Inf];
+%! y = [0, 0, 1/4, 1/9, 0];
+%!assert (burrpdf (x, ones(1,5), ones (1,5), ones (1,5)), y)
+%!assert (burrpdf (x, 1, 1, 1), y)
+%!assert (burrpdf (x, [1, 1, NaN, 1, 1], 1, 1), [y(1:2), NaN, y(4:5)])
+%!assert (burrpdf (x, 1, [1, 1, NaN, 1, 1], 1), [y(1:2), NaN, y(4:5)])
+%!assert (burrpdf (x, 1, 1, [1, 1, NaN, 1, 1]), [y(1:2), NaN, y(4:5)])
+%!assert (burrpdf ([x, NaN], 1, 1, 1), [y, NaN])
+
+## Test class of input preserved
+%!assert (burrpdf (single ([x, NaN]), 1, 1, 1), single ([y, NaN]))
+%!assert (burrpdf ([x, NaN], single (1), 1, 1), single ([y, NaN]))
+%!assert (burrpdf ([x, NaN], 1, single (1), 1), single ([y, NaN]))
+%!assert (burrpdf ([x, NaN], 1, 1, single (1)), single ([y, NaN]))
+
+## Test input validation
+%!error burrpdf ()
+%!error burrpdf (1)
+%!error burrpdf (1,2)
+%!error burrpdf (1,2,3)
+%!error burrpdf (1,2,3,4,5)
+%!error burrpdf (ones (3), ones (2), ones(2), ones(2))
+%!error burrpdf (ones (2), ones (3), ones(2), ones(2))
+%!error burrpdf (ones (2), ones (2), ones(3), ones(2))
+%!error burrpdf (ones (2), ones (2), ones(2), ones(3))
+%!error burrpdf (i, 2, 2, 2)
+%!error burrpdf (2, i, 2, 2)
+%!error burrpdf (2, 2, i, 2)
+%!error burrpdf (2, 2, 2, i)
+
diff -r 1449e3b98941 -r f047e12820ae scripts/statistics/distributions/burrrnd.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/statistics/distributions/burrrnd.m	Fri May 13 23:49:24 2016 +0200
@@ -0,0 +1,141 @@
+## Copyright (C) 2016 Dag Lyberg
+## Copyright (C) 1995-2015 Kurt Hornik
+##
+## This file is part of Octave.
+##
+## Octave 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  {} {} burrrnd (@var{alpha}, @var{c}, @var{k})
+## @deftypefnx {} {} burrrnd (@var{alpha}, @var{c}, @var{k}, @var{r})
+## @deftypefnx {} {} burrrnd (@var{alpha}, @var{c}, @var{k}, @var{r}, @var{c}, @dots{})
+## @deftypefnx {} {} burrrnd (@var{alpha}, @var{c}, @var{k}, [@var{sz}])
+## Return a matrix of random samples from the generalized Pareto distribution
+## with scale parameter @var{alpha} and shape parameters @var{c} and @var{k}.
+##
+## When called with a single size argument, return 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 of dimensions @var{sz}.
+##
+## If no size arguments are given then the result matrix is the common size of
+## @var{alpha}, @var{c} and @var{k}.
+## @end deftypefn
+
+## Author: Dag Lyberg <daglyberg80@gmail.com>
+## Description: Random deviates from the generalized extreme value (GEV) distribution
+
+function rnd = burrrnd (alpha, c, k, varargin)
+
+  if (nargin < 3)
+    print_usage ();
+  endif
+
+  if (! isscalar (alpha) || ! isscalar (c) || ! isscalar (k))
+    [retval, alpha, c, k] = common_size (alpha, c, k);
+    if (retval > 0)
+      error ("burrrnd: ALPHA, C and K must be of common size or scalars");
+    endif
+  endif
+
+  if (iscomplex (alpha) || iscomplex (c) || iscomplex (k))
+    error ("burrrnd: ALPHA, C and K must not be complex");
+  endif
+
+  if (nargin == 3)
+    sz = size (alpha);
+  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 ("burrrnd: dimension vector must be row vector of non-negative integers");
+    endif
+  elseif (nargin > 4)
+    if (any (cellfun (@(x) (! isscalar (x) || x < 0), varargin)))
+      error ("burrrnd: dimensions must be non-negative integers");
+    endif
+    sz = [varargin{:}];
+  endif
+
+  if (! isscalar (alpha) && ! isequal (size (c), sz) && ! isequal (size (k), sz))
+    error ("burrrnd: ALPHA, C and K must be scalar or of size SZ");
+  endif
+
+  if (isa (alpha, "single") || isa (c, "single") || isa (k, "single"))
+    cls = "single";
+  else
+    cls = "double";
+  endif
+
+  if (isscalar (alpha) && isscalar (c) && isscalar(k))
+    if ((0 < alpha) && (alpha < Inf) && (0 < c) && (c < Inf) ...
+        && (0 < k) && (k < Inf))
+      rnd = rand (sz, cls);
+      rnd(:) = ((1 - rnd(:) / alpha).^(-1 / k) - 1).^(1 / c);
+    else
+      rnd = NaN (sz, cls);
+    endif
+  else
+    rnd = NaN (sz, cls);
+
+    j = (0 < alpha) && (alpha < Inf) && (0 < c) && (c < Inf) ...
+        && (0 < k) && (k < Inf);
+    rnd(k) = rand(sum(j(:)),1);
+    rnd(k) = ((1 - rnd(j) / alpha(j)).^(-1 ./ k(j)) - 1).^(1 ./ c(j));
+  endif
+
+endfunction
+
+
+%!assert (size (burrrnd (1, 1, 1)), [1 1])
+%!assert (size (burrrnd (ones (2,1), 1, 1)), [2, 1])
+%!assert (size (burrrnd (ones (2,2), 1, 1)), [2, 2])
+%!assert (size (burrrnd (1, ones (2,1), 1)), [2, 1])
+%!assert (size (burrrnd (1, ones (2,2), 1)), [2, 2])
+%!assert (size (burrrnd (1, 1, ones (2,1))), [2, 1])
+%!assert (size (burrrnd (1, 1, ones (2,2))), [2, 2])
+%!assert (size (burrrnd (1, 1, 1, 3)), [3, 3])
+%!assert (size (burrrnd (1, 1, 1, [4 1])), [4, 1])
+%!assert (size (burrrnd (1, 1, 1, 4, 1)), [4, 1])
+
+## Test class of input preserved
+%!assert (class (burrrnd (1,1,1)), "double")
+%!assert (class (burrrnd (single (1),1,1)), "single")
+%!assert (class (burrrnd (single ([1 1]),1,1)), "single")
+%!assert (class (burrrnd (1,single (1),1)), "single")
+%!assert (class (burrrnd (1,single ([1 1]),1)), "single")
+%!assert (class (burrrnd (1,1,single (1))), "single")
+%!assert (class (burrrnd (1,1,single ([1 1]))), "single")
+
+## Test input validation
+%!error burrrnd ()
+%!error burrrnd (1)
+%!error burrrnd (1,2)
+%!error burrrnd (ones (3), ones (2), ones (2), 2)
+%!error burrrnd (ones (2), ones (3), ones (2), 2)
+%!error burrrnd (ones (2), ones (2), ones (3), 2)
+%!error burrrnd (i, 2, 2)
+%!error burrrnd (2, i, 2)
+%!error burrrnd (2, 2, i)
+%!error burrrnd (4,2,2, -1)
+%!error burrrnd (4,2,2, ones (2))
+%!error burrrnd (4,2,2, [2 -1 2])
+%!error burrrnd (4*ones (2),2,2, 3)
+%!error burrrnd (4*ones (2),2,2, [3, 2])
+%!error burrrnd (4*ones (2),2,2, 3, 2)
+
diff -r 1449e3b98941 -r f047e12820ae scripts/statistics/distributions/gevcdf.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/statistics/distributions/gevcdf.m	Fri May 13 23:49:24 2016 +0200
@@ -0,0 +1,169 @@
+## Copyright (C) 2016 Dag Lyberg
+## Copyright (C) 1995-2015 Kurt Hornik
+##
+## This file is part of Octave.
+##
+## Octave 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 {} {} gevcdf (@var{x}, @var{a}, @var{b}, @var{c})
+## Compute the cumulative distribution function (CDF) at @var{x} of the
+## generalized extreme value (GEV) distribution with parameters @var{location},
+## @var{scale}, and @var{shape}.
+## @end deftypefn
+
+## Author: Dag Lyberg <daglyberg80@gmail.com>
+## Description: CDF of the generalized extreme value (GEV) distribution.
+
+function cdf = gevcdf (x, location, scale, shape)
+
+  if (nargin != 4)
+    print_usage ();
+  endif
+
+  if (! isscalar (location) | ! isscalar (shape) | ! isscalar (scale))
+    [retval, x, location, shape, scale] = ...
+        common_size (x, location, shape, scale);
+    if (retval > 0)
+      error ("gevcdf: X, LOCATION, SCALE and SHAPE must be of common size or scalars");
+    endif
+  endif
+
+  if (iscomplex (x) || iscomplex (location) || iscomplex (shape) ...
+      || iscomplex (scale))
+    error ("gevcdf: X, LOCATION, SCALE and SHAPE must not be complex");
+  endif
+
+  if (isa (x, "single") || isa (location, "single") || isa (shape, "single") ...
+      || isa (scale, "single"))
+    cdf = zeros (size (x), "single");
+  else
+    cdf = zeros (size (x));
+  endif
+
+  k = isnan (x) | ! (-Inf < location) | ! (location < Inf) ...
+      | ! (scale > 0) | ! (scale < Inf) | ! (-Inf < shape) | ! (shape < Inf);
+  cdf(k) = NaN;
+
+  # Gumbel, type I extreme value distribution
+  k = (shape == 0) & (-Inf <= x) & (x <= Inf) ...
+      & (-Inf < location) & (location < Inf) & (scale > 0) & (scale < Inf);
+  if (sum(k) > 0)
+    if (isscalar (location) && isscalar (shape) && isscalar (scale))
+      z = (x - location) / scale;
+      cdf(k) = exp (-exp (-z(k)));
+    elseif (sum(k) > 0)
+      z = (x - location) ./ scale;
+      cdf(k) = exp (-exp (-z(k)));
+    endif
+  endif
+
+  # Frechet, type II extreme value distribution
+  # Reversed Weibull, type III extreme value distribution
+  k = (shape > 0 | shape < 0) & (-Inf <= x) & (x <= Inf) ...
+      & (-Inf < location) & (location < Inf) & (scale > 0);
+  if (any (k))
+    if (isscalar (location) && isscalar (scale) && isscalar (shape))
+      z = shape * (x - location) / scale;
+      j = k & (z > -1);
+      d = 1/shape;
+      cdf(j) = exp (-(1 + z(j)).^(-d));
+
+      # Reversed Weibull distribution
+      j = k & (shape < 0) & (z <= -1);
+      cdf(j) = 1;
+    else
+      z = shape .* (x - location) ./ scale;
+      j = k & (z > -1);
+      d = 1./shape;
+      cdf(j) =  exp (-(1 + z(j)).^(-d(j)));
+
+      # Reversed Weibull distribution
+      j = k & (shape < 0) & (z <= -1);
+      cdf(j) = 1;
+    endif
+  endif
+endfunction
+
+
+%!shared x,yg,yf,yw
+%! x = [-Inf, -1, 0, 1, Inf];
+%! yg = [0, 0.065988035845312543, 0.36787944117144233, 0.69220062755534639, 1];
+%! yf = [0, 0, 0.36787944117144233, 0.60653065971263342, 1];
+%! yw = [0, 0.1353352832366127, 0.36787944117144233, 1, 1];
+%!assert (gevcdf (x, zeros (1,5), ones (1,5), zeros (1,5)), yg)
+%!assert (gevcdf (x, zeros (1,5), 1, 0), yg)
+%!assert (gevcdf (x, 0, ones (1,5), 0), yg)
+%!assert (gevcdf (x, 0, 1, zeros (1,5)), yg)
+%!assert (gevcdf (x, 0, 1, 0), yg)
+%!assert (gevcdf (x, [0, 0, NaN, 0, 0], 1, 0), [yg(1:2), NaN, yg(4:5)])
+%!assert (gevcdf (x, 0, [1, 1, NaN, 1, 1], 0), [yg(1:2), NaN, yg(4:5)])
+%!assert (gevcdf (x, 0, 1, [0, 0, NaN, 0, 0]), [yg(1:2), NaN, yg(4:5)])
+%!assert (gevcdf ([x, NaN], 0, 1, 0), [yg, NaN])
+
+%!assert (gevcdf (x, zeros (1,5), ones (1,5), ones (1,5)), yf)
+%!assert (gevcdf (x, zeros (1,5), 1, 1), yf)
+%!assert (gevcdf (x, 0, ones (1,5), 1), yf)
+%!assert (gevcdf (x, 0, 1, ones (1,5)), yf)
+%!assert (gevcdf (x, 0, 1, 1), yf)
+%!assert (gevcdf (x, [0, 0, NaN, 0, 0], 1, 1), [yf(1:2), NaN, yf(4:5)])
+%!assert (gevcdf (x, 0, [1, 1, NaN, 1, 1], 1), [yf(1:2), NaN, yf(4:5)])
+%!assert (gevcdf (x, 0, 1, [1, 1, NaN, 1, 1]), [yf(1:2), NaN, yf(4:5)])
+%!assert (gevcdf ([x, NaN], 0, 1, 1), [yf, NaN])
+
+%!assert (gevcdf (x, zeros (1,5), ones (1,5), -ones (1,5)), yw)
+%!assert (gevcdf (x, zeros (1,5), 1, -1), yw)
+%!assert (gevcdf (x, 0, ones (1,5), -1), yw)
+%!assert (gevcdf (x, 0, 1, -ones (1,5)), yw)
+%!assert (gevcdf (x, 0, 1, -1), yw)
+%!assert (gevcdf (x, [0, 0, NaN, 0, 0], 1, -1), [yw(1:2), NaN, yw(4:5)])
+%!assert (gevcdf (x, 0, [1, 1, NaN, 1, 1], -1), [yw(1:2), NaN, yw(4:5)])
+%!assert (gevcdf (x, 0, 1, -[1, 1, NaN, 1, 1]), [yw(1:2), NaN, yw(4:5)])
+%!assert (gevcdf ([x, NaN], 0, 1, -1), [yw, NaN])
+
+## Test class of input preserved
+%!assert (gevcdf ([x, NaN], 0,1,0), [yg, NaN])
+%!assert (gevcdf (single ([x, NaN]), 0,1,0), single ([yg, NaN]), eps('single'))
+%!assert (gevcdf ([x, NaN], single (0),1,0), single ([yg, NaN]), eps('single'))
+%!assert (gevcdf ([x, NaN], 0,single (1),0), single ([yg, NaN]), eps('single'))
+%!assert (gevcdf ([x, NaN], 0,1,single (0)), single ([yg, NaN]), eps('single'))
+
+%!assert (gevcdf ([x, NaN], 0,1,1), [yf, NaN])
+%!assert (gevcdf (single ([x, NaN]), 0,1,1), single ([yf, NaN]), eps('single'))
+%!assert (gevcdf ([x, NaN], single (0),1,1), single ([yf, NaN]), eps('single'))
+%!assert (gevcdf ([x, NaN], 0,single (1),1), single ([yf, NaN]), eps('single'))
+%!assert (gevcdf ([x, NaN], 0,1,single (1)), single ([yf, NaN]), eps('single'))
+
+%!assert (gevcdf ([x, NaN], 0,1,-1), [yw, NaN])
+%!assert (gevcdf (single ([x, NaN]), 0,1,-1), single ([yw, NaN]), eps('single'))
+%!assert (gevcdf ([x, NaN], single (0),1,-1), single ([yw, NaN]), eps('single'))
+%!assert (gevcdf ([x, NaN], 0,single (1),-1), single ([yw, NaN]), eps('single'))
+%!assert (gevcdf ([x, NaN], 0,1,single (-1)), single ([yw, NaN]), eps('single'))
+
+## Test input validation
+%!error gevcdf ()
+%!error gevcdf (1)
+%!error gevcdf (1,2)
+%!error gevcdf (1,2,3)
+%!error gevcdf (1,2,3,4,5)
+%!error gevcdf (ones (3), ones (2), ones (2), ones (2))
+%!error gevcdf (ones (2), ones (3), ones (2), ones (2))
+%!error gevcdf (ones (2), ones (2), ones (3), ones (2))
+%!error gevcdf (ones (2), ones (2), ones (2), ones (3))
+%!error gevcdf (i, 2, 2, 2)
+%!error gevcdf (2, i, 2, 2)
+%!error gevcdf (2, 2, i, 2)
+%!error gevcdf (2, 2, 2, i)
+
diff -r 1449e3b98941 -r f047e12820ae scripts/statistics/distributions/gevinv.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/statistics/distributions/gevinv.m	Fri May 13 23:49:24 2016 +0200
@@ -0,0 +1,162 @@
+## Copyright (C) 2016 Dag Lyberg
+## Copyright (C) 1995-2015 Kurt Hornik
+##
+## This file is part of Octave.
+##
+## Octave 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 {} {} gevinv (@var{x}, @var{location}, @var{scale}, @var{shape})
+## For each element of @var{x}, compute the quantile (the inverse of the CDF)
+## at @var{x} of the generalized extreme value (GEV) distribution with
+## parameters @var{location}, @var{scale}, and @var{shape}.
+## @end deftypefn
+
+## Author: Dag Lyberg <daglyberg80@gmail.com>
+## Description: Quantile function of the generalized extreme value (GEV) distribution
+
+function inv = gevinv (x, location, scale, shape)
+
+  if (nargin != 4)
+    print_usage ();
+  endif
+
+  if (! isscalar (location) | ! isscalar (scale) | ! isscalar (shape))
+    [retval, x, location, scale, shape] = ...
+        common_size (x, location, scale, shape);
+    if (retval > 0)
+      error ("gevinv: X, LOCATION, SCALE and SHAPE must be of common size or scalars");
+    endif
+  endif
+
+  if (iscomplex (x) || iscomplex (location) || iscomplex (scale) || iscomplex (shape))
+    error ("gevinv: X, LOCATION, SCALE and SHAPE must not be complex");
+  endif
+
+  if (isa (x, "single") || isa (location, "single") ...
+      || isa (scale, "single") || isa (shape, "single"))
+    inv = nan (size (x), "single");
+  else
+    inv = nan (size (x));
+  endif
+
+  # Gumbel, type I extreme value distribution
+  k = (shape == 0) & (0 <= x) & (x <= 1) & ...
+    (-Inf < location) & (location < Inf) & (scale > 0);
+  if (any (k))
+    if (isscalar (location) && isscalar (scale) && isscalar (shape))
+      inv(k) = -scale * log (-log(x(k))) + location;
+    else
+      inv(k) = -scale(k) .* log (-log(x(k))) + location(k);
+    endif
+  endif
+  
+  # Frechet, type II extreme value distribution
+  # Reversed Weibull, type III extreme value distribution
+  k = (shape > 0 | shape < 0) & (0 <= x) & (x <= 1) & ...
+    (scale > 0) & (scale < Inf) & (-Inf < location) & (location < Inf);
+  if (any (k))
+    if (isscalar (location) && isscalar (scale) && isscalar (shape))
+      inv(k) = (-log (x(k))).^(-shape) - 1;
+      inv(k) = (scale / shape) .* inv(k) + location;
+    else
+      inv(k) = (-log (x(k))).^(-shape(k)) - 1;
+      inv(k) = (scale(k) ./ shape(k)) .* inv(k) + location(k);
+    endif
+
+    # Frechet distribution
+    j = k & (shape > 0) & (x == 0);
+    inv(j) = NaN;
+    j = k & (shape > 0) & (x == 1);
+    inv(j) = Inf;
+
+    # Reversed Weibull distribution
+    j = k & (shape < 0) & (x == 0);
+    inv(j) = -Inf;
+    j = k & (shape < 0) & (x == 1);
+    inv(j) = NaN;
+  endif
+endfunction
+
+
+%!shared x,yg,yf,yw
+%! x = [-Inf, -1, 0, 0.5, 1, 2, Inf];
+%! yg = [NaN, NaN, -Inf, 0.36651292058166435, Inf, NaN, NaN];
+%! yf = [NaN, NaN, NaN, 0.44269504088896344, Inf, NaN, NaN];
+%! yw = [NaN, NaN, -Inf, 0.30685281944005471, NaN, NaN, NaN];
+%!assert (gevinv (x, zeros (1,7), ones (1,7), zeros (1,7)), yg, eps)
+%!assert (gevinv (x, zeros (1,7), 1, 0), yg, eps)
+%!assert (gevinv (x, 0, ones (1,7), 0), yg, eps)
+%!assert (gevinv (x, 0, 1, zeros (1,7)), yg, eps)
+%!assert (gevinv (x, 0, 1, 0), yg, eps)
+%!assert (gevinv (x, [0, 0, 0, NaN, 0, 0, 0], 1, 0), [yg(1:3), NaN, yg(5:7)], eps)
+%!assert (gevinv (x, 0, [1, 1, 1, NaN, 1, 1, 1], 0), [yg(1:3), NaN, yg(5:7)], eps)
+%!assert (gevinv (x, 0, 1, [0, 0, 0, NaN, 0, 0, 0]), [yg(1:3), NaN, yg(5:7)], eps)
+%!assert (gevinv ([x, NaN], 0, 1, 0), [yg, NaN], eps)
+
+%!assert (gevinv (x, zeros (1,7), ones (1,7), ones (1,7)), yf, eps)
+%!assert (gevinv (x, zeros (1,7), 1, 1), yf, eps)
+%!assert (gevinv (x, 0, ones (1,7), 1), yf, eps)
+%!assert (gevinv (x, 0, 1, ones (1,7)), yf, eps)
+%!assert (gevinv (x, 0, 1, 1), yf, eps)
+%!assert (gevinv (x, [0, 0, 0, NaN, 0, 0, 0], 1, 1), [yf(1:3), NaN, yf(5:7)], eps)
+%!assert (gevinv (x, 0, [1, 1, 1, NaN, 1, 1, 1], 1), [yf(1:3), NaN, yf(5:7)], eps)
+%!assert (gevinv (x, 0, 1, [1, 1, 1, NaN, 1, 1, 1]), [yf(1:3), NaN, yf(5:7)], eps)
+%!assert (gevinv ([x, NaN], 0, 1, 1), [yf, NaN], eps)
+
+%!assert (gevinv (x, zeros (1,7), ones (1,7), -ones (1,7)), yw, eps)
+%!assert (gevinv (x, zeros (1,7), 1, -1), yw, eps)
+%!assert (gevinv (x, 0, ones (1,7), -1), yw, eps)
+%!assert (gevinv (x, 0, 1, -ones (1,7)), yw, eps)
+%!assert (gevinv (x, 0, 1, -1), yw, eps)
+%!assert (gevinv (x, [0, 0, 0, NaN, 0, 0, 0], 1, -1), [yw(1:3), NaN, yw(5:7)], eps)
+%!assert (gevinv (x, 0, [1, 1, 1, NaN, 1, 1, 1], -1), [yw(1:3), NaN, yw(5:7)], eps)
+%!assert (gevinv (x, 0, 1, -[1, 1, 1, NaN, 1, 1, 1]), [yw(1:3), NaN, yw(5:7)], eps)
+%!assert (gevinv ([x, NaN], 0, 1, -1), [yw, NaN], eps)
+
+## Test class of input preserved
+%!assert (gevinv ([x, NaN], 0,1,0), [yg, NaN], eps('double'))
+%!assert (gevinv (single ([x, NaN]), 0,1,0), single ([yg, NaN]), eps('single'))
+%!assert (gevinv ([x, NaN], single (0),1,0), single ([yg, NaN]), eps('single'))
+%!assert (gevinv ([x, NaN], 0,single (1),0), single ([yg, NaN]), eps('single'))
+%!assert (gevinv ([x, NaN], 0,1,single (0)), single ([yg, NaN]), eps('single'))
+
+%!assert (gevinv ([x, NaN], 0,1,1), [yf, NaN], eps('double'))
+%!assert (gevinv (single ([x, NaN]), 0,1,1), single ([yf, NaN]), eps('single'))
+%!assert (gevinv ([x, NaN], single (0),1,1), single ([yf, NaN]), eps('single'))
+%!assert (gevinv ([x, NaN], 0,single (1),1), single ([yf, NaN]), eps('single'))
+%!assert (gevinv ([x, NaN], 0,1,single (1)), single ([yf, NaN]), eps('single'))
+
+%!assert (gevinv ([x, NaN], 0,1,-1), [yw, NaN], eps('double'))
+%!assert (gevinv (single ([x, NaN]), 0,1,-1), single ([yw, NaN]), eps('single'))
+%!assert (gevinv ([x, NaN], single (0),1,-1), single ([yw, NaN]), eps('single'))
+%!assert (gevinv ([x, NaN], 0,single (1),-1), single ([yw, NaN]), eps('single'))
+%!assert (gevinv ([x, NaN], 0,1,single (-1)), single ([yw, NaN]), eps('single'))
+
+## Test input validation
+%!error gevinv ()
+%!error gevinv (1)
+%!error gevinv (1,2)
+%!error gevinv (1,2,3)
+%!error gevinv (1,2,3,4,5)
+%!error gevinv (ones (3), ones (2), ones (2), ones (2))
+%!error gevinv (ones (2), ones (3), ones (2), ones (2))
+%!error gevinv (ones (2), ones (2), ones (3), ones (2))
+%!error gevinv (ones (2), ones (2), ones (2), ones (3))
+%!error gevinv (i, 2, 2, 2)
+%!error gevinv (2, i, 2, 2)
+%!error gevinv (2, 2, i, 2)
+%!error gevinv (2, 2, 2, i)
+
diff -r 1449e3b98941 -r f047e12820ae scripts/statistics/distributions/gevpdf.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/statistics/distributions/gevpdf.m	Fri May 13 23:49:24 2016 +0200
@@ -0,0 +1,161 @@
+## Copyright (C) 2016 Dag Lyberg
+## Copyright (C) 1995-2015 Kurt Hornik
+##
+## This file is part of Octave.
+##
+## Octave 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 {} {} gevpdf (@var{x}, @var{location}, @var{scale}, @var{shape})
+## Compute the probability density function (PDF) at @var{x} of the
+## generalized extreme value (GEV) distribution with parameters @var{location},
+## @var{scale}, and @var{shape}.
+## @end deftypefn
+
+## Author: Dag Lyberg <daglyberg80@gmail.com>
+## Description: PDF of the generalized extreme value (GEV) distribution.
+
+function pdf = gevpdf (x, location, scale, shape)
+
+  if (nargin != 4)
+    print_usage ();
+  endif
+
+  if (! isscalar (location) | ! isscalar (shape) | ! isscalar (scale))
+    [retval, x, location, shape, scale] = ...
+        common_size (x, location, shape, scale);
+    if (retval > 0)
+      error ("gevpdf: X, LOCATION, SCALE and SHAPE must be of common size or scalars");
+    endif
+  endif
+
+  if (iscomplex (x) || iscomplex (location) || iscomplex (shape) ...
+      || iscomplex (scale))
+    error ("gevpdf: X, LOCATION, SCALE and SHAPE must not be complex");
+  endif
+
+  if (isa (x, "single") || isa (location, "single") || isa (shape, "single") ...
+      || isa (scale, "single"))
+    pdf = zeros (size (x), "single");
+  else
+    pdf = zeros (size (x));
+  endif
+
+  k = isnan (x) | ! (-Inf < location) | ! (location < Inf) ...
+      | ! (scale > 0) | ! (scale < Inf) | ! (-Inf < shape) | ! (shape < Inf);
+  pdf(k) = NaN;
+
+  # Gumbel, type I extreme value distribution
+  k = (shape == 0) & (-Inf < x) & (x < Inf) ...
+      & (-Inf < location) & (location < Inf) & (scale > 0) & (scale < Inf);
+  if (any (k))
+    if (isscalar (location) && isscalar (shape) && isscalar (scale))
+      z = (x - location) / scale;
+      pdf(k) = (1/scale) * exp (-z(k) - exp (-z(k)));
+    else
+      z = (x - location) ./ scale;
+      pdf(k) = (1./scale(k)) .* exp (-z(k) - exp (-z(k)));
+    endif
+  endif
+  
+  # Frechet, type II extreme value distribution
+  # Reversed Weibull, type III extreme value distribution
+  k = ((shape < 0) | (shape > 0)) & (-Inf < x) & (x < Inf) ...
+      & (-Inf < location) & (location < Inf) & (scale > 0) & (scale < Inf);
+  if (any (k))
+    if (isscalar (location) && isscalar (shape) && isscalar (scale))
+      z = shape * (x - location) / scale;
+      j = k & (z > -1);
+      d = 1 / shape;
+      pdf(j) =  (1/scale) .* (1+z(j)).^(-d-1) .* exp (-(1+z(j)).^(-d)) ;
+    else
+      z = shape .* (x - location) ./ scale;
+      j = k & (z > -1);
+      d = 1 ./ shape;
+      pdf(j) =  (1./scale(j)) .* (1+z(j)).^(-d(j)-1) .* exp (-(1+z(j)).^(-d(j)));
+    endif
+  endif
+  
+endfunction
+
+
+%!shared x,yg,yf,yw
+%! x = [-Inf, -1, 0, 1, Inf];
+%! yg = [0, 0.1793740787340172, 0.36787944117144233, 0.25464638004358248, 0];
+%! yf = [0, 0, 0.36787944117144233, 0.15163266492815833, 0];
+%! yw = [0, 0.1353352832366127, 0.36787944117144233, 0, 0];
+%!assert (gevpdf (x, zeros (1,5), ones (1,5), zeros (1,5)), yg)
+%!assert (gevpdf (x, zeros (1,5), 1, 0), yg)
+%!assert (gevpdf (x, 0, ones (1,5), 0), yg)
+%!assert (gevpdf (x, 0, 1, zeros (1,5)), yg)
+%!assert (gevpdf (x, 0, 1, 0), yg)
+%!assert (gevpdf (x, [0, 0, NaN, 0, 0], 1, 0), [yg(1:2), NaN, yg(4:5)])
+%!assert (gevpdf (x, 0, [1, 1, NaN, 1, 1], 0), [yg(1:2), NaN, yg(4:5)])
+%!assert (gevpdf (x, 0, 1, [0, 0, NaN, 0, 0]), [yg(1:2), NaN, yg(4:5)])
+%!assert (gevpdf ([x, NaN], 0, 1, 0), [yg, NaN])
+
+%!assert (gevpdf (x, zeros (1,5), ones (1,5), ones (1,5)), yf, eps)
+%!assert (gevpdf (x, zeros (1,5), 1, 1), yf, eps)
+%!assert (gevpdf (x, 0, ones (1,5), 1), yf, eps)
+%!assert (gevpdf (x, 0, 1, ones (1,5)), yf, eps)
+%!assert (gevpdf (x, 0, 1, 1), yf, eps)
+%!assert (gevpdf (x, [0, 0, NaN, 0, 0], 1, 1), [yf(1:2), NaN, yf(4:5)], eps)
+%!assert (gevpdf (x, 0, [1, 1, NaN, 1, 1], 1), [yf(1:2), NaN, yf(4:5)], eps)
+%!assert (gevpdf (x, 0, 1, [1, 1, NaN, 1, 1]), [yf(1:2), NaN, yf(4:5)], eps)
+%!assert (gevpdf ([x, NaN], 0, 1, 1), [yf, NaN], eps)
+
+%!assert (gevpdf (x, zeros (1,5), ones (1,5), -ones (1,5)), yw, eps)
+%!assert (gevpdf (x, zeros (1,5), 1, -1), yw, eps)
+%!assert (gevpdf (x, 0, ones (1,5), -1), yw, eps)
+%!assert (gevpdf (x, 0, 1, -ones (1,5)), yw, eps)
+%!assert (gevpdf (x, 0, 1, -1), yw, eps)
+%!assert (gevpdf (x, [0, 0, NaN, 0, 0], 1, -1), [yw(1:2), NaN, yw(4:5)], eps)
+%!assert (gevpdf (x, 0, [1, 1, NaN, 1, 1], -1), [yw(1:2), NaN, yw(4:5)], eps)
+%!assert (gevpdf (x, 0, 1, [-1, -1, NaN, -1, -1]), [yw(1:2), NaN, yw(4:5)], eps)
+%!assert (gevpdf ([x, NaN], 0, 1, -1), [yw, NaN], eps)
+
+## Test class of input preserved
+%!assert (gevpdf ([x, NaN], 0,1,0), [yg, NaN], eps)
+%!assert (gevpdf (single ([x, NaN]), 0,1,0), single ([yg, NaN]), eps('single'))
+%!assert (gevpdf ([x, NaN], single (0),1,0), single ([yg, NaN]), eps('single'))
+%!assert (gevpdf ([x, NaN], 0,single (1),0), single ([yg, NaN]), eps('single'))
+%!assert (gevpdf ([x, NaN], 0,1,single (0)), single ([yg, NaN]), eps('single'))
+
+%!assert (gevpdf ([x, NaN], 0,1,1), [yf, NaN], eps)
+%!assert (gevpdf (single ([x, NaN]), 0,1,1), single ([yf, NaN]), eps('single'))
+%!assert (gevpdf ([x, NaN], single (0),1,1), single ([yf, NaN]), eps('single'))
+%!assert (gevpdf ([x, NaN], 0,single (1),1), single ([yf, NaN]), eps('single'))
+%!assert (gevpdf ([x, NaN], 0,1,single (1)), single ([yf, NaN]), eps('single'))
+
+%!assert (gevpdf ([x, NaN], 0,1,-1), [yw, NaN], eps)
+%!assert (gevpdf (single ([x, NaN]), 0,1,-1), single ([yw, NaN]), eps('single'))
+%!assert (gevpdf ([x, NaN], single (0),1,-1), single ([yw, NaN]), eps('single'))
+%!assert (gevpdf ([x, NaN], 0,single (1),-1), single ([yw, NaN]), eps('single'))
+%!assert (gevpdf ([x, NaN], 0,1,single (-1)), single ([yw, NaN]), eps('single'))
+
+## Test input validation
+%!error gevpdf ()
+%!error gevpdf (1)
+%!error gevpdf (1,2)
+%!error gevpdf (1,2,3)
+%!error gevpdf (1,2,3,4,5)
+%!error gevpdf (ones (3), ones (2), ones (2), ones (2))
+%!error gevpdf (ones (2), ones (3), ones (2), ones (2))
+%!error gevpdf (ones (2), ones (2), ones (3), ones (2))
+%!error gevpdf (ones (2), ones (2), ones (2), ones (3))
+%!error gevpdf (i, 2, 2)
+%!error gevpdf (2, i, 2)
+%!error gevpdf (2, 2, i)
+
diff -r 1449e3b98941 -r f047e12820ae scripts/statistics/distributions/gevrnd.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/statistics/distributions/gevrnd.m	Fri May 13 23:49:24 2016 +0200
@@ -0,0 +1,211 @@
+## Copyright (C) 2016 Dag Lyberg
+## Copyright (C) 1995-2015 Kurt Hornik
+##
+## This file is part of Octave.
+##
+## Octave 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  {} {} gevrnd (@var{location}, @var{scale}, @var{shape})
+## @deftypefnx {} {} gevrnd (@var{location}, @var{scale}, @var{shape}, @var{r})
+## @deftypefnx {} {} gevrnd (@var{location}, @var{scale}, @var{shape}, @var{r}, @var{c}, @dots{})
+## @deftypefnx {} {} gevrnd (@var{location}, @var{scale}, @var{shape}, [@var{sz}])
+## Return a matrix of random samples from the generalized extreme value (GEV)
+## distribution.
+##
+## When called with a single size argument, return 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 of dimensions @var{sz}.
+##
+## If no size arguments are given then the result matrix is the common size of
+## @var{location}, @var{scale} and @var{shape}.
+## @end deftypefn
+
+## Author: Dag Lyberg <daglyberg80@gmail.com>
+## Description: Random deviates from the generalized extreme value (GEV) distribution
+
+function rnd = gevrnd (location, scale, shape, varargin)
+
+  if (nargin < 3)
+    print_usage ();
+  endif
+
+  if (! isscalar (location) || ! isscalar (scale) || ! isscalar (shape))
+    [retval, location, scale, shape] = common_size (location, scale, shape);
+    if (retval > 0)
+      error ("gevrnd: LOCATION, SCALE and SHAPE must be of common size or scalars");
+    endif
+  endif
+
+  if (iscomplex (location) || iscomplex (scale) || iscomplex (shape))
+    error ("gevrnd: LOCATION, SCALE and SHAPE must not be complex");
+  endif
+
+  if (nargin == 3)
+    sz = size (location);
+  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 (location) && ! isequal (size (location), sz))
+    error ("gevrnd: LOCATION, SCALE and SHAPE must be scalar or of size SZ");
+  endif
+
+  if (isa (location, "single") || isa (scale, "single") || isa (shape, "single"))
+    cls = "single";
+  else
+    cls = "double";
+  endif
+
+  
+  if (isscalar (location) && isscalar (scale) && isscalar (shape))
+    rnd = rand (sz, cls);
+    # Gumbel, type I extreme value distribution
+    if ((shape == 0) && ! isinf(shape) && ! isinf(location) ...
+        && (scale > 0) && (scale < Inf))
+      rnd(:) = -scale * log (-log(rnd(:))) + location;
+    # Frechet, type II extreme value distribution
+    elseif ((shape > 0) && ! isinf(shape) && ! isinf(location) ...
+        && (scale > 0) && (scale < Inf))
+      rnd(:) = (-log (rnd(:))).^(-shape) - 1;
+      rnd = (scale / shape) .* rnd + location;
+    # Reversed Weibull, type III extreme value distribution
+    elseif ((shape < 0) && ! isinf(shape) && ! isinf(location) ...
+        && (scale > 0) && (scale < Inf))
+      rnd(:) = (-log (rnd(:))).^(-shape) - 1;
+      rnd = (scale / shape) .* rnd + location;
+    else
+      rnd = NaN (sz, cls);
+    endif
+  else
+    rnd = rand (sz, cls);
+
+    # Gumbel, type I extreme value distribution
+    k = (shape == 0) & ! isinf(shape) & ! isinf(location) ...
+        & (scale > 0) & (scale < Inf);
+    if (any (k))
+      rnd = -scale .* log (-log(rnd)) + location;
+    endif
+
+    # Frechet, type II extreme value distribution
+    k = (shape > 0) & ! isinf(shape) & ! isinf(location) ...
+        & (scale > 0) & (scale < Inf);
+    if (any (k))
+      rnd = (scale./shape) .* ((-log (rnd)).^(-shape) - 1) + location;
+    endif
+
+    # Reversed Weibull, type III extreme value distribution
+    k = (shape < 0) & ! isinf(shape) & ! isinf(location) & ...
+        (scale > 0) & (scale < Inf);
+    if (any (k))
+      rnd = (scale./shape) .* ((-log (rnd)).^(-shape) - 1) + location;
+    endif
+
+    k = ! (-Inf < location) | ! (location < Inf) ...
+      | ! (scale > 0) | ! (scale < Inf) | ! (-Inf < shape) | ! (shape < Inf);
+    rnd(k) = NaN;
+  endif
+
+endfunction
+
+
+%!assert (size (gevrnd (0, 1, 0)), [1 1])
+%!assert (size (gevrnd (zeros (2,1), 1,0)), [2, 1])
+%!assert (size (gevrnd (zeros (2,2), 1,0)), [2, 2])
+%!assert (size (gevrnd (0, ones (2,1), 0)), [2, 1])
+%!assert (size (gevrnd (0, ones (2,2), 0)), [2, 2])
+%!assert (size (gevrnd (0, 1, zeros (2,1))), [2, 1])
+%!assert (size (gevrnd (0, 1, zeros (2,2))), [2, 2])
+%!assert (size (gevrnd (0, 1, 0, 3)), [3, 3])
+%!assert (size (gevrnd (0, 1, 0, [4 1])), [4, 1])
+%!assert (size (gevrnd (0, 1, 0, 4, 1)), [4, 1])
+
+%!assert (size (gevrnd (0, 1, 1)), [1 1])
+%!assert (size (gevrnd (zeros (2,1), 1,1)), [2, 1])
+%!assert (size (gevrnd (zeros (2,2), 1,1)), [2, 2])
+%!assert (size (gevrnd (0, ones (2,1), 1)), [2, 1])
+%!assert (size (gevrnd (0, ones (2,2), 1)), [2, 2])
+%!assert (size (gevrnd (0, 1, ones (2,1))), [2, 1])
+%!assert (size (gevrnd (0, 1, ones (2,2))), [2, 2])
+%!assert (size (gevrnd (0, 1, 1, 3)), [3, 3])
+%!assert (size (gevrnd (0, 1, 1, [4 1])), [4, 1])
+%!assert (size (gevrnd (0, 1, 1, 4, 1)), [4, 1])
+
+%!assert (size (gevrnd (0, 1, -1)), [1 1])
+%!assert (size (gevrnd (zeros (2,1), 1,-1)), [2, 1])
+%!assert (size (gevrnd (zeros (2,2), 1,-1)), [2, 2])
+%!assert (size (gevrnd (0, ones (2,1), -1)), [2, 1])
+%!assert (size (gevrnd (0, ones (2,2), -1)), [2, 2])
+%!assert (size (gevrnd (0, 1, -ones (2,1))), [2, 1])
+%!assert (size (gevrnd (0, 1, -ones (2,2))), [2, 2])
+%!assert (size (gevrnd (0, 1, -1, 3)), [3, 3])
+%!assert (size (gevrnd (0, 1, -1, [4 1])), [4, 1])
+%!assert (size (gevrnd (0, 1, -1, 4, 1)), [4, 1])
+
+## Test class of input preserved
+%!assert (class (gevrnd (0,1,0)), "double")
+%!assert (class (gevrnd (single (0),1,0)), "single")
+%!assert (class (gevrnd (single ([0 0]),1,0)), "single")
+%!assert (class (gevrnd (0,single (1),0)), "single")
+%!assert (class (gevrnd (0,single ([1 1]),0)), "single")
+%!assert (class (gevrnd (0,1,single (0))), "single")
+%!assert (class (gevrnd (0,1,single ([0 0]))), "single")
+
+%!assert (class (gevrnd (0,1,1)), "double")
+%!assert (class (gevrnd (single (0),1,1)), "single")
+%!assert (class (gevrnd (single ([0 0]),1,1)), "single")
+%!assert (class (gevrnd (0,single (1),1)), "single")
+%!assert (class (gevrnd (0,single ([1 1]),1)), "single")
+%!assert (class (gevrnd (0,1,single (1))), "single")
+%!assert (class (gevrnd (0,1,single ([1 1]))), "single")
+
+%!assert (class (gevrnd (0,1,-1)), "double")
+%!assert (class (gevrnd (single (0),1,-1)), "single")
+%!assert (class (gevrnd (single ([0 0]),1,-1)), "single")
+%!assert (class (gevrnd (0,single (1),-1)), "single")
+%!assert (class (gevrnd (0,single ([1 1]),-1)), "single")
+%!assert (class (gevrnd (0,1,single (-1))), "single")
+%!assert (class (gevrnd (0,1,single ([-1 -1]))), "single")
+
+## Test input validation
+%!error gevrnd ()
+%!error gevrnd (1)
+%!error gevrnd (1,2)
+%!error gevrnd (ones (3), ones (2), ones (2), 2)
+%!error gevrnd (ones (2), ones (3), ones (2), 2)
+%!error gevrnd (ones (2), ones (2), ones (3), 2)
+%!error gevrnd (i, 2, 3)
+%!error gevrnd (1, i, 3)
+%!error gevrnd (1, 2, i)
+%!error gevrnd (1,2,3, -1)
+%!error gevrnd (1,2,3, ones (2))
+%!error gevrnd (1,2,3, [2 -1 2])
+%!error gevrnd (ones (2),1,2, 3)
+%!error gevrnd (ones (2),1,2, [3, 2])
+%!error gevrnd (ones (2),1,2, 3, 2)
+
diff -r 1449e3b98941 -r f047e12820ae scripts/statistics/distributions/gpcdf.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/statistics/distributions/gpcdf.m	Fri May 13 23:49:24 2016 +0200
@@ -0,0 +1,175 @@
+## Copyright (C) 2016 Dag Lyberg
+## Copyright (C) 1997-2015 Kurt Hornik
+##
+## This file is part of Octave.
+##
+## Octave 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 {} {} gpcdf (@var{x}, @var{location}, @var{scale}, @var{shape})
+## Compute the cumulative distribution function (CDF) at @var{x} of the
+## generalized Pareto distribution with parameters @var{location}, @var{scale},
+## and @var{shape}.
+## @end deftypefn
+
+## Author: Dag Lyberg <daglyberg80@gmail.com>
+## Description: PDF of the generalized Pareto distribution
+
+function cdf = gpcdf (x, location, scale, shape)
+
+  if (nargin != 4)
+    print_usage ();
+  endif
+
+  if (! isscalar (location) || ! isscalar (scale) || ! isscalar (shape))
+    [retval, x, location, scale, shape] = ...
+      common_size (x, location, scale, shape);
+    if (retval > 0)
+      error ("gpcdf: X, LOCATION, SCALE and SHAPE must be of common size or scalars");
+    endif
+  endif
+
+  if (iscomplex (x) || iscomplex (location) || iscomplex (scale) ...
+      || iscomplex (shape))
+    error ("gpcdf: X, LOCATION, SCALE and SHAPE must not be complex");
+  endif
+
+  if (isa (x, "single") || isa (location, "single") || isa (scale, "single") ...
+      || isa (shape, "single"))
+    cdf = zeros (size (x), "single");
+  else
+    cdf = zeros (size (x));
+  endif
+
+  k = isnan (x) | ! (-Inf < location) | ! (location < Inf) | ! (scale > 0) ...
+    | ! (-Inf < shape) | ! (shape < Inf);
+  cdf(k) = NaN;
+
+  k = (x == Inf) & (-Inf < location) & (location < Inf) & (scale > 0) ...
+    & (-Inf < shape) & (shape < Inf);
+  cdf(k) = 1;
+
+  k = (-Inf < x) & (x < Inf) & (-Inf < location) & (location < Inf) ...
+      & (scale > 0) & (scale < Inf) & (-Inf < shape) & (shape < Inf);
+  if (isscalar (location) && isscalar (scale) && isscalar (shape))
+    z = (x - location) / scale;
+
+    j = k & (shape == 0) & (z >= 0);
+    if (any (j))
+      cdf(j) = 1 - exp (-z(j));
+    endif
+
+    j = k & (shape > 0) & (z >= 0);
+    if (any (j))
+      cdf(j) = 1 - (shape * z(j) + 1).^(-1 / shape);
+    endif
+
+    if (shape < 0)
+      j = k & (shape < 0) & (0 <= z) & (z <= -1 ./ shape);
+      if (any (j))
+        cdf(j) = 1 - (shape * z(j) + 1).^(-1 / shape);
+      endif
+    endif
+  else
+    z = (x - location) ./ scale;
+
+    j = k & (shape == 0) & (z >= 0);
+    if (any (j))
+      cdf(j) = 1 - exp (-z(j));
+    endif
+
+    j = k & (shape > 0) & (z >= 0);
+    if (any (j))
+      cdf(j) = 1 - (shape(j) .* z(j) + 1).^(-1 ./ shape(j));
+    endif
+    
+    if (any (shape < 0))
+      j = k & (shape < 0) & (0 <= z) & (z <= -1 ./ shape);
+      if (any (j))
+        cdf(j) = 1 - (shape(j) .* z(j) + 1).^(-1 ./ shape(j));
+      endif
+    endif
+  endif
+
+endfunction
+
+
+%!shared x,y1,y2,y3
+%! x = [-Inf, -1, 0, 1/2, 1, Inf];
+%! y1 = [0, 0, 0, 0.3934693402873666, 0.6321205588285577, 1];
+%! y2 = [0, 0, 0, 1/3, 1/2, 1];
+%! y3 = [0, 0, 0, 1/2, 1, 1];
+%! seps = eps('single')*5;
+%!assert (gpcdf (x, zeros (1,6), ones (1,6), zeros (1,6)), y1, eps)
+%!assert (gpcdf (x, zeros (1,6), 1, 0), y1, eps)
+%!assert (gpcdf (x, 0, ones (1,6), 0), y1, eps)
+%!assert (gpcdf (x, 0, 1, zeros (1,6)), y1, eps)
+%!assert (gpcdf (x, 0, 1, 0), y1, eps)
+%!assert (gpcdf (x, [0, 0, 0, NaN, 0, 0], 1, 0), [y1(1:3), NaN, y1(5:6)], eps)
+%!assert (gpcdf (x, 0, [1, 1, 1, NaN, 1, 1], 0), [y1(1:3), NaN, y1(5:6)], eps)
+%!assert (gpcdf (x, 0, 1, [0, 0, 0, NaN, 0, 0]), [y1(1:3), NaN, y1(5:6)], eps)
+%!assert (gpcdf ([x(1:3), NaN, x(5:6)], 0, 1, 0), [y1(1:3), NaN, y1(5:6)], eps)
+
+%!assert (gpcdf (x, zeros (1,6), ones (1,6), ones (1,6)), y2, eps)
+%!assert (gpcdf (x, zeros (1,6), 1, 1), y2, eps)
+%!assert (gpcdf (x, 0, ones (1,6), 1), y2, eps)
+%!assert (gpcdf (x, 0, 1, ones (1,6)), y2, eps)
+%!assert (gpcdf (x, 0, 1, 1), y2, eps)
+%!assert (gpcdf (x, [0, 0, 0, NaN, 0, 0], 1, 1), [y2(1:3), NaN, y2(5:6)], eps)
+%!assert (gpcdf (x, 0, [1, 1, 1, NaN, 1, 1], 1), [y2(1:3), NaN, y2(5:6)], eps)
+%!assert (gpcdf (x, 0, 1, [1, 1, 1, NaN, 1, 1]), [y2(1:3), NaN, y2(5:6)], eps)
+%!assert (gpcdf ([x(1:3), NaN, x(5:6)], 0, 1, 1), [y2(1:3), NaN, y2(5:6)], eps)
+
+%!assert (gpcdf (x, zeros (1,6), ones (1,6), -ones (1,6)), y3, eps)
+%!assert (gpcdf (x, zeros (1,6), 1, -1), y3, eps)
+%!assert (gpcdf (x, 0, ones (1,6), -1), y3, eps)
+%!assert (gpcdf (x, 0, 1, -ones (1,6)), y3, eps)
+%!assert (gpcdf (x, 0, 1, -1), y3, eps)
+%!assert (gpcdf (x, [0, 0, 0, NaN, 0, 0], 1, -1), [y1(1:3), NaN, y3(5:6)], eps)
+%!assert (gpcdf (x, 0, [1, 1, 1, NaN, 1, 1], -1), [y1(1:3), NaN, y3(5:6)], eps)
+%!assert (gpcdf (x, 0, 1, [-1, -1, -1, NaN, -1, -1]), [y1(1:3), NaN, y3(5:6)], eps)
+%!assert (gpcdf ([x(1:3), NaN, x(5:6)], 0, 1, -1), [y1(1:3), NaN, y3(5:6)], eps)
+
+## Test class of input preserved
+%!assert (gpcdf (single ([x, NaN]), 0, 1, 0), single ([y1, NaN]), eps('single'))
+%!assert (gpcdf ([x, NaN], single (0), 1, 0), single ([y1, NaN]), eps('single'))
+%!assert (gpcdf ([x, NaN], 0, single (1), 0), single ([y1, NaN]), eps('single'))
+%!assert (gpcdf ([x, NaN], 0, 1, single (0)), single ([y1, NaN]), eps('single'))
+
+%!assert (gpcdf (single ([x, NaN]), 0, 1, 1), single ([y2, NaN]), eps('single'))
+%!assert (gpcdf ([x, NaN], single (0), 1, 1), single ([y2, NaN]), eps('single'))
+%!assert (gpcdf ([x, NaN], 0, single (1), 1), single ([y2, NaN]), eps('single'))
+%!assert (gpcdf ([x, NaN], 0, 1, single (1)), single ([y2, NaN]), eps('single'))
+
+%!assert (gpcdf (single ([x, NaN]), 0, 1, -1), single ([y3, NaN]), eps('single'))
+%!assert (gpcdf ([x, NaN], single (0), 1, -1), single ([y3, NaN]), eps('single'))
+%!assert (gpcdf ([x, NaN], 0, single (1), -1), single ([y3, NaN]), eps('single'))
+%!assert (gpcdf ([x, NaN], 0, 1, single (-1)), single ([y3, NaN]), eps('single'))
+
+## Test input validation
+%!error gpcdf ()
+%!error gpcdf (1)
+%!error gpcdf (1,2)
+%!error gpcdf (1,2,3)
+%!error gpcdf (1,2,3,4,5)
+%!error gpcdf (ones (3), ones (2), ones (2), ones (2))
+%!error gpcdf (ones (2), ones (3), ones (2), ones (2))
+%!error gpcdf (ones (2), ones (2), ones (3), ones (2))
+%!error gpcdf (ones (2), ones (2), ones (2), ones (3))
+%!error gpcdf (i, 2, 2, 2)
+%!error gpcdf (2, i, 2, 2)
+%!error gpcdf (2, 2, i, 2)
+%!error gpcdf (2, 2, 2, i)
+
diff -r 1449e3b98941 -r f047e12820ae scripts/statistics/distributions/gpinv.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/statistics/distributions/gpinv.m	Fri May 13 23:49:24 2016 +0200
@@ -0,0 +1,161 @@
+## Copyright (C) 2016 Dag Lyberg
+## Copyright (C) 1997-2015 Kurt Hornik
+##
+## This file is part of Octave.
+##
+## Octave 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 {} {} gpinv (@var{x}, @var{location}, @var{scale}, @var{shape})
+## For each element of @var{x}, compute the quantile (the inverse of the CDF)
+## at @var{x} of the generalized Pareto distribution with parameters
+## @var{location}, @var{scale}, and @var{shape}.
+## @end deftypefn
+
+## Author: Dag Lyberg <daglyberg80@gmail.com>
+## Description: Quantile function of the generalized Pareto distribution
+
+function inv = gpinv (x, location, scale, shape)
+  if (nargin != 4)
+    print_usage ();
+  endif
+
+  if (! isscalar (location) || ! isscalar (scale) || ! isscalar (shape))
+    [retval, x, location, scale, shape] = ...
+        common_size (x, location, scale, shape);
+    if (retval > 0)
+      error ("gpinv: X, LOCATION, SCALE and SHAPE must be of common size or scalars");
+    endif
+  endif
+
+  if (iscomplex (x) || iscomplex (location) ...
+      || iscomplex (scale) || iscomplex (shape))
+    error ("gpinv: X, LOCATION, SCALE and SHAPE must not be complex");
+  endif
+
+  if (isa (x, "single") || isa (location, "single") ...
+      || isa (scale, "single") || isa (shape, "single"))
+    inv = zeros (size (x), "single");
+  else
+    inv = zeros (size (x));
+  endif
+
+  k = isnan (x) | ! (0 <= x) | ! (x <= 1) ...
+      | ! (-Inf < location) | ! (location < Inf) ...
+      | ! (scale > 0) | ! (scale < Inf) ...
+      | ! (-Inf < shape) | ! (shape < Inf);
+  inv(k) = NaN;
+
+  k = (0 <= x) & (x <= 1) & (-Inf < location) & (location < Inf) ...
+      & (scale > 0) & (scale < Inf) & (-Inf < shape) & (shape < Inf);
+  if (isscalar (location) && isscalar (scale) && isscalar (shape))
+    if (shape == 0)
+      inv(k) = -log(1 - x(k));
+      inv(k) = scale * inv(k) + location;
+    elseif (shape > 0)
+      inv(k) = (1 - x(k)).^(-shape) - 1;
+      inv(k) = (scale / shape) * inv(k) + location;
+    elseif (shape < 0)
+      inv(k) = (1 - x(k)).^(-shape) - 1;
+      inv(k) = (scale / shape) * inv(k)  + location;
+    end
+  else
+    j = k & (shape == 0);
+    if (any (j))
+      inv(j) = -log (1 - x(j));
+      inv(j) = scale(j) .* inv(j) + location(j);
+    endif
+    
+    j = k & (shape > 0);
+    if (any (j))
+      inv(j) = (1 - x(j)).^(-shape(j)) - 1;
+      inv(j) = (scale(j) ./ shape(j)) .* inv(j) + location(j);
+    endif
+    
+    j = k & (shape < 0);
+    if (any (j))
+      inv(j) = (1 - x(j)).^(-shape(j)) - 1;
+      inv(j) = (scale(j) ./ shape(j)) .* inv(j) + location(j);
+    endif
+  endif
+endfunction
+
+
+%!shared x,y1,y2,y3
+%! x = [-1, 0, 1/2, 1, 2];
+%! y1 = [NaN, 0, 0.6931471805599453, Inf, NaN];
+%! y2 = [NaN, 0, 1, Inf, NaN];
+%! y3 = [NaN, 0, 1/2, 1, NaN];
+%!assert (gpinv (x, zeros (1,5), ones (1,5), zeros (1,5)), y1)
+%!assert (gpinv (x, zeros (1,5), 1, 0), y1)
+%!assert (gpinv (x, 0, ones (1,5), 0), y1)
+%!assert (gpinv (x, 0, 1, zeros (1,5)), y1)
+%!assert (gpinv (x, 0, 1, 0), y1)
+%!assert (gpinv (x, [0, 0, NaN, 0, 0], 1, 0), [y1(1:2), NaN, y1(4:5)])
+%!assert (gpinv (x, 0, [1, 1, NaN, 1, 1], 0), [y1(1:2), NaN, y1(4:5)])
+%!assert (gpinv (x, 0, 1, [0, 0, NaN, 0, 0]), [y1(1:2), NaN, y1(4:5)])
+%!assert (gpinv ([x(1:2), NaN, x(4:5)], 0, 1, 0), [y1(1:2), NaN, y1(4:5)])
+
+%!assert (gpinv (x, zeros (1,5), ones (1,5), ones (1,5)), y2)
+%!assert (gpinv (x, zeros (1,5), 1, 1), y2)
+%!assert (gpinv (x, 0, ones (1,5), 1), y2)
+%!assert (gpinv (x, 0, 1, ones (1,5)), y2)
+%!assert (gpinv (x, 0, 1, 1), y2)
+%!assert (gpinv (x, [0, 0, NaN, 0, 0], 1, 1), [y2(1:2), NaN, y2(4:5)])
+%!assert (gpinv (x, 0, [1, 1, NaN, 1, 1], 1), [y2(1:2), NaN, y2(4:5)])
+%!assert (gpinv (x, 0, 1, [1, 1, NaN, 1, 1]), [y2(1:2), NaN, y2(4:5)])
+%!assert (gpinv ([x(1:2), NaN, x(4:5)], 0, 1, 1), [y2(1:2), NaN, y2(4:5)])
+
+%!assert (gpinv (x, zeros (1,5), ones (1,5), -ones (1,5)), y3)
+%!assert (gpinv (x, zeros (1,5), 1, -1), y3)
+%!assert (gpinv (x, 0, ones (1,5), -1), y3)
+%!assert (gpinv (x, 0, 1, -ones (1,5)), y3)
+%!assert (gpinv (x, 0, 1, -1), y3)
+%!assert (gpinv (x, [0, 0, NaN, 0, 0], 1, -1), [y3(1:2), NaN, y3(4:5)])
+%!assert (gpinv (x, 0, [1, 1, NaN, 1, 1], -1), [y3(1:2), NaN, y3(4:5)])
+%!assert (gpinv (x, 0, 1, -[1, 1, NaN, 1, 1]), [y3(1:2), NaN, y3(4:5)])
+%!assert (gpinv ([x(1:2), NaN, x(4:5)], 0, 1, -1), [y3(1:2), NaN, y3(4:5)])
+
+## Test class of input preserved
+%!assert (gpinv (single ([x, NaN]), 0, 1, 0), single ([y1, NaN]))
+%!assert (gpinv ([x, NaN], single (0), 1, 0), single ([y1, NaN]))
+%!assert (gpinv ([x, NaN], 0, single (1), 0), single ([y1, NaN]))
+%!assert (gpinv ([x, NaN], 0, 1, single (0)), single ([y1, NaN]))
+
+%!assert (gpinv (single ([x, NaN]), 0, 1, 1), single ([y2, NaN]))
+%!assert (gpinv ([x, NaN], single (0), 1, 1), single ([y2, NaN]))
+%!assert (gpinv ([x, NaN], 0, single (1), 1), single ([y2, NaN]))
+%!assert (gpinv ([x, NaN], 0, 1, single (1)), single ([y2, NaN]))
+
+%!assert (gpinv (single ([x, NaN]), 0, 1, -1), single ([y3, NaN]))
+%!assert (gpinv ([x, NaN], single (0), 1, -1), single ([y3, NaN]))
+%!assert (gpinv ([x, NaN], 0, single (1), -1), single ([y3, NaN]))
+%!assert (gpinv ([x, NaN], 0, 1, single (-1)), single ([y3, NaN]))
+
+## Test input validation
+%!error gpinv ()
+%!error gpinv (1)
+%!error gpinv (1,2)
+%!error gpinv (1,2,3)
+%!error gpinv (1,2,3,4,5)
+%!error gpinv (ones (3), ones (2), ones (2), ones (2))
+%!error gpinv (ones (2), ones (3), ones (2), ones (2))
+%!error gpinv (ones (2), ones (2), ones (3), ones (2))
+%!error gpinv (ones (2), ones (2), ones (2), ones (3))
+%!error gpinv (i, 2, 2, 2)
+%!error gpinv (2, i, 2, 2)
+%!error gpinv (2, 2, i, 2)
+%!error gpinv (2, 2, 2, i)
+
diff -r 1449e3b98941 -r f047e12820ae scripts/statistics/distributions/gppdf.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/statistics/distributions/gppdf.m	Fri May 13 23:49:24 2016 +0200
@@ -0,0 +1,169 @@
+## Copyright (C) 2016 Dag Lyberg
+## Copyright (C) 1997-2015 Kurt Hornik
+##
+## This file is part of Octave.
+##
+## Octave 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 {} {} gppdf (@var{x}, @var{location}, @var{scale}, @var{shape})
+## Compute the probability density function (PDF) at @var{x} of the
+## generalized Pareto distribution with parameters @var{location}, @var{scale},
+## and @var{shape}.
+## @end deftypefn
+
+## Author: Dag Lyberg <daglyberg80@gmail.com>
+## Description: PDF of the generalized Pareto distribution
+
+function pdf = gppdf (x, location, scale, shape)
+
+  if (nargin != 4)
+    print_usage ();
+  endif
+
+  if (! isscalar (location) || ! isscalar (scale) || ! isscalar (shape))
+    [retval, x, location, scale, shape] = ...
+      common_size (x, location, scale, shape);
+    if (retval > 0)
+      error ("gppdf: X, LOCATION, SCALE and SHAPE must be of common size or scalars");
+    endif
+  endif
+
+  if (iscomplex (x) || iscomplex (location) || iscomplex (scale) ...
+      || iscomplex (shape))
+    error ("gppdf: X, LOCATION, SCALE and SHAPE must not be complex");
+  endif
+
+  if (isa (x, "single") || isa (location, "single") || isa (scale, "single") ...
+      || isa (shape, "single"))
+    pdf = zeros (size (x), "single");
+  else
+    pdf = zeros (size (x));
+  endif
+
+  k = isnan (x) | ! (-Inf < location) | ! (location < Inf) | ...
+      ! (scale > 0) | ! (scale < Inf) | ! (-Inf < shape) | ! (shape < Inf);
+  pdf(k) = NaN;
+
+  k = (-Inf < x) & (x < Inf) & (-Inf < location) & (location < Inf) & ...
+        (scale > 0) & (scale < Inf) & (-Inf < shape) & (shape < Inf);
+  if (isscalar (location) && isscalar (scale) && isscalar (shape))
+    z = (x - location) / scale;
+
+    j = k & (shape == 0) & (z >= 0);
+    if (any (j))
+      pdf(j) = exp (-z(j));
+    endif
+
+    j = k & (shape > 0) & (z >= 0);
+    if (any (j))
+      pdf(j) = (shape * z(j) + 1).^(-(shape + 1) / shape);
+    endif
+
+    if (shape < 0)
+      j = k & (shape < 0) & (0 <= z) & (z <= -1. / shape);
+      if (any (j))
+        pdf(j) = (shape * z(j) + 1).^(-(shape + 1) / shape);
+      endif
+    endif
+  else
+    z = (x - location) ./ scale;
+
+    j = k & (shape == 0) & (z >= 0);
+    if (any (j))
+      pdf(j) = exp( -z(j));
+    endif
+
+    j = k & (shape > 0) & (z >= 0);
+    if (any (j))
+      pdf(j) = (shape(j) .* z(j) + 1).^(-(shape(j) + 1) ./ shape(j));
+    endif
+
+    if (any (shape < 0))
+      j = k & (shape < 0) & (0 <= z) & (z <= -1 ./ shape);
+      if (any (j))
+        pdf(j) = (shape(j) .* z(j) + 1).^(-(shape(j) + 1) ./ shape(j));
+      endif
+    endif
+  endif
+
+endfunction
+
+
+%!shared x,y1,y2,y3
+%! x = [-Inf, -1, 0, 1/2, 1, Inf];
+%! y1 = [0, 0, 1, 0.6065306597126334, 0.36787944117144233, 0];
+%! y2 = [0, 0, 1, 4/9, 1/4, 0];
+%! y3 = [0, 0, 1, 1, 1, 0];
+%!assert (gppdf (x, zeros (1,6), ones (1,6), zeros (1,6)), y1, eps)
+%!assert (gppdf (x, zeros (1,6), 1, 0), y1, eps)
+%!assert (gppdf (x, 0, ones (1,6), 0), y1, eps)
+%!assert (gppdf (x, 0, 1, zeros (1,6)), y1, eps)
+%!assert (gppdf (x, 0, 1, 0), y1, eps)
+%!assert (gppdf (x, [0, 0, 0, NaN, 0, 0], 1, 0), [y1(1:3), NaN, y1(5:6)])
+%!assert (gppdf (x, 0, [1, 1, 1, NaN, 1, 1], 0), [y1(1:3), NaN, y1(5:6)])
+%!assert (gppdf (x, 0, 1, [0, 0, 0, NaN, 0, 0]), [y1(1:3), NaN, y1(5:6)])
+%!assert (gppdf ([x(1:3), NaN, x(5:6)], 0, 1, 0), [y1(1:3), NaN, y1(5:6)])
+
+%!assert (gppdf (x, zeros (1,6), ones (1,6), ones (1,6)), y2, eps)
+%!assert (gppdf (x, zeros (1,6), 1, 1), y2, eps)
+%!assert (gppdf (x, 0, ones (1,6), 1), y2, eps)
+%!assert (gppdf (x, 0, 1, ones (1,6)), y2, eps)
+%!assert (gppdf (x, 0, 1, 1), y2, eps)
+%!assert (gppdf (x, [0, 0, 0, NaN, 0, 0], 1, 1), [y2(1:3), NaN, y2(5:6)])
+%!assert (gppdf (x, 0, [1, 1, 1, NaN, 1, 1], 1), [y2(1:3), NaN, y2(5:6)])
+%!assert (gppdf (x, 0, 1, [1, 1, 1, NaN, 1, 1]), [y2(1:3), NaN, y2(5:6)])
+%!assert (gppdf ([x(1:3), NaN, x(5:6)], 0, 1, 1), [y2(1:3), NaN, y2(5:6)])
+
+%!assert (gppdf (x, zeros (1,6), ones (1,6), -ones (1,6)), y3, eps)
+%!assert (gppdf (x, zeros (1,6), 1, -1), y3, eps)
+%!assert (gppdf (x, 0, ones (1,6), -1), y3, eps)
+%!assert (gppdf (x, 0, 1, -ones (1,6)), y3, eps)
+%!assert (gppdf (x, 0, 1, -1), y3, eps)
+%!assert (gppdf (x, [0, 0, 0, NaN, 0, 0], 1, -1), [y3(1:3), NaN, y3(5:6)])
+%!assert (gppdf (x, 0, [1, 1, 1, NaN, 1, 1], -1), [y3(1:3), NaN, y3(5:6)])
+%!assert (gppdf (x, 0, 1, [-1, -1, -1, NaN, -1, -1]), [y3(1:3), NaN, y3(5:6)])
+%!assert (gppdf ([x(1:3), NaN, x(5:6)], 0, 1, -1), [y3(1:3), NaN, y3(5:6)])
+
+## Test class of input preserved
+%!assert (gppdf (single ([x, NaN]), 0, 1, 0), single ([y1, NaN]))
+%!assert (gppdf ([x, NaN], single (0), 1, 0), single ([y1, NaN]))
+%!assert (gppdf ([x, NaN], 0, single (1), 0), single ([y1, NaN]))
+%!assert (gppdf ([x, NaN], 0, 1, single (0)), single ([y1, NaN]))
+
+%!assert (gppdf (single ([x, NaN]), 0, 1, 1), single ([y2, NaN]))
+%!assert (gppdf ([x, NaN], single (0), 1, 1), single ([y2, NaN]))
+%!assert (gppdf ([x, NaN], 0, single (1), 1), single ([y2, NaN]))
+%!assert (gppdf ([x, NaN], 0, 1, single (1)), single ([y2, NaN]))
+
+%!assert (gppdf (single ([x, NaN]), 0, 1, -1), single ([y3, NaN]))
+%!assert (gppdf ([x, NaN], single (0), 1, -1), single ([y3, NaN]))
+%!assert (gppdf ([x, NaN], 0, single (1), -1), single ([y3, NaN]))
+%!assert (gppdf ([x, NaN], 0, 1, single (-1)), single ([y3, NaN]))
+
+## Test input validation
+%!error gppdf ()
+%!error gppdf (1)
+%!error gppdf (1,2)
+%!error gppdf (1,2,3)
+%!error gppdf (1,2,3,4,5)
+%!error gppdf (1, ones (3), ones (2), ones (2))
+%!error gppdf (1, ones (2), ones (3), ones (2))
+%!error gppdf (1, ones (2), ones (2), ones (3))
+%!error gppdf (i, 2, 2, 2)
+%!error gppdf (2, i, 2, 2)
+%!error gppdf (2, 2, i, 2)
+%!error gppdf (2, 2, 2, i)
+
diff -r 1449e3b98941 -r f047e12820ae scripts/statistics/distributions/gprnd.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/statistics/distributions/gprnd.m	Fri May 13 23:49:24 2016 +0200
@@ -0,0 +1,174 @@
+## Copyright (C) 2016 Dag Lyberg
+## Copyright (C) 1995-2015 Kurt Hornik
+##
+## This file is part of Octave.
+##
+## Octave 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  {} {} gprnd (@var{location}, @var{scale}, @var{shape})
+## @deftypefnx {} {} gprnd (@var{location}, @var{scale}, @var{shape}, @var{r})
+## @deftypefnx {} {} gprnd (@var{location}, @var{scale}, @var{shape}, @var{r}, @var{c}, @dots{})
+## @deftypefnx {} {} gprnd (@var{location}, @var{scale}, @var{shape}, [@var{sz}])
+## Return a matrix of random samples from the generalized Pareto distribution
+## with parameters @var{location}, @var{scale} and @var{shape}.
+##
+## When called with a single size argument, return 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 of dimensions @var{sz}.
+##
+## If no size arguments are given then the result matrix is the common size of
+## @var{location}, @var{scale} and @var{shape}.
+## @end deftypefn
+
+## Author: Dag Lyberg <daglyberg80@gmail.com>
+## Description: Random deviates from the generalized Pareto distribution
+
+function rnd = gprnd (location, scale, shape, varargin)
+
+  if (nargin < 3)
+    print_usage ();
+  endif
+
+  if (! isscalar (location) || ! isscalar (scale) || ! isscalar (shape))
+    [retval, location, scale, shape] = common_size (location, scale, shape);
+    if (retval > 0)
+      error ("gpgrnd: LOCATION, SCALE and SHAPE must be of common size or scalars");
+    endif
+  endif
+
+  if (iscomplex (location) || iscomplex (scale) || iscomplex (shape))
+    error ("gprnd: LOCATION, SCALE and SHAPE must not be complex");
+  endif
+
+  if (nargin == 3)
+    sz = size (location);
+  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 ("gprnd: dimension vector must be row vector of non-negative integers");
+    endif
+  elseif (nargin > 4)
+    if (any (cellfun (@(x) (! isscalar (x) || x < 0), varargin)))
+      error ("gprnd: dimensions must be non-negative integers");
+    endif
+    sz = [varargin{:}];
+  endif
+
+  if (! isscalar (location) && ! isequal (size (location), sz))
+    error ("gprnd: LOCATION, SCALE and SHAPE must be scalar or of size SZ");
+  endif
+
+  if (isa (location, "single") || isa (scale, "single") || isa (shape, "single"))
+    cls = "single";
+  else
+    cls = "double";
+  endif
+
+  if (isscalar (location) && isscalar (scale) && isscalar (shape))
+    if ((-Inf < location) && (location < Inf) && (0 < scale) && (scale < Inf) ...
+          && (-Inf < shape) && (shape < Inf))
+      rnd = rand(sz,cls);
+      if (shape == 0)
+        rnd = -log(1 - rnd);
+        rnd = scale * rnd + location;
+      elseif ((shape < 0) || (shape > 0))
+        rnd = (1 - rnd).^(-shape) - 1;
+        rnd = (scale / shape) * rnd + location;
+      end
+    else
+      rnd = NaN (sz, cls);
+    endif
+  else
+    rnd = NaN (sz, cls);
+    
+    k = (-Inf < location) & (location < Inf) & (scale > 0) ...
+        & (-Inf < shape) & (shape < Inf);
+    rnd(k(:)) = rand (1, sum(k(:)), cls);
+    if (any (shape == 0))
+        rnd(k) = -log(1 - rnd(k));
+        rnd(k) = scale(k) .* rnd(k) + location(k);
+    elseif (any (shape < 0 | shape > 0))
+      rnd(k) = (1 - rnd(k)).^(-shape(k)) - 1;
+      rnd(k) = (scale(k) ./ shape(k)) .* rnd(k) + location(k);
+    end
+  endif
+endfunction
+
+
+%!assert (size (gprnd (0,1,0)), [1, 1])
+%!assert (size (gprnd (zeros (2,1), 1,0)), [2, 1])
+%!assert (size (gprnd (zeros (2,2), 1,0)), [2, 2])
+%!assert (size (gprnd (0, ones (2,1), 0)), [2, 1])
+%!assert (size (gprnd (0, ones (2,2), 0)), [2, 2])
+%!assert (size (gprnd (0,1, zeros (2,1))), [2, 1])
+%!assert (size (gprnd (0,1, zeros (2,2))), [2, 2])
+%!assert (size (gprnd (0,1, 0, 3)), [3, 3])
+%!assert (size (gprnd (0,1, 0, [4 1])), [4, 1])
+%!assert (size (gprnd (0,1, 0, 4, 1)), [4, 1])
+
+%!assert (size (gprnd (0,1,1)), [1, 1])
+%!assert (size (gprnd (zeros (2,1), 1,1)), [2, 1])
+%!assert (size (gprnd (zeros (2,2), 1,1)), [2, 2])
+%!assert (size (gprnd (0, ones (2,1), 1)), [2, 1])
+%!assert (size (gprnd (0, ones (2,2), 1)), [2, 2])
+%!assert (size (gprnd (0,1, ones (2,1))), [2, 1])
+%!assert (size (gprnd (0,1, ones (2,2))), [2, 2])
+%!assert (size (gprnd (0,1, 1, 3)), [3, 3])
+%!assert (size (gprnd (0,1, 1, [4 1])), [4, 1])
+%!assert (size (gprnd (0,1, 1, 4, 1)), [4, 1])
+
+%!assert (size (gprnd (0,1,-1)), [1, 1])
+%!assert (size (gprnd (zeros (2,1), -1,1)), [2, 1])
+%!assert (size (gprnd (zeros (2,2), -1,1)), [2, 2])
+%!assert (size (gprnd (0, ones (2,1), -1)), [2, 1])
+%!assert (size (gprnd (0, ones (2,2), -1)), [2, 2])
+%!assert (size (gprnd (0,1, -ones (2,1))), [2, 1])
+%!assert (size (gprnd (0,1, -ones (2,2))), [2, 2])
+%!assert (size (gprnd (0,1, -1, 3)), [3, 3])
+%!assert (size (gprnd (0,1, -1, [4 1])), [4, 1])
+%!assert (size (gprnd (0,1, -1, 4, 1)), [4, 1])
+
+## Test class of input preserved
+%!assert (class (gprnd (0,1,0)), "double")
+%!assert (class (gprnd (single (0),1,0)), "single")
+%!assert (class (gprnd (single ([0 0]),1,0)), "single")
+%!assert (class (gprnd (0,single (1),0)), "single")
+%!assert (class (gprnd (0,single ([1 1]),0)), "single")
+%!assert (class (gprnd (0,1,single (0))), "single")
+%!assert (class (gprnd (0,1,single ([0 0]))), "single")
+
+## Test input validation
+%!error gprnd ()
+%!error gprnd (1)
+%!error gprnd (1,2)
+%!error gprnd (zeros (3), ones (2), zeros (2))
+%!error gprnd (zeros (2), ones (3), zeros (2))
+%!error gprnd (zeros (2), ones (2), zeros (3))
+%!error gprnd (i, 1, 0)
+%!error gprnd (0, i, 0)
+%!error gprnd (0, 1, i)
+%!error gprnd (0,1,0, -1)
+%!error gprnd (0,1,0, ones (2))
+%!error gprnd (0,1,0, [2 -1 2])
+%!error gprnd (zeros (2),1,0, 3)
+%!error gprnd (zeros (2),1,0, [3, 2])
+%!error gprnd (zeros (2),1,0, 3, 2)
+
diff -r 1449e3b98941 -r f047e12820ae scripts/statistics/distributions/nakacdf.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/statistics/distributions/nakacdf.m	Fri May 13 23:49:24 2016 +0200
@@ -0,0 +1,103 @@
+## Copyright (C) 2016 Dag Lyberg
+## Copyright (C) 1995-2015 Kurt Hornik
+##
+## This file is part of Octave.
+##
+## Octave 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 {} {} nakacdf (@var{x}, @var{m}, @var{w})
+## For each element of @var{x}, compute the cumulative distribution function
+## (CDF) at @var{x} of the Nakagami distribution with shape parameter @var{m}
+## and scale parameter @var{w}.
+##
+## @end deftypefn
+
+## Author: Dag Lyberg <daglyberg80@gmail.com>
+## Description: CDF of the Nakagami distribution
+
+function cdf = nakacdf (x, m, w)
+
+  if (nargin != 3)
+    print_usage ();
+  endif
+
+  if (! isscalar (m) || ! isscalar (w))
+    [retval, x, m, w] = common_size (x, m, w);
+    if (retval > 0)
+      error ("nakacdf: X, M and W must be of common size or scalars");
+    endif
+  endif
+
+  if (iscomplex (x) || iscomplex (m) || iscomplex (w))
+    error ("nakacdf: X, M and W must not be complex");
+  endif
+
+  if (isa (x, "single") || isa (m, "single") || isa (w, "single"))
+    inv = zeros (size (x), "single");
+  else
+    inv = zeros (size (x));
+  endif
+
+  k = isnan (x) | ! (m > 0) | ! (w > 0);
+  cdf(k) = NaN;
+
+  k = (x == Inf) & (0 < m) & (m < Inf) & (0 < w) & (w < Inf);
+  cdf(k) = 1;
+  
+  k = (0 < x) & (x < Inf) & (0 < m) & (m < Inf) & (0 < w) & (w < Inf);
+  if (isscalar(x) && isscalar (m) && isscalar(w))
+    left = m;
+    right = (m/w) * x^2;
+    cdf(k) = gammainc(right, left);
+  elseif (isscalar (m) && isscalar(w))
+    left = m * ones(size(x));
+    right = (m/w) * x.^2;
+    cdf(k) = gammainc(right(k), left(k));
+  else
+    left = m .* ones(size(x));
+    right = (m./w) .* x.^2;
+    cdf(k) = gammainc(right(k), left(k));
+  endif
+
+endfunction
+
+
+%!shared x,y
+%! x = [-1, 0, 1, 2, Inf];
+%! y = [0, 0, 0.63212055882855778, 0.98168436111126578, 1];
+%!assert (nakacdf (x, ones (1,5), ones (1,5)), y)
+%!assert (nakacdf (x, 1, 1), y)
+%!assert (nakacdf (x, [1, 1, NaN, 1, 1], 1), [y(1:2), NaN, y(4:5)])
+%!assert (nakacdf (x, 1, [1, 1, NaN, 1, 1]), [y(1:2), NaN, y(4:5)])
+%!assert (nakacdf ([x, NaN], 1, 1), [y, NaN])
+
+## Test class of input preserved
+%!assert (nakacdf (single ([x, NaN]), 1, 1), single ([y, NaN]), eps('single'))
+%!assert (nakacdf ([x, NaN], single (1), 1), single ([y, NaN]), eps('single'))
+%!assert (nakacdf ([x, NaN], 1, single (1)), single ([y, NaN]), eps('single'))
+
+## Test input validation
+%!error nakacdf ()
+%!error nakacdf (1)
+%!error nakacdf (1,2)
+%!error nakacdf (1,2,3,4)
+%!error nakacdf (ones (3), ones (2), ones(2))
+%!error nakacdf (ones (2), ones (3), ones(2))
+%!error nakacdf (ones (2), ones (2), ones(3))
+%!error nakacdf (i, 2, 2)
+%!error nakacdf (2, i, 2)
+%!error nakacdf (2, 2, i)
+
diff -r 1449e3b98941 -r f047e12820ae scripts/statistics/distributions/nakainv.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/statistics/distributions/nakainv.m	Fri May 13 23:49:24 2016 +0200
@@ -0,0 +1,102 @@
+## Copyright (C) 2016 Dag Lyberg
+## Copyright (C) 1995-2015 Kurt Hornik
+##
+## This file is part of Octave.
+##
+## Octave 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 {} {} nakainv (@var{x}, @var{m}, @var{w})
+## For each element of @var{x}, compute the quantile (the inverse of the CDF)
+## at @var{x} of the Nakagami distribution with shape parameter @var{m} and
+## scale parameter @var{w}.
+##
+## @end deftypefn
+
+## Author: Dag Lyberg <daglyberg80@gmail.com>
+## Description: Quantile function of the Nakagami distribution
+
+function inv = nakainv (x, m, w)
+
+  if (nargin != 3)
+    print_usage ();
+  endif
+
+  if (! isscalar (m) || ! isscalar (w))
+    [retval, x, m, w] = common_size (x, m, w);
+    if (retval > 0)
+      error ("nakainv: X, M and W must be of common size or scalars");
+    endif
+  endif
+
+  if (iscomplex (x) || iscomplex (m) || iscomplex (w))
+    error ("nakainv: X, M, and W must not be complex");
+  endif
+
+  if (isa (x, "single") || isa (m, "single") || isa (w, "single"))
+    inv = zeros (size (x), "single");
+  else
+    inv = zeros (size (x));
+  endif
+
+  k = isnan (x) | ! (0 <= x) | ! (x <= 1) | ! (-Inf < m) | ! (m < Inf) ...
+    | ! (0 < w) | ! (w < Inf);
+  inv(k) = NaN;
+
+  k = (x == 1) & (-Inf < m) & (m < Inf) & (0 < w) & (w < Inf);
+  inv(k) = Inf;
+
+  k = (0 < x) & (x < 1) & (0 < m) & (m < Inf) & (0 < w) & (w < Inf);
+  if (isscalar (m) && isscalar(w))
+    m_gamma = m;
+    w_gamma = w/m;
+    inv(k) = gaminv(x(k), m_gamma, w_gamma);
+    inv(k) = sqrt(inv(k));
+  else
+    m_gamma = m;
+    w_gamma = w./m;
+    inv(k) = gaminv(x(k), m_gamma(k), w_gamma(k));
+    inv(k) = sqrt(inv(k));
+  endif
+
+endfunction
+
+
+%!shared x,y
+%! x = [-Inf, -1, 0, 1/2, 1, 2, Inf];
+%! y = [NaN, NaN, 0, 0.83255461115769769, Inf, NaN, NaN];
+%!assert (nakainv (x, ones (1,7), ones (1,7)), y, eps)
+%!assert (nakainv (x, 1, 1), y, eps)
+%!assert (nakainv (x, [1, 1, 1, NaN, 1, 1, 1], 1), [y(1:3), NaN, y(5:7)], eps)
+%!assert (nakainv (x, 1, [1, 1, 1, NaN, 1, 1, 1]), [y(1:3), NaN, y(5:7)], eps)
+%!assert (nakainv ([x, NaN], 1, 1), [y, NaN], eps)
+
+## Test class of input preserved
+%!assert (nakainv (single ([x, NaN]), 1, 1), single ([y, NaN]))
+%!assert (nakainv ([x, NaN], single (1), 1), single ([y, NaN]))
+%!assert (nakainv ([x, NaN], 1, single (1)), single ([y, NaN]))
+
+## Test input validation
+%!error nakainv ()
+%!error nakainv (1)
+%!error nakainv (1,2)
+%!error nakainv (1,2,3,4)
+%!error nakainv (ones (3), ones (2), ones(2))
+%!error nakainv (ones (2), ones (3), ones(2))
+%!error nakainv (ones (2), ones (2), ones(3))
+%!error nakainv (i, 2, 2)
+%!error nakainv (2, i, 2)
+%!error nakainv (2, 2, i)
+
diff -r 1449e3b98941 -r f047e12820ae scripts/statistics/distributions/nakapdf.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/statistics/distributions/nakapdf.m	Fri May 13 23:49:24 2016 +0200
@@ -0,0 +1,95 @@
+## Copyright (C) 2016 Dag Lyberg
+## Copyright (C) 1995-2015 Kurt Hornik
+##
+## This file is part of Octave.
+##
+## Octave 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 {} {} nakapdf (@var{x}, @var{m}, @var{w})
+## For each element of @var{x}, compute the probability density function (PDF)
+## at @var{x} of the Nakagami distribution with shape parameter @var{m} and
+## scale parameter @var{w}.
+## @end deftypefn
+
+## Author: Dag Lyberg <daglyberg80@gmail.com>
+## Description: PDF of the Nakagami distribution
+
+function pdf = nakapdf (x, m, w)
+
+  if (nargin != 3)
+    print_usage ();
+  endif
+
+  if (! isscalar (m) || ! isscalar (w))
+    [retval, x, m, w] = common_size (x, m, w);
+    if (retval > 0)
+      error ("nakapdf: X, M and W must be of common size or scalars");
+    endif
+  endif
+
+  if (iscomplex (x) || iscomplex (m) || iscomplex (w))
+    error ("nakapdf: X, M and W must not be complex");
+  endif
+
+  if (isa (x, "single") || isa (m, "single") || isa (w, "single"))
+    pdf = zeros (size (x), "single");
+  else
+    pdf = zeros (size (x));
+  endif
+
+  k = isnan (x) | ! (m > 0.5) | ! (w > 0);
+  pdf(k) = NaN;
+
+  k = (0 < x) & (x < Inf) & (0 < m) & (m < Inf) & (0 < w) & (w < Inf);
+  if (isscalar (m) && isscalar(w))
+    pdf(k) = exp (log (2) + m*log (m) - log (gamma (m)) - ...
+               m*log (w) + (2*m-1) * ...
+               log (x(k)) - (m/w) * x(k).^2);
+  else
+    pdf(k) = exp(log(2) + m(k).*log (m(k)) - log (gamma (m(k))) - ...
+               m(k).*log (w(k)) + (2*m(k)-1) ...
+               .* log (x(k)) - (m(k)./w(k)) .* x(k).^2);
+  endif
+
+endfunction
+
+
+%!shared x,y
+%! x = [-1, 0, 1, 2, Inf];
+%! y = [0, 0, 0.73575888234288467, 0.073262555554936715, 0];
+%!assert (nakapdf (x, ones (1,5), ones (1,5)), y, eps)
+%!assert (nakapdf (x, 1, 1), y, eps)
+%!assert (nakapdf (x, [1, 1, NaN, 1, 1], 1), [y(1:2), NaN, y(4:5)], eps)
+%!assert (nakapdf (x, 1, [1, 1, NaN, 1, 1]), [y(1:2), NaN, y(4:5)], eps)
+%!assert (nakapdf ([x, NaN], 1, 1), [y, NaN], eps)
+
+## Test class of input preserved
+%!assert (nakapdf (single ([x, NaN]), 1, 1), single ([y, NaN]))
+%!assert (nakapdf ([x, NaN], single (1), 1), single ([y, NaN]))
+%!assert (nakapdf ([x, NaN], 1, single (1)), single ([y, NaN]))
+
+## Test input validation
+%!error nakapdf ()
+%!error nakapdf (1)
+%!error nakapdf (1,2)
+%!error nakapdf (1,2,3,4)
+%!error nakapdf (ones (3), ones (2), ones(2))
+%!error nakapdf (ones (2), ones (3), ones(2))
+%!error nakapdf (ones (2), ones (2), ones(3))
+%!error nakapdf (i, 2, 2)
+%!error nakapdf (2, i, 2)
+%!error nakapdf (2, 2, i)
+
diff -r 1449e3b98941 -r f047e12820ae scripts/statistics/distributions/nakarnd.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/statistics/distributions/nakarnd.m	Fri May 13 23:49:24 2016 +0200
@@ -0,0 +1,137 @@
+## Copyright (C) 2016 Dag Lyberg
+## Copyright (C) 1995-2015 Kurt Hornik
+##
+## This file is part of Octave.
+##
+## Octave 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  {} {} nakarnd (@var{m}, @var{w})
+## @deftypefnx {} {} nakarnd (@var{m}, @var{w}, @var{r})
+## @deftypefnx {} {} nakarnd (@var{m}, @var{w}, @var{r}, @var{c}, @dots{})
+## @deftypefnx {} {} nakarnd (@var{m}, @var{w}, [@var{sz}])
+## Return a matrix of random samples from the Nakagami distribution with
+## shape parameter @var{m} and scale @var{w}.
+##
+## When called with a single size argument, return 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 of dimensions @var{sz}.
+##
+## If no size arguments are given then the result matrix is the common size of
+## @var{m} and @var{w}.
+## @end deftypefn
+
+## Author: Dag Lyberg <daglyberg80@gmail.com>
+## Description: Random deviates from the Nakagami distribution
+
+function rnd = nakarnd (m, w, varargin)
+
+  if (nargin < 2)
+    print_usage ();
+  endif
+
+  if (! isscalar (m) || ! isscalar (w))
+    [retval, m, w] = common_size (m, w);
+    if (retval > 0)
+      error ("nakarnd: M and W must be of common size or scalars");
+    endif
+  endif
+
+  if (iscomplex (m) || iscomplex (w))
+    error ("nakarnd: M and W must not be complex");
+  endif
+
+  if (nargin == 2)
+    sz = size (m);
+  elseif (nargin == 3)
+    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 ("nakarnd: dimension vector must be row vector of non-negative integers");
+    endif
+  elseif (nargin > 3)
+    if (any (cellfun (@(x) (! isscalar (x) || x < 0), varargin)))
+      error ("nakarnd: dimensions must be non-negative integers");
+    endif
+    sz = [varargin{:}];
+  endif
+
+  if (! isscalar (m) && ! isequal (size (w), sz))
+    error ("nakagrnd: M and W must be scalar or of size SZ");
+  endif
+
+  if (isa (m, "single") || isa (w, "single"))
+    cls = "single";
+  else
+    cls = "double";
+  endif
+
+  if (isscalar (m) && isscalar (w))
+    if ((0 < m) && (m < Inf) && (0 < w) && (w < Inf))
+      m_gamma = m;
+      w_gamma = w/m;
+      rnd = gamrnd(m_gamma, w_gamma, sz);
+      rnd = sqrt(rnd);
+    else
+      rnd = NaN (sz, cls);
+    endif
+  else
+    rnd = NaN (sz, cls);
+    k = (0 < m) & (m < Inf) & (0 < w) & (w < Inf);
+    m_gamma = m;
+    w_gamma = w./m;
+    rnd(k) = gamrnd(m_gamma(k), w_gamma(k));
+    rnd(k) = sqrt(rnd(k));
+  endif
+
+endfunction
+
+
+%!assert (size (nakarnd (1,1)), [1, 1])
+%!assert (size (nakarnd (ones (2,1), 1)), [2, 1])
+%!assert (size (nakarnd (ones (2,2), 1)), [2, 2])
+%!assert (size (nakarnd (1, ones (2,1))), [2, 1])
+%!assert (size (nakarnd (1, ones (2,2))), [2, 2])
+%!assert (size (nakarnd (1,1, 3)), [3, 3])
+%!assert (size (nakarnd (1,1, [4 1])), [4, 1])
+%!assert (size (nakarnd (1,1, 4, 1)), [4, 1])
+
+## Test class of input preserved
+%!assert (class (nakarnd (1,1)), "double")
+%!assert (class (nakarnd (single (1),1)), "single")
+%!assert (class (nakarnd (single ([1 1]),1)), "single")
+%!assert (class (nakarnd (1,single (1))), "single")
+%!assert (class (nakarnd (1,single ([1 1]))), "single")
+
+## Test input validation
+%!error nakarnd ()
+%!error nakarnd (1)
+%!error nakarnd (zeros (3), ones (2))
+%!error nakarnd (zeros (2), ones (3))
+%!error nakarnd (i, 2)
+%!error nakarnd (1, i)
+%!error nakarnd (1,2, -1)
+%!error nakarnd (1,2, ones (2))
+%!error nakarnd (1, 2, [2 -1 2])
+%!error nakarnd (1,2, 1, ones (2))
+%!error nakarnd (1,2, 1, -1)
+%!error nakarnd (ones (2,2), 2, 3)
+%!error nakarnd (ones (2,2), 2, [3, 2])
+%!error nakarnd (ones (2,2), 2, 2, 3)
+
diff -r 1449e3b98941 -r f047e12820ae scripts/statistics/distributions/raylcdf.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/statistics/distributions/raylcdf.m	Fri May 13 23:49:24 2016 +0200
@@ -0,0 +1,89 @@
+## Copyright (C) 2016 Dag Lyberg
+## Copyright (C) 1995-2015 Kurt Hornik
+##
+## This file is part of Octave.
+##
+## Octave 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 {} {} raylcdf (@var{x}, @var{scale})
+## For each element of @var{x}, compute the cumulative distribution function
+## (CDF) at @var{x} of the Rayleigh distribution with parameter @var{scale}.
+## @end deftypefn
+
+## Author: Dag Lyberg <daglyberg80@gmail.com>
+## Description: CDF of the Rayleigh distribution
+
+function cdf = raylcdf (x, scale)
+
+  if (nargin != 2)
+    print_usage ();
+  endif
+
+  if (! isscalar (scale))
+    [retval, x, scale] = common_size (x, scale);
+    if (retval > 0)
+      error ("raylcdf: X and SCALE must be of common size or scalars");
+    endif
+  endif
+
+  if (iscomplex (x) || iscomplex (scale))
+    error ("raylcdf: X and SCALE must not be complex");
+  endif
+
+  if (isa (x, "single") || isa (scale, "single"))
+    cdf = zeros (size (x), "single");
+  else
+    cdf = zeros (size (x));
+  endif
+
+  k = isnan (x) | ! (scale > 0);
+  cdf(k) = NaN;
+
+  k = (x == Inf) & (scale > 0);
+  cdf(k) = 1;
+
+  k = (0 <= x) & (x < Inf) & (0 < scale) & (scale < Inf);
+  if (isscalar(scale))
+    cdf(k) = 1 - exp(-x(k).^2/(2*scale^2));
+  else
+    cdf(k) = 1 - exp(-x(k).^2./(2*scale(k).^2));
+  endif
+
+endfunction
+
+
+%!shared x,y
+%! x = [-1, 0, 1, 2, Inf];
+%! y = [0, 0, 0.39346934028736658, 0.8646647167633873, 1];
+%!assert (raylcdf (x, ones (1,5)), y)
+%!assert (raylcdf (x, 1), y)
+%!assert (raylcdf (x, 1*[1 1 NaN 1 1]), [y(1:2) NaN y(4:5)])
+%!assert (raylcdf ([x, NaN], 1), [y, NaN])
+
+## Test class of input preserved
+%!assert (raylcdf (x, 1), double (y))
+%!assert (raylcdf (single ([x, NaN]), 1), single ([y, NaN]))
+%!assert (raylcdf ([x, NaN], single (1)), single ([y, NaN]))
+
+## Test input validation
+%!error raylcdf ()
+%!error raylcdf (1)
+%!error raylcdf (1,2,3)
+%!error raylcdf (ones (3), ones (2))
+%!error raylcdf (ones (2), ones (3))
+%!error raylcdf (i, 2)
+%!error raylcdf (2, i)
+
diff -r 1449e3b98941 -r f047e12820ae scripts/statistics/distributions/raylinv.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/statistics/distributions/raylinv.m	Fri May 13 23:49:24 2016 +0200
@@ -0,0 +1,89 @@
+## Copyright (C) 2016 Dag Lyberg
+## Copyright (C) 1995-2015 Kurt Hornik
+##
+## This file is part of Octave.
+##
+## Octave 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 {} {} raylinv (@var{x}, @var{scale})
+## For each element of @var{x}, compute the quantile (the inverse of the CDF)
+## at @var{x} of the Rayleigh distribution with parameter @var{scale}.
+## @end deftypefn
+
+## Author: Dag Lyberg <daglyberg80@gmail.com>
+## Description: CDF of the Rayleigh distribution
+
+function inv = raylinv (x, scale)
+
+  if (nargin != 2)
+    print_usage ();
+  endif
+
+  if (! isscalar (scale))
+    [retval, x, scale] = common_size (x, scale);
+    if (retval > 0)
+      error ("raylinv: X and SCALE must be of common size or scalars");
+    endif
+  endif
+
+  if (iscomplex (x) || iscomplex (scale))
+    error ("raylinv: X and SCALE must not be complex");
+  endif
+
+  if (isa (x, "single") || isa (scale, "single"))
+    inv = zeros (size (x), "single");
+  else
+    inv = zeros (size (x));
+  endif
+
+  k = isnan (x) | ! (0 <= x) | ! (x <= 1) | ! (scale > 0);
+  inv(k) = NaN;
+
+  k = (x == 1) & (scale > 0);
+  inv(k) = Inf;
+
+  k = (0 <= x) & (x < 1) & (0 < scale) & (scale < Inf);
+  if (isscalar(scale))
+    inv(k) = sqrt (-2*scale.*log(1-x(k)));
+  else
+    inv(k) = sqrt (-2*scale(k).*log(1-x(k)));
+  endif
+
+endfunction
+
+
+%!shared x,y
+%! x = [-1, 0, 0.5, 1, 2];
+%! y = [NaN, 0, 1.1774100225154747, Inf, NaN];
+%!assert (raylinv (x, ones (1,5)), y)
+%!assert (raylinv (x, 1), y)
+%!assert (raylinv (x, [1, 1, NaN, 1, 1]), [y(1:2), NaN, y(4:5)])
+%!assert (raylinv ([x, NaN], 1), [y, NaN])
+
+## Test class of input preserved
+%!assert (raylinv (x, 1), double (y))
+%!assert (raylinv (single ([x, NaN]), 1), single ([y, NaN]))
+%!assert (raylinv ([x, NaN], single (1)), single ([y, NaN]))
+
+## Test input validation
+%!error raylinv ()
+%!error raylinv (1)
+%!error raylinv (1,2,3)
+%!error raylinv (ones (3), ones (2))
+%!error raylinv (ones (2), ones (3))
+%!error raylinv (i, 2)
+%!error raylinv (2, i)
+
diff -r 1449e3b98941 -r f047e12820ae scripts/statistics/distributions/raylpdf.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/statistics/distributions/raylpdf.m	Fri May 13 23:49:24 2016 +0200
@@ -0,0 +1,86 @@
+## Copyright (C) 2016 Dag Lyberg
+## Copyright (C) 1995-2015 Kurt Hornik
+##
+## This file is part of Octave.
+##
+## Octave 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 {} {} raylpdf (@var{x}, @var{scale})
+## For each element of @var{x}, compute the probability density function (PDF)
+## at @var{x} of the Rayleigh distribution with parameter @var{scale}.
+## @end deftypefn
+
+## Author: Dag Lyberg <daglyberg80@gmail.com>
+## Description: PDF of the Rayleigh distribution
+
+function pdf = raylpdf (x, scale)
+
+  if (nargin != 2)
+    print_usage ();
+  endif
+
+  if (! isscalar (scale))
+    [retval, x, scale] = common_size (x, scale);
+    if (retval > 0)
+      error ("raylpdf: X and SCALE must be of common size or scalars");
+    endif
+  endif
+
+  if (iscomplex (x) || iscomplex (scale))
+    error ("raylpdf: X and SCALE must not be complex");
+  endif
+
+  if (isa (x, "single") || isa (scale, "single"))
+    pdf = zeros (size (x), "single");
+  else
+    pdf = zeros (size (x));
+  endif
+
+  k = isnan (x) | ! (scale > 0);
+  pdf(k) = NaN;
+
+  k = (0 <= x) & (x < Inf) & (0 < scale) & (scale < Inf);
+  if (isscalar(scale))
+    pdf(k) = (x(k)/scale^2) .* exp(-x(k).^2/(2*scale^2));
+  else
+    pdf(k) = (x(k)./scale(k).^2) .* exp(-x(k).^2./(2*scale(k).^2));
+  endif
+
+endfunction
+
+
+%!shared x,y
+%! x = [-1, 0, 1, 2, Inf];
+%! y = [0, 0, 0.60653065971263342, 0.2706705664732254, 0];
+%!assert (raylpdf (x, ones (1,5)), y)
+%!assert (raylpdf (x, 1), y)
+%!assert (raylpdf (x, [1 1 NaN 1 1]), [y(1:2) NaN y(4:5)])
+%!assert (raylpdf ([x, NaN], 1), [y, NaN])
+
+## Test class of input preserved
+%!assert (raylpdf (x, 1), double (y))
+%!assert (raylpdf (single ([x, NaN]), 1), single ([y, NaN]))
+%!assert (raylpdf ([x, NaN], single (1)), single ([y, NaN]))
+
+## Test input validation
+%!error raylpdf ()
+%!error raylpdf (1)
+%!error raylpdf (1,2,3)
+%!error raylpdf (ones (3), ones (2))
+%!error raylpdf (ones (2), ones (3))
+%!error raylpdf (i, 2)
+%!error raylpdf (2, i)
+
diff -r 1449e3b98941 -r f047e12820ae scripts/statistics/distributions/raylrnd.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/statistics/distributions/raylrnd.m	Fri May 13 23:49:24 2016 +0200
@@ -0,0 +1,118 @@
+## Copyright (C) 2016 Dag Lyberg
+## Copyright (C) 1995-2015 Kurt Hornik
+##
+## This file is part of Octave.
+##
+## Octave 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  {} {} trirnd (@var{scale})
+## @deftypefnx {} {} trirnd (@var{scale}, @var{r})
+## @deftypefnx {} {} trirnd (@var{scale}, @var{r}, @var{c}, @dots{})
+## @deftypefnx {} {} trirnd (@var{scale}, [@var{sz}])
+## Return a matrix of random samples from the Rayleigh distribution.
+##
+## When called with a single size argument, return 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 of dimensions @var{sz}.
+##
+## If no size arguments are given then the result matrix is the size of
+## @var{scale}.
+## @end deftypefn
+
+## Author: Dag Lyberg <daglyberg80@gmail.com>
+## Description: Random deviates from the Rayleigh distribution
+
+function rnd = raylrnd (scale, varargin)
+
+  if (nargin < 1)
+    print_usage ();
+  endif
+
+  if (iscomplex (scale))
+    error ("raylrnd: SCALE must not be complex");
+  endif
+
+  if (nargin == 1)
+    sz = size (scale);
+  elseif (nargin == 2)
+    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 ("raylrnd: dimension vector must be row vector of non-negative integers");
+    endif
+  elseif (nargin > 2)
+    if (any (cellfun (@(x) (! isscalar (x) || x < 0), varargin)))
+      error ("raylrnd: dimensions must be non-negative integers");
+    endif
+    sz = [varargin{:}];
+  endif
+
+  if (! isscalar (scale) && ! isequal (size (scale), sz))
+    error ("raylrnd: SCALE must be scalar or of size SZ");
+  endif
+
+  if (isa (scale, "single"))
+    cls = "single";
+  else
+    cls = "double";
+  endif
+
+  if (isscalar (scale))
+    if ((0 < scale) && (scale < Inf))
+      rnd  = rand(sz,cls);
+      rnd(:) = scale * sqrt (-2 * log (1 - rnd(:)));
+    else
+      rnd = NaN (sz, cls);
+    endif
+  else
+    rnd = NaN (sz, cls);
+
+    k = (0 < scale) & (scale < Inf);
+    rnd(k) = rand(sum(k(:)),1);
+    rnd(k) = scale(k) .* sqrt (-2 * log (1 - rnd(k)));
+  endif
+
+endfunction
+
+
+%!assert (size (raylrnd (1)), [1 1])
+%!assert (size (raylrnd (ones (2,1))), [2, 1])
+%!assert (size (raylrnd (ones (2,2))), [2, 2])
+%!assert (size (raylrnd (1, 3)), [3, 3])
+%!assert (size (raylrnd (1, [4 1])), [4, 1])
+%!assert (size (raylrnd (1, 4, 1)), [4, 1])
+
+## Test class of input preserved
+%!assert (class (raylrnd (1)), "double")
+%!assert (class (raylrnd (single (1))), "single")
+%!assert (class (raylrnd (single ([1 1]))), "single")
+
+## Test input validation
+%!error raylrnd ()
+%!error raylrnd (1, -1)
+%!error raylrnd (1, ones (2))
+%!error raylrnd (i)
+%!error raylrnd (1, [2 -1 2])
+%!error raylrnd (1, 2, -1)
+%!error raylrnd (1, 2, ones (2))
+%!error raylrnd (ones (2,2), 3)
+%!error raylrnd (ones (2,2), [3, 2])
+%!error raylrnd (ones (2,2), 2, 3)
+
diff -r 1449e3b98941 -r f047e12820ae scripts/statistics/distributions/tricdf.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/statistics/distributions/tricdf.m	Fri May 13 23:49:24 2016 +0200
@@ -0,0 +1,131 @@
+## Copyright (C) 2016 Dag Lyberg
+## Copyright (C) 1997-2015 Kurt Hornik
+##
+## This file is part of Octave.
+##
+## Octave 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 {} {} tricdf (@var{x}, @var{a}, @var{b}, @var{c})
+## Compute the cumulative distribution function (CDF) at @var{x} of the
+## triangular distribution with parameters @var{a}, @var{b}, and @var{c}
+## on the interval [@var{a}, @var{b}].
+## @end deftypefn
+
+## Author: Dag Lyberg <daglyberg80@gmail.com>
+## Description: CDF of the triangle distribution
+
+function cdf = tricdf (x, a, b, c)
+
+  if (nargin != 4)
+    print_usage ();
+  endif
+
+  if (! isscalar (a) || ! isscalar (b) || ! isscalar (c))
+    [retval, x, a, b, c] = common_size (x, a, b, c);
+    if (retval > 0)
+      error ("tricdf: X, A, B, and C must be of common size or scalars");
+    endif
+  endif
+
+  if (iscomplex (x) || iscomplex (a) || iscomplex (b) || iscomplex (c))
+    error ("tricdf: X, A, B, and C must not be complex");
+  endif
+
+  if (isa (x, "single") || isa (a, "single")
+      || isa (b, "single") || isa (c, "single"))
+    cdf = zeros (size (x), "single");
+  else
+    cdf = zeros (size (x));
+  endif
+
+  k = isnan (x) | !(a < b) | !(c >= a) | !(c <= b) ;
+  cdf(k) = NaN;
+
+  k = (x > a) & (a < b) & (a <= c) & (c <= b);
+  if (isscalar (a) && isscalar (b) && isscalar (c))
+    h = 2 / (b-a);
+
+    k_temp = k & (c <= x);
+    full_area = (c-a) * h / 2;
+    cdf(k_temp) += full_area;
+
+    k_temp = k & (a < x) & (x < c);
+    area = (x(k_temp) - a).^2 * h;
+    cdf(k_temp) += area;
+
+    k_temp = k & (b <= x);
+    full_area = (b-c) * h / 2;
+    cdf(k_temp) += full_area;
+
+    k_temp = k & (c < x) & (x < b);
+    area = (b-x(k_temp)).^2 * h;
+    cdf(k_temp) += full_area - area;
+  else
+    h = 2 ./ (b-a);
+
+    k_temp = k & (c <= x);
+    full_area = (c(k_temp)-a(k_temp)) .* h(k_temp) / 2;
+    cdf(k_temp) += full_area;
+
+    k_temp = k & (a <= x) & (x < c);
+    area = (x(k_temp) - a(k_temp)).^2 .* h(k_temp);
+    cdf(k_temp) += area;
+
+    k_temp = k & (b <= x);
+    full_area = (b(k_temp)-c(k_temp)) .* h(k_temp) / 2;
+    cdf(k_temp) += full_area;
+
+    k_temp = k & (c <= x) & (x < b);
+    area = (b(k_temp)-x(k_temp)).^2 .* h(k_temp);
+    cdf(k_temp) += full_area - area;
+  endif
+
+endfunction
+
+
+%!shared x,y
+%! x = [-1, 0, 0.1, 0.5, 0.9, 1, 2] + 1;
+%! y = [0, 0, 0.02, 0.5, 0.98, 1 1];
+%!assert (tricdf (x, ones (1,7), 2*ones (1,7), 1.5*ones (1,7)), y, eps)
+%!assert (tricdf (x, 1*ones (1,7), 2, 1.5), y, eps)
+%!assert (tricdf (x, 1, 2*ones (1,7), 1.5), y, eps)
+%!assert (tricdf (x, 1, 2, 1.5*ones (1,7)), y, eps)
+%!assert (tricdf (x, 1, 2, 1.5), y, eps)
+%!assert (tricdf (x, [1, 1, NaN, 1, 1, 1, 1], 2, 1.5), [y(1:2), NaN, y(4:7)], eps)
+%!assert (tricdf (x, 1, 2*[1, 1, NaN, 1, 1, 1, 1], 1.5), [y(1:2), NaN, y(4:7)], eps)
+%!assert (tricdf (x, 1, 2, 1.5*[1, 1, NaN, 1, 1, 1, 1]), [y(1:2), NaN, y(4:7)], eps)
+%!assert (tricdf ([x, NaN], 1, 2, 1.5), [y, NaN], eps)
+
+## Test class of input preserved
+%!assert (tricdf (single ([x, NaN]), 1, 2, 1.5), single ([y, NaN]), eps('single'))
+%!assert (tricdf ([x, NaN], single (1), 2, 1.5), single ([y, NaN]), eps('single'))
+%!assert (tricdf ([x, NaN], 1, single (2), 1.5), single ([y, NaN]), eps('single'))
+%!assert (tricdf ([x, NaN], 1, 2, single (1.5)), single ([y, NaN]), eps('single'))
+
+## Test input validation
+%!error tricdf ()
+%!error tricdf (1)
+%!error tricdf (1,2)
+%!error tricdf (1,2,3)
+%!error tricdf (1,2,3,4,5)
+%!error tricdf (1, ones (3), ones (2), ones (2))
+%!error tricdf (1, ones (2), ones (3), ones (2))
+%!error tricdf (1, ones (2), ones (2), ones (3))
+%!error tricdf (i, 2, 2, 2)
+%!error tricdf (2, i, 2, 2)
+%!error tricdf (2, 2, i, 2)
+%!error tricdf (2, 2, 2, i)
+
diff -r 1449e3b98941 -r f047e12820ae scripts/statistics/distributions/triinv.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/statistics/distributions/triinv.m	Fri May 13 23:49:24 2016 +0200
@@ -0,0 +1,119 @@
+## Copyright (C) 2016 Dag Lyberg
+## Copyright (C) 1995-2015 Kurt Hornik
+##
+## This file is part of Octave.
+##
+## Octave 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 {} {} triinv (@var{x}, @var{a}, @var{b}, @var{c})
+## For each element of @var{x}, compute the quantile (the inverse of the CDF)
+## at @var{x} of the triangular distribution with parameters
+## @var{a}, @var{b}, and @var{c} on the interval [@var{a}, @var{b}].
+## @end deftypefn
+
+## Author: Dag Lyberg <daglyberg80@gmail.com>
+## Description: Quantile function of the triangular distribution
+
+function inv = triinv (x, a, b, c)
+
+  if (nargin != 4)
+    print_usage ();
+  endif
+
+  if (! isscalar (a) || ! isscalar (b) || ! isscalar (c))
+    [retval, x, a, b, c] = common_size (x, a, b, c);
+    if (retval > 0)
+      error ("triinv: X, A, B, and C must be of common size or scalars");
+    endif
+  endif
+
+  if (iscomplex (x) || iscomplex (a) || iscomplex (b) || iscomplex (c))
+    error ("triinv: X, A, B, and C must not be complex");
+  endif
+
+  if (isa (x, "single") || isa (a, "single") || isa (b, "single"))
+    inv = NaN (size (x), "single");
+  else
+    inv = NaN (size (x));
+  endif
+
+  k = (x >= 0) & (x <= 1) & (a < b) & (a <= c) & (c <= b);
+  inv(k) = 0;
+  
+  if (isscalar (a) && isscalar (b) && isscalar(c))
+    h = 2 / (b-a);
+    w = c-a;
+    area1 = h * w / 2;
+    j = k & (x <= area1);
+    inv(j) += (x(j) * (h/2) * w).^0.5 + a;
+
+    w = b-c;
+    j = k & (area1 < x) & (x < 1);
+    inv(j) += b - ((1-x(j)) * (h/2) * w).^0.5;
+
+    j = k & (x == 1);
+    inv(j) = b;
+  else
+    h = 2 ./ (b-a);
+    w = c-a;
+    area1 = h .* w / 2;
+    j = k & (x <= area1);
+    inv(j) += (x(j) .* (h(j)/2) .* w(j)).^0.5 + a(j);
+
+    w = b-c;
+    j = k & (area1 < x) & (x < 1);
+    inv(j) += b(j) - ((1-x(j)) .* (h(j)/2) .* w(j)).^0.5;
+
+    j = k & (x == 1);
+    inv(j) = b(j);
+  endif
+
+endfunction
+
+
+%!shared x,y
+%! x = [-1, 0, 0.02, 0.5, 0.98, 1, 2];
+%! y = [NaN, 0, 0.1, 0.5, 0.9, 1, NaN] + 1;
+%!assert (triinv (x, ones (1,7), 2*ones (1,7), 1.5*ones (1,7)), y, eps)
+%!assert (triinv (x, 1*ones (1,7), 2, 1.5), y, eps)
+%!assert (triinv (x, 1, 2*ones (1,7), 1.5), y, eps)
+%!assert (triinv (x, 1, 2, 1.5*ones (1,7)), y, eps)
+%!assert (triinv (x, 1, 2, 1.5), y, eps)
+%!assert (triinv (x, [1, 1, NaN, 1, 1, 1, 1], 2, 1.5), [y(1:2), NaN, y(4:7)], eps)
+%!assert (triinv (x, 1, 2*[1, 1, NaN, 1, 1, 1, 1], 1.5), [y(1:2), NaN, y(4:7)], eps)
+%!assert (triinv (x, 1, 2, 1.5*[1, 1, NaN, 1, 1, 1, 1]), [y(1:2), NaN, y(4:7)], eps)
+%!assert (triinv ([x, NaN], 1, 2, 1.5), [y, NaN], eps)
+
+## Test class of input preserved
+%!assert (triinv (single ([x, NaN]), 1, 2, 1.5), single ([y, NaN]), eps('single'))
+%!assert (triinv ([x, NaN], single (1), 2, 1.5), single ([y, NaN]), eps('single'))
+%!assert (triinv ([x, NaN], 1, single (2), 1.5), single ([y, NaN]), eps('single'))
+%!assert (triinv ([x, NaN], 1, 2, single (1.5)), single ([y, NaN]), eps('single'))
+
+## Test input validation
+%!error triinv ()
+%!error triinv (1)
+%!error triinv (1,2)
+%!error triinv (1,2,3)
+%!error triinv (1,2,3,4,5)
+%!error triinv (1, ones (3), ones (2), ones (2))
+%!error triinv (1, ones (2), ones (3), ones (2))
+%!error triinv (1, ones (2), ones (2), ones (3))
+%!error triinv (i, 2, 2, 2)
+%!error triinv (2, i, 2, 2)
+%!error triinv (2, 2, i, 2)
+%!error triinv (2, 2, 2, i)
+
diff -r 1449e3b98941 -r f047e12820ae scripts/statistics/distributions/tripdf.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/statistics/distributions/tripdf.m	Fri May 13 23:49:24 2016 +0200
@@ -0,0 +1,111 @@
+## Copyright (C) 2016 Dag Lyberg
+## Copyright (C) 1997-2015 Kurt Hornik
+##
+## This file is part of Octave.
+##
+## Octave 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 {} {} tripdf (@var{x}, @var{a}, @var{b}, @var{c})
+## Compute the probability density function (PDF) at @var{x} of the triangular
+## distribution with parameters @var{a}, @var{b}, and @var{c} on the interval
+## [@var{a}, @var{b}].
+## @end deftypefn
+
+## Author: Dag Lyberg <daglyberg80@gmail.com>
+## Description: PDF of the triangular distribution
+
+function pdf = tripdf (x, a, b, c)
+
+  if (nargin != 4)
+    print_usage ();
+  endif
+
+  if (! isscalar (a) || ! isscalar (b) || ! isscalar (c))
+    [retval, x, a, b, c] = common_size (x, a, b, c);
+    if (retval > 0)
+      error ("tripdf: X, A, B, and C must be of common size or scalars");
+    endif
+  endif
+
+  if (iscomplex (x) || iscomplex (a) || iscomplex (b) || iscomplex (c))
+    error ("tripdf: X, A, B, and C must not be complex");
+  endif
+
+  if (isa (x, "single") || isa (a, "single") ...
+      || isa (b, "single") || isa (c, "single"))
+    pdf = zeros (size (x), "single");
+  else
+    pdf = zeros (size (x));
+  endif
+
+  k = isnan (x) | !(a < b) | !(c >= a) | !(c <= b) ;
+  pdf(k) = NaN;
+
+  k = (x >= a) & (x <= b) & (a < b) & (a <= c) & (c <= b);
+  h = 2 ./ (b-a);
+  if (isscalar (a) && isscalar (b) && isscalar (c))
+    j = k & (a <= x) & (x < c);
+    pdf(j) = h * (x(j)-a) / (c-a);
+    j = k & (x == c);
+    pdf(j) = h;
+    j = k & (c < x) & (x <= b);
+    pdf(j) = h * (b-x(j)) / (b-c);
+  else
+    j = k & (a <= x) & (x < c);
+    pdf(j) = h(j) .* (x(j)-a(j)) ./ (c(j)-a(j));
+    j = k & (x == c);
+    pdf(j) = h(j);
+    j = k & (c < x) & (x <= b);
+    pdf(j) = h(j) .* (b(j)-x(j)) ./ (b(j)-c(j));
+  endif
+
+endfunction
+
+
+%!shared x,y,deps
+%! x = [-1, 0, 0.1, 0.5, 0.9, 1, 2] + 1;
+%! y = [0, 0, 0.4, 2, 0.4, 0, 0];
+%! deps = 2*eps;
+%!assert (tripdf (x, ones (1,7), 2*ones (1,7), 1.5*ones (1,7)), y, deps)
+%!assert (tripdf (x, 1*ones (1,7), 2, 1.5), y, deps)
+%!assert (tripdf (x, 1, 2*ones (1,7), 1.5), y, deps)
+%!assert (tripdf (x, 1, 2, 1.5*ones (1,7)), y, deps)
+%!assert (tripdf (x, 1, 2, 1.5), y, deps)
+%!assert (tripdf (x, [1, 1, NaN, 1, 1, 1, 1], 2, 1.5), [y(1:2), NaN, y(4:7)], deps)
+%!assert (tripdf (x, 1, 2*[1, 1, NaN, 1, 1, 1, 1], 1.5), [y(1:2), NaN, y(4:7)], deps)
+%!assert (tripdf (x, 1, 2, 1.5*[1, 1, NaN, 1, 1, 1, 1]), [y(1:2), NaN, y(4:7)], deps)
+%!assert (tripdf ([x, NaN], 1, 2, 1.5), [y, NaN], deps)
+
+## Test class of input preserved
+%!assert (tripdf (single ([x, NaN]), 1, 2, 1.5), single ([y, NaN]), eps('single'))
+%!assert (tripdf ([x, NaN], single (1), 2, 1.5), single ([y, NaN]), eps('single'))
+%!assert (tripdf ([x, NaN], 1, single (2), 1.5), single ([y, NaN]), eps('single'))
+%!assert (tripdf ([x, NaN], 1, 2, single (1.5)), single ([y, NaN]), eps('single'))
+
+## Test input validation
+%!error tripdf ()
+%!error tripdf (1)
+%!error tripdf (1,2)
+%!error tripdf (1,2,3)
+%!error tripdf (1,2,3,4,5)
+%!error tripdf (1, ones (3), ones (2), ones (2))
+%!error tripdf (1, ones (2), ones (3), ones (2))
+%!error tripdf (1, ones (2), ones (2), ones (3))
+%!error tripdf (i, 2, 2, 2)
+%!error tripdf (2, i, 2, 2)
+%!error tripdf (2, 2, i, 2)
+%!error tripdf (2, 2, 2, i)
+
diff -r 1449e3b98941 -r f047e12820ae scripts/statistics/distributions/trirnd.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/statistics/distributions/trirnd.m	Fri May 13 23:49:24 2016 +0200
@@ -0,0 +1,153 @@
+## Copyright (C) 2016 Dag Lyberg
+## Copyright (C) 1995-2015 Kurt Hornik
+##
+## This file is part of Octave.
+##
+## Octave 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  {} {} trirnd (@var{a}, @var{b}, @var{c})
+## @deftypefnx {} {} trirnd (@var{a}, @var{b}, @var{c}, @var{r})
+## @deftypefnx {} {} trirnd (@var{a}, @var{b}, @var{c}, @var{r}, @var{c}, @dots{})
+## @deftypefnx {} {} trirnd (@var{a}, @var{b}, @var{c}, [@var{sz}])
+## Return a matrix of random samples from the rectangular distribution with
+## parameters @var{a}, @var{b}, and @var{c} on the interval [@var{a}, @var{b}].
+##
+## When called with a single size argument, return 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 of dimensions @var{sz}.
+##
+## If no size arguments are given then the result matrix is the common size of
+## @var{a}, @var{b} and @var{c}.
+## @end deftypefn
+
+## Author: Dag Lyberg <daglyberg80@gmail.com>
+## Description: Random deviates from the triangular distribution
+
+function rnd = trirnd (a, b, c, varargin)
+
+  if (nargin < 3)
+    print_usage ();
+  endif
+
+  if (! isscalar (a) || ! isscalar (b) || ! isscalar (c))
+    [retval, a, b, c] = common_size (a, b, c);
+    if (retval > 0)
+      error ("trirnd: A, B, and C must be of common size or scalars");
+    endif
+  endif
+
+  if (iscomplex (a) || iscomplex (b) || iscomplex (c))
+    error ("trirnd: A, B, and C must not be complex");
+  endif
+
+  if (nargin == 3)
+    sz = size (a);
+  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 ("trirnd: dimension vector must be row vector of non-negative integers");
+    endif
+  elseif (nargin > 4)
+    if (any (cellfun (@(x) (! isscalar (x) || x < 0), varargin)))
+      error ("trirnd: dimensions must be non-negative integers");
+    endif
+    sz = [varargin{:}];
+  endif
+
+  if (! isscalar (a) && ! isequal (size (b), sz))
+    error ("trirnd: A, B, and C must be scalar or of size SZ");
+  endif
+
+  if (isa (a, "single") || isa (b, "single") || isa (c, "single"))
+    cls = "single";
+  else
+    cls = "double";
+  endif
+
+  if (isscalar (a) && isscalar (b) && isscalar (c))
+    if ((-Inf < a) && (a < b) && (a <= c) && (c <= b) && (b < Inf))
+      w = b-a;
+      left_width = c-a;
+      right_width = b-c;
+      h = 2 / w;
+      left_area = h * left_width / 2;
+      rnd = rand (sz, cls);
+      idx = rnd < left_area;
+      rnd(idx) = a + (rnd(idx) * w * left_width).^0.5;
+      rnd(~idx) = b - ((1-rnd(~idx)) * w * right_width).^0.5;
+    else
+      rnd = NaN (sz, cls);
+    endif
+  else
+    w = b-a;
+    left_width = c-a;
+    right_width = b-c;
+    h = 2 ./ w;
+    left_area = h .* left_width / 2;
+    rnd = rand (sz, cls);
+    k = rnd < left_area;
+    rnd(k) = a(k) + (rnd(k) .* w(k) .* left_width(k)).^0.5;
+    rnd(~k) = b(~k) - ((1-rnd(~k)) .* w(~k) .* right_width(~k)).^0.5;
+
+    k = ! (-Inf < a) | ! (a < b) | ! (a <= c) | ! (c <= b) | ! (b < Inf);
+    rnd(k) = NaN;
+  endif
+
+endfunction
+
+
+%!assert (size (trirnd (1,2,1.5)), [1, 1])
+%!assert (size (trirnd (1*ones (2,1), 2,1.5)), [2, 1])
+%!assert (size (trirnd (1*ones (2,2), 2,1.5)), [2, 2])
+%!assert (size (trirnd (1, 2*ones (2,1), 1.5)), [2, 1])
+%!assert (size (trirnd (1, 2*ones (2,2), 1.5)), [2, 2])
+%!assert (size (trirnd (1, 2, 1.5*ones (2,1))), [2, 1])
+%!assert (size (trirnd (1, 2, 1.5*ones (2,2))), [2, 2])
+%!assert (size (trirnd (1, 2, 1.5, 3)), [3, 3])
+%!assert (size (trirnd (1, 2, 1.5, [4 1])), [4, 1])
+%!assert (size (trirnd (1, 2, 1.5, 4, 1)), [4, 1])
+
+## Test class of input preserved
+%!assert (class (trirnd (1,2,1.5)), "double")
+%!assert (class (trirnd (single (1),2,1.5)), "single")
+%!assert (class (trirnd (single ([1 1]),2,1.5)), "single")
+%!assert (class (trirnd (1,single (2),1.5)), "single")
+%!assert (class (trirnd (1,single ([2 2]),1.5)), "single")
+%!assert (class (trirnd (1,2,single (1.5))), "single")
+%!assert (class (trirnd (1,2,single ([1.5 1.5]))), "single")
+
+## Test input validation
+%!error trirnd ()
+%!error trirnd (1)
+%!error trirnd (1,2)
+%!error trirnd (ones (3), 2*ones (2), 1.5*ones (2), 2)
+%!error trirnd (ones (2), 2*ones (3), 1.5*ones (2), 2)
+%!error trirnd (ones (2), 2*ones (2), 1.5*ones (3), 2)
+%!error trirnd (i, 2, 1.5)
+%!error trirnd (1, i, 1.5)
+%!error trirnd (1, 2, i)
+%!error trirnd (1,2,1.5, -1)
+%!error trirnd (1,2,1.5, ones (2))
+%!error trirnd (1,2,1.5, [2 -1 2])
+%!error trirnd (1*ones (2),2,1.5, 3)
+%!error trirnd (1*ones (2),2,1.5, [3, 2])
+%!error trirnd (1*ones (2),2,1.5, 3, 2)
+
