Mercurial > octave-nkf

comparison scripts/statistics/distributions/binoinv.m @ 5410:56e066f5efc1

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 2005-07-13 17:43:35 by jwe]
author jwe
date Wed, 13 Jul 2005 17:43:35 +0000
parents
children bee21f388110
comparison
equal deleted inserted replaced
5409:fda074a55b5c 5410:56e066f5efc1
1 ## Copyright (C) 1995, 1996, 1997 Kurt Hornik
2 ##
3 ## This file is part of Octave.
4 ##
5 ## Octave is free software; you can redistribute it and/or modify it
6 ## under the terms of the GNU General Public License as published by
7 ## the Free Software Foundation; either version 2, or (at your option)
8 ## any later version.
9 ##
10 ## Octave is distributed in the hope that it will be useful, but
11 ## WITHOUT ANY WARRANTY; without even the implied warranty of
12 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
13 ## General Public License for more details.
14 ##
15 ## You should have received a copy of the GNU General Public License
16 ## along with Octave; see the file COPYING. If not, write to the Free
17 ## Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
18 ## 02110-1301, USA.
19
20 ## -*- texinfo -*-
21 ## @deftypefn {Function File} {} binomial_inv (@var{x}, @var{n}, @var{p})
22 ## For each element of @var{x}, compute the quantile at @var{x} of the
23 ## binomial distribution with parameters @var{n} and @var{p}.
24 ## @end deftypefn
25
26 ## Author: KH <Kurt.Hornik@ci.tuwien.ac.at>
27 ## Description: Quantile function of the binomial distribution
28
29 function inv = binomial_inv (x, n, p)
30
31 if (nargin != 3)
32 usage ("binomial_inv (x, n, p)");
33 endif
34
35 if (!isscalar (n) || !isscalar (p))
36 [retval, x, n, p] = common_size (x, n, p);
37 if (retval > 0)
38 error ("binomial_inv: x, n and p must be of common size or scalars");
39 endif
40 endif
41
42 sz = size (x);
43 inv = zeros (sz);
44
45 k = find (!(x >= 0) | !(x <= 1) | !(n >= 0) | (n != round (n))
46 | !(p >= 0) | !(p <= 1));
47 if (any (k))
48 inv(k) = NaN;
49 endif
50
51 k = find ((x >= 0) & (x <= 1) & (n >= 0) & (n == round (n))
52 & (p >= 0) & (p <= 1));
53 if (any (k))
54 if (isscalar (n) && isscalar (p))
55 cdf = binomial_pdf (0, n, p) * ones (size(k));
56 while (any (inv(k) < n))
57 m = find (cdf < x(k));
58 if (any (m))
59 inv(k(m)) = inv(k(m)) + 1;
60 cdf(m) = cdf(m) + binomial_pdf (inv(k(m)), n, p);
61 else
62 break;
63 endif
64 endwhile
65 else
66 cdf = binomial_pdf (0, n(k), p(k));
67 while (any (inv(k) < n(k)))
68 m = find (cdf < x(k));
69 if (any (m))
70 inv(k(m)) = inv(k(m)) + 1;
71 cdf(m) = cdf(m) + binomial_pdf (inv(k(m)), n(k(m)), p(k(m)));
72 else
73 break;
74 endif
75 endwhile
76 endif
77 endif
78
79 endfunction