← Back to branch summary

~gabriel1984sibiu/octave/octave

« back to all changes in this revision

Viewing changes to scripts/statistics/distributions/binoinv.m

  • Committer: Grevutiu Gabriel
  • Date: 2014-01-02 13:05:54 UTC
  • Revision ID: gabriel1984sibiu@gmail.com-20140102130554-3r7ivdjln1ni6kcg
New version (3.8.0) from upstream.

Show diffs side-by-side

added added

removed removed

Lines of Context:
scripts/statistics/distributions/binoinv.m
 
1
## Copyright (C) 2012 Rik Wehbring
 
2
## Copyright (C) 1995-2013 Kurt Hornik
 
3
##
 
4
## This file is part of Octave.
 
5
##
 
6
## Octave is free software; you can redistribute it and/or modify it
 
7
## under the terms of the GNU General Public License as published by
 
8
## the Free Software Foundation; either version 3 of the License, or (at
 
9
## your option) any later version.
 
10
##
 
11
## Octave is distributed in the hope that it will be useful, but
 
12
## WITHOUT ANY WARRANTY; without even the implied warranty of
 
13
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 
14
## General Public License for more details.
 
15
##
 
16
## You should have received a copy of the GNU General Public License
 
17
## along with Octave; see the file COPYING.  If not, see
 
18
## <http://www.gnu.org/licenses/>.
 
19
 
 
20
## -*- texinfo -*-
 
21
## @deftypefn {Function File} {} binoinv (@var{x}, @var{n}, @var{p})
 
22
## For each element of @var{x}, compute the quantile (the inverse of
 
23
## the CDF) at @var{x} of the binomial distribution with parameters 
 
24
## @var{n} and @var{p}, where @var{n} is the number of trials and
 
25
## @var{p} is the probability of success.
 
26
## @end deftypefn
 
27
 
 
28
## Author: KH <Kurt.Hornik@wu-wien.ac.at>
 
29
## Description: Quantile function of the binomial distribution
 
30
 
 
31
function inv = binoinv (x, n, p)
 
32
 
 
33
  if (nargin != 3)
 
34
    print_usage ();
 
35
  endif
 
36
 
 
37
  if (!isscalar (n) || !isscalar (p))
 
38
    [retval, x, n, p] = common_size (x, n, p);
 
39
    if (retval > 0)
 
40
      error ("binoinv: X, N, and P must be of common size or scalars");
 
41
    endif
 
42
  endif
 
43
 
 
44
  if (iscomplex (x) || iscomplex (n) || iscomplex (p))
 
45
    error ("binoinv: X, N, and P must not be complex");
 
46
  endif
 
47
 
 
48
  if (isa (x, "single") || isa (n, "single") || isa (p, "single"));
 
49
    inv = zeros (size (x), "single");
 
50
  else
 
51
    inv = zeros (size (x));
 
52
  endif
 
53
 
 
54
  k = (!(x >= 0) | !(x <= 1) | !(n >= 0) | (n != fix (n)) |
 
55
       !(p >= 0) | !(p <= 1));
 
56
  inv(k) = NaN;
 
57
 
 
58
  k = find ((x >= 0) & (x <= 1) & (n >= 0) & (n == fix (n)
 
59
             & (p >= 0) & (p <= 1)));
 
60
  if (any (k))
 
61
    if (isscalar (n) && isscalar (p))
 
62
      cdf = binopdf (0, n, p) * ones (size (k));
 
63
      while (any (inv(k) < n))
 
64
        m = find (cdf < x(k));
 
65
        if (any (m))
 
66
          inv(k(m)) = inv(k(m)) + 1;
 
67
          cdf(m) = cdf(m) + binopdf (inv(k(m)), n, p);
 
68
        else
 
69
          break;
 
70
        endif
 
71
      endwhile
 
72
    else
 
73
      cdf = binopdf (0, n(k), p(k));
 
74
      while (any (inv(k) < n(k)))
 
75
        m = find (cdf < x(k));
 
76
        if (any (m))
 
77
          inv(k(m)) = inv(k(m)) + 1;
 
78
          cdf(m) = cdf(m) + binopdf (inv(k(m)), n(k(m)), p(k(m)));
 
79
        else
 
80
          break;
 
81
        endif
 
82
      endwhile
 
83
    endif
 
84
  endif
 
85
 
 
86
endfunction
 
87
 
 
88
 
 
89
%!shared x
 
90
%! x = [-1 0 0.5 1 2];
 
91
%!assert (binoinv (x, 2*ones (1,5), 0.5*ones (1,5)), [NaN 0 1 2 NaN])
 
92
%!assert (binoinv (x, 2, 0.5*ones (1,5)), [NaN 0 1 2 NaN])
 
93
%!assert (binoinv (x, 2*ones (1,5), 0.5), [NaN 0 1 2 NaN])
 
94
%!assert (binoinv (x, 2*[0 -1 NaN 1.1 1], 0.5), [NaN NaN NaN NaN NaN])
 
95
%!assert (binoinv (x, 2, 0.5*[0 -1 NaN 3 1]), [NaN NaN NaN NaN NaN])
 
96
%!assert (binoinv ([x(1:2) NaN x(4:5)], 2, 0.5), [NaN 0 NaN 2 NaN])
 
97
 
 
98
%% Test class of input preserved
 
99
%!assert (binoinv ([x, NaN], 2, 0.5), [NaN 0 1 2 NaN NaN])
 
100
%!assert (binoinv (single ([x, NaN]), 2, 0.5), single ([NaN 0 1 2 NaN NaN]))
 
101
%!assert (binoinv ([x, NaN], single (2), 0.5), single ([NaN 0 1 2 NaN NaN]))
 
102
%!assert (binoinv ([x, NaN], 2, single (0.5)), single ([NaN 0 1 2 NaN NaN]))
 
103
 
 
104
%% Test input validation
 
105
%!error binoinv ()
 
106
%!error binoinv (1)
 
107
%!error binoinv (1,2)
 
108
%!error binoinv (1,2,3,4)
 
109
%!error binoinv (ones (3), ones (2), ones (2))
 
110
%!error binoinv (ones (2), ones (3), ones (2))
 
111
%!error binoinv (ones (2), ones (2), ones (3))
 
112
%!error binoinv (i, 2, 2)
 
113
%!error binoinv (2, i, 2)
 
114
%!error binoinv (2, 2, i)
 
115