Mercurial > octave-nkf
comparison scripts/statistics/distributions/gaminv.m @ 13171:19b9f17d22af
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 | fd0a3ac60b0e |
children | 0c15fece33ad |
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 |
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} {} gaminv (@var{x}, @var{a}, @var{b}) | 21 ## @deftypefn {Function File} {} gaminv (@var{x}, @var{a}, @var{b}) |
21 ## For each component of @var{x}, compute the quantile (the inverse of | 22 ## For each element of @var{x}, compute the quantile (the inverse of |
22 ## the CDF) at @var{x} of the Gamma distribution with parameters @var{a} | 23 ## the CDF) at @var{x} of the Gamma distribution with parameters @var{a} |
23 ## and @var{b}. | 24 ## and @var{b}. |
24 ## @seealso{gamma, gammaln, gammainc, gampdf, gamcdf, gamrnd} | |
25 ## @end deftypefn | 25 ## @end deftypefn |
26 | 26 |
27 ## Author: KH <Kurt.Hornik@wu-wien.ac.at> | 27 ## Author: KH <Kurt.Hornik@wu-wien.ac.at> |
28 ## Description: Quantile function of the Gamma distribution | 28 ## Description: Quantile function of the Gamma distribution |
29 | 29 |
31 | 31 |
32 if (nargin != 3) | 32 if (nargin != 3) |
33 print_usage (); | 33 print_usage (); |
34 endif | 34 endif |
35 | 35 |
36 if (!isscalar (a) || !isscalar(b)) | 36 if (!isscalar (a) || !isscalar (b)) |
37 [retval, x, a, b] = common_size (x, a, b); | 37 [retval, x, a, b] = common_size (x, a, b); |
38 if (retval > 0) | 38 if (retval > 0) |
39 error ("gaminv: X, A and B must be of common size or scalars"); | 39 error ("gaminv: X, A, and B must be of common size or scalars"); |
40 endif | 40 endif |
41 endif | 41 endif |
42 | 42 |
43 sz = size (x); | 43 if (iscomplex (x) || iscomplex (a) || iscomplex (b)) |
44 inv = zeros (sz); | 44 error ("gaminv: X, A, and B must not be complex"); |
45 | |
46 k = find ((x < 0) | (x > 1) | isnan (x) | !(a > 0) | !(b > 0)); | |
47 if (any (k)) | |
48 inv (k) = NaN; | |
49 endif | 45 endif |
50 | 46 |
51 k = find ((x == 1) & (a > 0) & (b > 0)); | 47 if (isa (x, "single") || isa (a, "single") || isa (b, "single")) |
52 if (any (k)) | 48 inv = zeros (size (x), "single"); |
53 inv (k) = Inf; | 49 else |
50 inv = zeros (size (x)); | |
54 endif | 51 endif |
55 | 52 |
56 k = find ((x > 0) & (x < 1) & (a > 0) & (b > 0)); | 53 k = ((x < 0) | (x > 1) | isnan (x) |
54 | !(a > 0) | !(a < Inf) | !(b > 0) | !(b < Inf)); | |
55 inv(k) = NaN; | |
56 | |
57 k = (x == 1) & (a > 0) & (a < Inf) & (b > 0) & (b < Inf); | |
58 inv(k) = Inf; | |
59 | |
60 k = find ((x > 0) & (x < 1) & (a > 0) & (a < Inf) & (b > 0) & (b < Inf)); | |
57 if (any (k)) | 61 if (any (k)) |
58 if (!isscalar(a) || !isscalar(b)) | 62 if (!isscalar (a) || !isscalar (b)) |
59 a = a (k); | 63 a = a(k); |
60 b = b (k); | 64 b = b(k); |
61 y = a .* b; | 65 y = a .* b; |
62 else | 66 else |
63 y = a * b * ones (size (k)); | 67 y = a * b * ones (size (k)); |
64 endif | 68 endif |
65 x = x (k); | 69 x = x(k); |
66 | 70 |
67 if (isa (x, "single")) | 71 if (isa (x, "single")) |
68 myeps = eps ("single"); | 72 myeps = eps ("single"); |
69 else | 73 else |
70 myeps = eps; | 74 myeps = eps; |
88 break; | 92 break; |
89 endif | 93 endif |
90 y_old = y_new; | 94 y_old = y_new; |
91 endfor | 95 endfor |
92 | 96 |
93 inv (k) = y_new; | 97 inv(k) = y_new; |
94 endif | 98 endif |
95 | 99 |
96 endfunction | 100 endfunction |
101 | |
102 | |
103 %!shared x | |
104 %! x = [-1 0 0.63212055882855778 1 2]; | |
105 %!assert(gaminv (x, ones(1,5), ones(1,5)), [NaN 0 1 Inf NaN], eps); | |
106 %!assert(gaminv (x, 1, ones(1,5)), [NaN 0 1 Inf NaN], eps); | |
107 %!assert(gaminv (x, ones(1,5), 1), [NaN 0 1 Inf NaN], eps); | |
108 %!assert(gaminv (x, [1 -Inf NaN Inf 1], 1), [NaN NaN NaN NaN NaN]); | |
109 %!assert(gaminv (x, 1, [1 -Inf NaN Inf 1]), [NaN NaN NaN NaN NaN]); | |
110 %!assert(gaminv ([x(1:2) NaN x(4:5)], 1, 1), [NaN 0 NaN Inf NaN]); | |
111 | |
112 %% Test class of input preserved | |
113 %!assert(gaminv ([x, NaN], 1, 1), [NaN 0 1 Inf NaN NaN], eps); | |
114 %!assert(gaminv (single([x, NaN]), 1, 1), single([NaN 0 1 Inf NaN NaN]), eps("single")); | |
115 %!assert(gaminv ([x, NaN], single(1), 1), single([NaN 0 1 Inf NaN NaN]), eps("single")); | |
116 %!assert(gaminv ([x, NaN], 1, single(1)), single([NaN 0 1 Inf NaN NaN]), eps("single")); | |
117 | |
118 %% Test input validation | |
119 %!error gaminv () | |
120 %!error gaminv (1) | |
121 %!error gaminv (1,2) | |
122 %!error gaminv (1,2,3,4) | |
123 %!error gaminv (ones(3),ones(2),ones(2)) | |
124 %!error gaminv (ones(2),ones(3),ones(2)) | |
125 %!error gaminv (ones(2),ones(2),ones(3)) | |
126 %!error gaminv (i, 2, 2) | |
127 %!error gaminv (2, i, 2) | |
128 %!error gaminv (2, 2, i) | |
129 |