Rayleigh distribution functions
Arno Onken <[email protected]> Sun, 01 Oct 2006 21:50:33 +0200
| Newsgroups | gmane.comp.gnu.octave.sources |
|---|---|
| Message-ID | <[email protected]> |
This is a multi-part message in MIME format.
--------------090603030001070304060002
Content-Type: text/plain; charset=ISO-8859-1; format=flowed
Content-Transfer-Encoding: 7bit
Hi,
I implemented a few functions for the Rayleigh distribution and
committed them to Octave-Forge. But I noticed that most of the
distribution functions are actually core functions, so Octave itself
might be the better place for them. The functions are compatible with
MATLAB R14 (using hist for comparing raylrnd).
I know, citing Wikipedia is generally not a good idea. But I think it is
useful in this case.
Regards,
Arno
--------------090603030001070304060002
Content-Type: text/x-objcsrc;
name="raylcdf.m"
Content-Transfer-Encoding: 7bit
Content-Disposition: inline;
filename="raylcdf.m"
## Copyright (C) 2006 Arno Onken
##
## This program is free software; you can redistribute it and/or modify
## it under the terms of the GNU General Public License as published by
## the Free Software Foundation; either version 2 of the License, or
## (at your option) any later version.
##
## This program is distributed in the hope that it will be useful,
## but WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
## GNU General Public License for more details.
##
## You should have received a copy of the GNU General Public License
## along with this program; if not, write to the Free Software
## Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
## -*- texinfo -*-
## @deftypefn {Function File} {@var{p} =} raylcdf (@var{x}, @var{sigma})
## Calculates the cumulative distribution function of the Rayleigh
## distribution.
##
## Arguments are
##
## @itemize
## @item
## @var{x} is the support. The elements of @var{x} must be non-negative.
##
## @item
## @var{sigma} is the parameter of the Rayleigh distribution. The elements
## of @var{sigma} must be positive.
## @end itemize
## @var{x} and @var{sigma} must be of common size or one of them must be
## scalar.
##
## Return values are
##
## @itemize
## @item
## @var{p} is the cumulative distribution of the Rayleigh distribution at
## each element of @var{x} and corresponding parameter @var{sigma}.
## @end itemize
##
## Examples:
##
## @example
## x = 0:0.5:2.5;
## sigma = 1:6;
## p = raylcdf (x, sigma)
##
## p = raylcdf (x, 0.5)
## @end example
##
## References:
##
## @enumerate
## @item
## W. L. Martinez and A. R. Martinez. @cite{Computational Statistics
## Handbook with MATLAB.} Chapman & Hall/CRC, pages 547-557, 2001.
##
## @item
## Wikipedia contributors. Rayleigh distribution. @cite{Wikipedia, The Free
## Encyclopedia.}
## @uref{http://en.wikipedia.org/w/index.php?title=Rayleigh_distribution&oldid=69294908},
## August 2006.
## @end enumerate
## @end deftypefn
## Author: Arno Onken <[email protected]>
## Description: CDF of the Rayleigh distribution
function p = raylcdf (x, sigma)
# Check arguments
if (nargin != 2)
usage ("p = raylcdf (x, sigma)");
endif
if (! isempty (x) && ! ismatrix (x))
error ("raylcdf: x must be a numeric matrix");
endif
if (! isempty (sigma) && ! ismatrix (sigma))
error ("raylcdf: sigma must be a numeric matrix");
endif
if (! isscalar (x) || ! isscalar (sigma))
[retval, x, sigma] = common_size (x, sigma);
if (retval > 0)
error ("raylcdf: x and sigma must be of common size or scalar");
endif
endif
# Calculate cdf
p = 1 - exp ((-x .^ 2) ./ (2 * sigma .^ 2));
# Continue argument check
k = find (! (x >= 0) | ! (x < Inf) | ! (sigma > 0));
if (any (k))
p (k) = NaN;
endif
endfunction
%!test
%! x = 0:0.5:2.5;
%! sigma = 1:6;
%! p = raylcdf (x, sigma);
%! expected_p = [0.0000, 0.0308, 0.0540, 0.0679, 0.0769, 0.0831];
%! assert (p, expected_p, 0.001);
%!test
%! x = 0:0.5:2.5;
%! p = raylcdf (x, 0.5);
%! expected_p = [0.0000, 0.3935, 0.8647, 0.9889, 0.9997, 1.0000];
%! assert (p, expected_p, 0.001);
--------------090603030001070304060002
Content-Type: text/x-objcsrc;
name="raylinv.m"
Content-Transfer-Encoding: 7bit
Content-Disposition: inline;
filename="raylinv.m"
## Copyright (C) 2006 Arno Onken
##
## This program is free software; you can redistribute it and/or modify
## it under the terms of the GNU General Public License as published by
## the Free Software Foundation; either version 2 of the License, or
## (at your option) any later version.
##
## This program is distributed in the hope that it will be useful,
## but WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
## GNU General Public License for more details.
##
## You should have received a copy of the GNU General Public License
## along with this program; if not, write to the Free Software
## Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
## -*- texinfo -*-
## @deftypefn {Function File} {@var{x} =} raylinv (@var{p}, @var{sigma})
## Calculates the quantile of the Rayleigh distribution. The quantile is the
## inverse of the cumulative distribution function.
##
## Arguments are
##
## @itemize
## @item
## @var{p} is the cumulative distribution. The elements of @var{p} must be
## probabilities.
##
## @item
## @var{sigma} is the parameter of the Rayleigh distribution. The elements
## of @var{sigma} must be positive.
## @end itemize
## @var{p} and @var{sigma} must be of common size or one of them must be
## scalar.
##
## Return values are
##
## @itemize
## @item
## @var{x} is the quantile of the Rayleigh distribution at each element of
## @var{p} and corresponding parameter @var{sigma}.
## @end itemize
##
## Examples:
##
## @example
## p = 0:0.1:0.5;
## sigma = 1:6;
## x = raylinv (p, sigma)
##
## x = raylinv (p, 0.5)
## @end example
##
## References:
##
## @enumerate
## @item
## W. L. Martinez and A. R. Martinez. @cite{Computational Statistics
## Handbook with MATLAB.} Chapman & Hall/CRC, pages 547-557, 2001.
##
## @item
## Wikipedia contributors. Rayleigh distribution. @cite{Wikipedia, The Free
## Encyclopedia.}
## @uref{http://en.wikipedia.org/w/index.php?title=Rayleigh_distribution&oldid=69294908},
## August 2006.
## @end enumerate
## @end deftypefn
## Author: Arno Onken <[email protected]>
## Description: Quantile of the Rayleigh distribution
function x = raylinv (p, sigma)
# Check arguments
if (nargin != 2)
usage ("x = raylinv (p, sigma)");
endif
if (! isempty (p) && ! ismatrix (p))
error ("raylinv: p must be a numeric matrix");
endif
if (! isempty (sigma) && ! ismatrix (sigma))
error ("raylinv: sigma must be a numeric matrix");
endif
if (! isscalar (p) || ! isscalar (sigma))
[retval, p, sigma] = common_size (p, sigma);
if (retval > 0)
error ("raylinv: p and sigma must be of common size or scalar");
endif
endif
# Calculate quantile
x = sqrt (-2 .* log (1 - p) .* sigma .^ 2);
k = find (p == 1);
if (any (k))
x (k) = Inf;
endif
# Continue argument check
k = find (! (p >= 0) | ! (p <= 1) | ! (sigma > 0));
if (any (k))
x (k) = NaN;
endif
endfunction
%!test
%! p = 0:0.1:0.5;
%! sigma = 1:6;
%! x = raylinv (p, sigma);
%! expected_x = [0.0000, 0.9181, 2.0041, 3.3784, 5.0538, 7.0645];
%! assert (x, expected_x, 0.001);
%!test
%! p = 0:0.1:0.5;
%! x = raylinv (p, 0.5);
%! expected_x = [0.0000, 0.2295, 0.3340, 0.4223, 0.5054, 0.5887];
%! assert (x, expected_x, 0.001);
--------------090603030001070304060002
Content-Type: text/x-objcsrc;
name="raylpdf.m"
Content-Transfer-Encoding: 7bit
Content-Disposition: inline;
filename="raylpdf.m"
## Copyright (C) 2006 Arno Onken
##
## This program is free software; you can redistribute it and/or modify
## it under the terms of the GNU General Public License as published by
## the Free Software Foundation; either version 2 of the License, or
## (at your option) any later version.
##
## This program is distributed in the hope that it will be useful,
## but WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
## GNU General Public License for more details.
##
## You should have received a copy of the GNU General Public License
## along with this program; if not, write to the Free Software
## Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
## -*- texinfo -*-
## @deftypefn {Function File} {@var{y} =} raylpdf (@var{x}, @var{sigma})
## Calculates the probability density function of the Rayleigh distribution.
##
## Arguments are
##
## @itemize
## @item
## @var{x} is the support. The elements of @var{x} must be non-negative.
##
## @item
## @var{sigma} is the parameter of the Rayleigh distribution. The elements
## of @var{sigma} must be positive.
## @end itemize
## @var{x} and @var{sigma} must be of common size or one of them must be
## scalar.
##
## Return values are
##
## @itemize
## @item
## @var{y} is the probability density of the Rayleigh distribution at each
## element of @var{x} and corresponding parameter @var{sigma}.
## @end itemize
##
## Examples:
##
## @example
## x = 0:0.5:2.5;
## sigma = 1:6;
## y = raylpdf (x, sigma)
##
## y = raylpdf (x, 0.5)
## @end example
##
## References:
##
## @enumerate
## @item
## W. L. Martinez and A. R. Martinez. @cite{Computational Statistics
## Handbook with MATLAB.} Chapman & Hall/CRC, pages 547-557, 2001.
##
## @item
## Wikipedia contributors. Rayleigh distribution. @cite{Wikipedia, The Free
## Encyclopedia.}
## @uref{http://en.wikipedia.org/w/index.php?title=Rayleigh_distribution&oldid=69294908},
## August 2006.
## @end enumerate
## @end deftypefn
## Author: Arno Onken <[email protected]>
## Description: PDF of the Rayleigh distribution
function y = raylpdf (x, sigma)
# Check arguments
if (nargin != 2)
usage ("y = raylpdf (x, sigma)");
endif
if (! isempty (x) && ! ismatrix (x))
error ("raylpdf: x must be a numeric matrix");
endif
if (! isempty (sigma) && ! ismatrix (sigma))
error ("raylpdf: sigma must be a numeric matrix");
endif
if (! isscalar (x) || ! isscalar (sigma))
[retval, x, sigma] = common_size (x, sigma);
if (retval > 0)
error ("raylpdf: x and sigma must be of common size or scalar");
endif
endif
# Calculate pdf
y = x .* exp ((-x .^ 2) ./ (2 .* sigma .^ 2)) ./ (sigma .^ 2);
# Continue argument check
k = find (! (x >= 0) | ! (x < Inf) | ! (sigma > 0));
if (any (k))
y (k) = NaN;
endif
endfunction
%!test
%! x = 0:0.5:2.5;
%! sigma = 1:6;
%! y = raylpdf (x, sigma);
%! expected_y = [0.0000, 0.1212, 0.1051, 0.0874, 0.0738, 0.0637];
%! assert (y, expected_y, 0.001);
%!test
%! x = 0:0.5:2.5;
%! y = raylpdf (x, 0.5);
%! expected_y = [0.0000, 1.2131, 0.5413, 0.0667, 0.0027, 0.0000];
%! assert (y, expected_y, 0.001);
--------------090603030001070304060002
Content-Type: text/x-objcsrc;
name="raylrnd.m"
Content-Transfer-Encoding: 7bit
Content-Disposition: inline;
filename="raylrnd.m"
## Copyright (C) 2006 Arno Onken
##
## This program is free software; you can redistribute it and/or modify
## it under the terms of the GNU General Public License as published by
## the Free Software Foundation; either version 2 of the License, or
## (at your option) any later version.
##
## This program is distributed in the hope that it will be useful,
## but WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
## GNU General Public License for more details.
##
## You should have received a copy of the GNU General Public License
## along with this program; if not, write to the Free Software
## Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
## -*- texinfo -*-
## @deftypefn {Function File} {@var{x} =} raylrnd (@var{sigma})
## @deftypefnx {Function File} {@var{x} =} raylrnd (@var{sigma}, @var{sz})
## @deftypefnx {Function File} {@var{x} =} raylrnd (@var{sigma}, @var{r}, @var{c})
## Returns a matrix of random samples from the Rayleigh distribution.
##
## Arguments are
##
## @itemize
## @item
## @var{sigma} is the parameter of the Rayleigh distribution. The elements
## of @var{sigma} must be positive.
##
## @item
## @var{sz} is the size of the matrix to be generated. @var{sz} must be a
## vector of non-negative integers.
##
## @item
## @var{r} is the number of rows of the matrix to be generated. @var{r} must
## be a non-negative integer.
##
## @item
## @var{c} is the number of columns of the matrix to be generated. @var{c}
## must be a non-negative integer.
## @end itemize
##
## Return values are
##
## @itemize
## @item
## @var{x} is a matrix of random samples from the Rayleigh distribution with
## corresponding parameter @var{sigma}. If neither @var{sz} nor @var{r} and
## @var{c} are specified, then @var{x} is of the same size as @{sigma}.
## @end itemize
##
## Examples:
##
## @example
## sigma = 1:6;
## x = raylrnd (sigma)
##
## sz = [2, 3];
## x = raylrnd (0.5, sz)
##
## r = 2;
## c = 3;
## x = raylrnd (0.5, r, c)
## @end example
##
## References:
##
## @enumerate
## @item
## W. L. Martinez and A. R. Martinez. @cite{Computational Statistics
## Handbook with MATLAB.} Chapman & Hall/CRC, pages 547-557, 2001.
##
## @item
## Wikipedia contributors. Rayleigh distribution. @cite{Wikipedia, The Free
## Encyclopedia.}
## @uref{http://en.wikipedia.org/w/index.php?title=Rayleigh_distribution&oldid=69294908},
## August 2006.
## @end enumerate
## @end deftypefn
## Author: Arno Onken <[email protected]>
## Description: Random samples from the Rayleigh distribution
function x = raylrnd (sigma, r, c)
# Check arguments
if (nargin == 1)
sz = size (sigma);
elseif (nargin == 2)
if (! isvector (r) || any ((r < 0) | round (r) != r))
error ("raylrnd: sz must be a vector of non-negative integers")
end
sz = r(:)';
if (! isscalar (sigma) && ! isempty (sigma) && (length (size (sigma)) != length (sz) || any (size (sigma) != sz)))
error ("raylrnd: sigma must be scalar or of size sz");
endif
elseif (nargin == 3)
if (! isscalar (r) || any ((r < 0) | round (r) != r))
error ("raylrnd: r must be a non-negative integer")
end
if (! isscalar (c) || any ((c < 0) | round (c) != c))
error ("raylrnd: c must be a non-negative integer")
end
sz = [r, c];
if (! isscalar (sigma) && ! isempty (sigma) && (length (size (sigma)) != length (sz) || any (size (sigma) != sz)))
error ("raylrnd: sigma must be scalar or of size [r, c]");
endif
else
usage ("x = raylrnd (sigma [, sz |, r, c])");
endif
if (! isempty (sigma) && ! ismatrix (sigma))
error ("raylrnd: sigma must be a numeric matrix");
endif
if (isempty (sigma))
x = [];
elseif (isscalar (sigma) && ! (sigma > 0))
x = NaN .* ones (sz);
else
# Draw random samples
x = sqrt (-2 .* log (1 - rand (sz)) .* sigma .^ 2);
# Continue argument check
k = find (! (sigma > 0));
if (any (k))
x (k) = NaN;
endif
endif
endfunction
%!test
%! sigma = 1:6;
%! x = raylrnd (sigma);
%! assert (size (x), size (sigma));
%! assert (all (x >= 0));
%!test
%! sigma = 0.5;
%! sz = [2, 3];
%! x = raylrnd (sigma, sz);
%! assert (size (x), sz);
%! assert (all (x >= 0));
%!test
%! sigma = 0.5;
%! r = 2;
%! c = 3;
%! x = raylrnd (sigma, r, c);
%! assert (size (x), [r c]);
%! assert (all (x >= 0));
--------------090603030001070304060002
Content-Type: text/plain; charset="us-ascii"
MIME-Version: 1.0
Content-Transfer-Encoding: 7bit
Content-Disposition: inline
_______________________________________________
Octave-sources mailing list
[email protected]
https://www.cae.wisc.edu/mailman/listinfo/octave-sources
--------------090603030001070304060002--