Mercurial > octave-nkf

comparison scripts/statistics/distributions/normpdf.m @ 13171:19b9f17d22af

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Overhaul of statistical distribution functions Support class "single" 75% reduction in memory usage More Matlab compatibility for corner cases * betacdf.m, betainv.m, betapdf.m, betarnd.m, binocdf.m, binoinv.m, binopdf.m, binornd.m, cauchy_cdf.m, cauchy_inv.m, cauchy_pdf.m, cauchy_rnd.m, chi2cdf.m, chi2inv.m, chi2pdf.m, chi2rnd.m, discrete_cdf.m, discrete_inv.m, discrete_pdf.m, discrete_rnd.m, empirical_cdf.m, empirical_inv.m, empirical_pdf.m, empirical_rnd.m, expcdf.m, expinv.m, exppdf.m, exprnd.m, fcdf.m, finv.m, fpdf.m, frnd.m, gamcdf.m, gaminv.m, gampdf.m, gamrnd.m, geocdf.m, geoinv.m, geopdf.m, geornd.m, hygecdf.m, hygeinv.m, hygepdf.m, hygernd.m, kolmogorov_smirnov_cdf.m, laplace_cdf.m, laplace_inv.m, laplace_pdf.m, laplace_rnd.m, logistic_cdf.m, logistic_inv.m, logistic_pdf.m, logistic_rnd.m, logncdf.m, logninv.m, lognpdf.m, lognrnd.m, nbincdf.m, nbininv.m, nbinpdf.m, nbinrnd.m, normcdf.m, norminv.m, normpdf.m, normrnd.m, poisscdf.m, poissinv.m, poisspdf.m, poissrnd.m, stdnormal_cdf.m, stdnormal_inv.m, stdnormal_pdf.m, stdnormal_rnd.m, tcdf.m, tinv.m, tpdf.m, trnd.m, unidcdf.m, unidinv.m, unidpdf.m, unidrnd.m, unifcdf.m, unifinv.m, unifpdf.m, unifrnd.m, wblcdf.m, wblinv.m, wblpdf.m, wblrnd.m: Return "single" outputs for "single" inputs, Use logical indexing rather than find() for 75% memory savings, Add tests for all functions, Use consistent documentation across all functions, More Matlab compatibilitcy for corner cases.
author Rik <octave@nomad.inbox5.com>
date Tue, 20 Sep 2011 12:13:13 -0700
parents e7cc2d4a6db3
children 72c96de7a403
comparison
equal deleted inserted replaced
13169:078729410a0d 13171:19b9f17d22af
1 ## Copyright (C) 2011 Rik Wehbring
1 ## Copyright (C) 1995-2011 Kurt Hornik 2 ## Copyright (C) 1995-2011 Kurt Hornik
2 ## 3 ##
3 ## This file is part of Octave. 4 ## This file is part of Octave.
4 ## 5 ##
5 ## Octave is free software; you can redistribute it and/or modify it 6 ## Octave is free software; you can redistribute it and/or modify it
15 ## You should have received a copy of the GNU General Public License 16 ## You should have received a copy of the GNU General Public License
16 ## along with Octave; see the file COPYING. If not, see 17 ## along with Octave; see the file COPYING. If not, see
17 ## <http://www.gnu.org/licenses/>. 18 ## <http://www.gnu.org/licenses/>.
18 19
19 ## -*- texinfo -*- 20 ## -*- texinfo -*-
20 ## @deftypefn {Function File} {} normpdf (@var{x}, @var{m}, @var{s}) 21 ## @deftypefn {Function File} {} normpdf (@var{x})
22 ## @deftypefnx {Function File} {} normpdf (@var{x}, @var{mu}, @var{sigma})
21 ## For each element of @var{x}, compute the probability density function 23 ## For each element of @var{x}, compute the probability density function
22 ## (PDF) at @var{x} of the normal distribution with mean @var{m} and 24 ## (PDF) at @var{x} of the normal distribution with mean @var{mu} and
23 ## standard deviation @var{s}. 25 ## standard deviation @var{sigma}.
24 ## 26 ##
25 ## Default values are @var{m} = 0, @var{s} = 1. 27 ## Default values are @var{mu} = 0, @var{sigma} = 1.
26 ## @end deftypefn 28 ## @end deftypefn
27 29
28 ## Author: TT <Teresa.Twaroch@ci.tuwien.ac.at> 30 ## Author: TT <Teresa.Twaroch@ci.tuwien.ac.at>
29 ## Description: PDF of the normal distribution 31 ## Description: PDF of the normal distribution
30 32
31 function pdf = normpdf (x, m, s) 33 function pdf = normpdf (x, mu = 0, sigma = 1)
32 34
33 if (nargin != 1 && nargin != 3) 35 if (nargin != 1 && nargin != 3)
34 print_usage (); 36 print_usage ();
35 endif 37 endif
36 38
37 if (nargin == 1) 39 if (!isscalar (mu) || !isscalar (sigma))
38 m = 0; 40 [retval, x, mu, sigma] = common_size (x, mu, sigma);
39 s = 1;
40 endif
41
42 if (!isscalar (m) || !isscalar (s))
43 [retval, x, m, s] = common_size (x, m, s);
44 if (retval > 0) 41 if (retval > 0)
45 error ("normpdf: X, M and S must be of common size or scalars"); 42 error ("normpdf: X, MU, and SIGMA must be of common size or scalars");
46 endif 43 endif
47 endif 44 endif
48 45
49 sz = size (x); 46 if (iscomplex (x) || iscomplex (mu) || iscomplex (sigma))
50 pdf = zeros (sz); 47 error ("normpdf: X, MU, and SIGMA must not be complex");
48 endif
51 49
52 if (isscalar (m) && isscalar (s)) 50 if (isa (x, "single") || isa (mu, "single") || isa (sigma, "single"))
53 if (find (isinf (m) | isnan (m) | !(s > 0) | !(s < Inf))) 51 pdf = zeros (size (x), "single");
54 pdf = NaN (sz); 52 else
53 pdf = zeros (size (x));
54 endif
55
56 if (isscalar (mu) && isscalar (sigma))
57 if (!isinf (mu) && !isnan (mu) && (sigma > 0) && (sigma < Inf))
58 pdf = stdnormal_pdf ((x - mu) / sigma) / sigma;
55 else 59 else
56 pdf = stdnormal_pdf ((x - m) ./ s) ./ s; 60 pdf = NaN (size (x), class (pdf));
57 endif 61 endif
58 else 62 else
59 k = find (isinf (m) | isnan (m) | !(s > 0) | !(s < Inf)); 63 k = isinf (mu) | !(sigma > 0) | !(sigma < Inf);
60 if (any (k)) 64 pdf(k) = NaN;
61 pdf(k) = NaN;
62 endif
63 65
64 k = find (!isinf (m) & !isnan (m) & (s > 0) & (s < Inf)); 66 k = !isinf (mu) & (sigma > 0) & (sigma < Inf);
65 if (any (k)) 67 pdf(k) = stdnormal_pdf ((x(k) - mu(k)) ./ sigma(k)) ./ sigma(k);
66 pdf(k) = stdnormal_pdf ((x(k) - m(k)) ./ s(k)) ./ s(k);
67 endif
68 endif 68 endif
69 69
70 pdf((s == 0) & (x == m)) = Inf; 70 endfunction
71 pdf((s == 0) & ((x < m) | (x > m))) = 0;
72 71
73 endfunction 72
73 %!shared x,y
74 %! x = [-Inf 1 2 Inf];
75 %! y = 1/sqrt(2*pi)*exp (-(x-1).^2/2);
76 %!assert(normpdf (x, ones(1,4), ones(1,4)), y);
77 %!assert(normpdf (x, 1, ones(1,4)), y);
78 %!assert(normpdf (x, ones(1,4), 1), y);
79 %!assert(normpdf (x, [0 -Inf NaN Inf], 1), [y(1) NaN NaN NaN]);
80 %!assert(normpdf (x, 1, [Inf NaN -1 0]), [NaN NaN NaN NaN]);
81 %!assert(normpdf ([x, NaN], 1, 1), [y, NaN]);
82
83 %% Test class of input preserved
84 %!assert(normpdf (single([x, NaN]), 1, 1), single([y, NaN]), eps("single"));
85 %!assert(normpdf ([x, NaN], single(1), 1), single([y, NaN]), eps("single"));
86 %!assert(normpdf ([x, NaN], 1, single(1)), single([y, NaN]), eps("single"));
87
88 %% Test input validation
89 %!error normpdf ()
90 %!error normpdf (1,2)
91 %!error normpdf (1,2,3,4)
92 %!error normpdf (ones(3),ones(2),ones(2))
93 %!error normpdf (ones(2),ones(3),ones(2))
94 %!error normpdf (ones(2),ones(2),ones(3))
95 %!error normpdf (i, 2, 2)
96 %!error normpdf (2, i, 2)
97 %!error normpdf (2, 2, i)
98