GNU Octave  4.4.1
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
__betainc__.cc
Go to the documentation of this file.
1 /*
2 
3 Copyright (C) 2018 Michele Ginesi
4 Copyright (C) 2018 Stefan Schlögl
5 
6 This file is part of Octave.
7 
8 Octave is free software: you can redistribute it and/or modify it
9 under the terms of the GNU General Public License as published by
10 the Free Software Foundation, either version 3 of the License, or
11 (at your option) any later version.
12 
13 Octave is distributed in the hope that it will be useful, but
14 WITHOUT ANY WARRANTY; without even the implied warranty of
15 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16 GNU General Public License for more details.
17 
18 You should have received a copy of the GNU General Public License
19 along with Octave; see the file COPYING. If not, see
20 <https://www.gnu.org/licenses/>.
21 
22 */
23 
24 #if defined (HAVE_CONFIG_H)
25 # include "config.h"
26 #endif
27 
28 #include "defun.h"
29 #include "dNDArray.h"
30 #include "fNDArray.h"
31 
32 DEFUN (__betainc__, args, ,
33  doc: /* -*- texinfo -*-
34 @deftypefn {} {@var{y} =} __betainc__ (@var{x}, @var{a}, @var{b})
35 Continued fraction for incomplete beta function.
36 @end deftypefn */)
37 {
38  int nargin = args.length ();
39 
40  if (nargin != 3)
41  print_usage ();
42 
43  bool is_single = (args(0).is_single_type () || args(1).is_single_type ()
44  || args(2).is_single_type ());
45 
46  // Total number of scenarios: get maximum of length of all vectors
47  int numel_x = args(0).numel ();
48  int numel_a = args(1).numel ();
49  int numel_b = args(2).numel ();
50  int len = std::max (std::max (numel_x, numel_a), numel_b);
51 
53  // Initialize output dimension vector
54  dim_vector output_dv (len, 1);
55 
56  // Lentz's algorithm in two cases: single and double precision
57  if (is_single)
58  {
59  // Initialize output and inputs
60  FloatColumnVector output (output_dv);
61  FloatNDArray x, a, b;
62 
63  if (numel_x == 1)
64  x = FloatNDArray (output_dv, args(0).float_scalar_value ());
65  else
66  x = args(0).float_array_value ();
67 
68 
69  if (numel_a == 1)
70  a = FloatNDArray (output_dv, args(1).float_scalar_value ());
71  else
72  a = args(1).float_array_value ();
73 
74  if (numel_b == 1)
75  b = FloatNDArray (output_dv, args(2).float_scalar_value ());
76  else
77  b = args(2).float_array_value ();
78 
79  // Initialize variables used in algorithm
80  static const float tiny = octave::math::exp2 (-50.0f);
81  static const float eps = std::numeric_limits<float>::epsilon ();
82  float xj, x2, y, Cj, Dj, aj, bj, Deltaj, alpha_j, beta_j;
83  int j, maxit;
84  maxit = 200;
85 
86  // Loop over all elements
87  for (octave_idx_type i = 0; i < len; ++i)
88  {
89  // Catch Ctrl+C
91 
92  // Variable initialization for the current element
93  xj = x(i);
94  y = tiny;
95  Cj = y;
96  Dj = 0;
97  aj = a(i);
98  bj = b(i);
99  Deltaj = 0;
100  alpha_j = 1;
101  beta_j = aj - (aj * (aj + bj)) / (aj + 1) * xj;
102  x2 = xj * xj;
103  j = 1;
104 
105  // Lentz's algorithm
106  while ((std::abs ((Deltaj - 1)) > eps) && (j < maxit))
107  {
108  Dj = beta_j + alpha_j * Dj;
109  if (Dj == 0)
110  Dj = tiny;
111  Cj = beta_j + alpha_j / Cj;
112  if (Cj == 0)
113  Cj = tiny;
114  Dj = 1 / Dj;
115  Deltaj = Cj * Dj;
116  y *= Deltaj;
117  alpha_j = ((aj + j - 1) * (aj + bj + j -1) * (bj - j) * j)
118  / ((aj + 2 * j - 1) * (aj + 2 * j - 1)) * x2;
119  beta_j = aj + 2 * j + ((j * (bj - j)) / (aj + 2 * j - 1)
120  - ((aj + j) * (aj + bj + j)) / (aj + 2 * j + 1)) * xj;
121  j++;
122  }
123 
124  output(i) = y;
125  }
126 
127  retval(0) = output;
128  }
129  else
130  {
131  // Initialize output and inputs
132  ColumnVector output (output_dv);
133  NDArray x, a, b;
134 
135  if (numel_x == 1)
136  x = NDArray (output_dv, args(0).scalar_value ());
137  else
138  x = args(0).array_value ();
139 
140  if (numel_a == 1)
141  a = NDArray (output_dv, args(1).scalar_value ());
142  else
143  a = args(1).array_value ();
144 
145  if (numel_b == 1)
146  b = NDArray (output_dv, args(2).scalar_value ());
147  else
148  b = args(2).array_value ();
149 
150  // Initialize variables used in algorithm
151  static const double tiny = octave::math::exp2 (-100.0);
152  static const double eps = std::numeric_limits<double>::epsilon ();
153  double xj, x2, y, Cj, Dj, aj, bj, Deltaj, alpha_j, beta_j;
154  int j, maxit;
155  maxit = 200;
156 
157  // Loop over all elements
158  for (octave_idx_type i = 0; i < len; ++i)
159  {
160  // Catch Ctrl+C
161  OCTAVE_QUIT;
162 
163  // Variable initialization for the current element
164  xj = x(i);
165  y = tiny;
166  Cj = y;
167  Dj = 0;
168  aj = a(i);
169  bj = b(i);
170  Deltaj = 0;
171  alpha_j = 1;
172  beta_j = aj - (aj * (aj + bj)) / (aj + 1) * xj;
173  x2 = xj * xj;
174  j = 1;
175 
176  // Lentz's algorithm
177  while ((std::abs ((Deltaj - 1)) > eps) && (j < maxit))
178  {
179  Dj = beta_j + alpha_j * Dj;
180  if (Dj == 0)
181  Dj = tiny;
182  Cj = beta_j + alpha_j / Cj;
183  if (Cj == 0)
184  Cj = tiny;
185  Dj = 1 / Dj;
186  Deltaj = Cj * Dj;
187  y *= Deltaj;
188  alpha_j = ((aj + j - 1) * (aj + bj + j - 1) * (bj - j) * j)
189  / ((aj + 2 * j - 1) * (aj + 2 * j - 1)) * x2;
190  beta_j = aj + 2 * j + ((j * (bj - j)) / (aj + 2 * j - 1)
191  - ((aj + j) * (aj + bj + j)) / (aj + 2 * j + 1)) * xj;
192  j++;
193  }
194 
195  output(i) = y;
196  }
197 
198  retval(0) = output;
199  }
200 
201  return retval;
202 }
OCTINTERP_API void print_usage(void)
Definition: defun.cc:54
F77_RET_T const F77_REAL const F77_REAL F77_REAL &F77_RET_T const F77_DBLE const F77_DBLE F77_DBLE &F77_RET_T const F77_DBLE F77_DBLE &F77_RET_T const F77_REAL F77_REAL &F77_RET_T const F77_DBLE const F77_DBLE * f
#define DEFUN(name, args_name, nargout_name, doc)
Macro to define a builtin function.
Definition: defun.h:53
static T abs(T x)
Definition: pr-output.cc:1696
calling an anonymous function involves an overhead quite comparable to the overhead of an m file function Passing a handle to a built in function is because the interpreter is not involved in the internal loop For a
Definition: cellfun.cc:400
octave_value retval
Definition: data.cc:6246
#define eps(C)
double exp2(double x)
Definition: lo-mappers.h:107
charNDArray max(char d, const charNDArray &m)
Definition: chNDArray.cc:227
the element is set to zero In other the statement xample y
Definition: data.cc:5264
b
Definition: cellfun.cc:400
args.length() nargin
Definition: file-io.cc:589
for i
Definition: data.cc:5264
#define OCTAVE_QUIT
Definition: quit.h:193
Vector representing the dimensions (size) of an Array.
Definition: dim-vector.h:87
F77_RET_T const F77_REAL const F77_REAL F77_REAL &F77_RET_T const F77_DBLE const F77_DBLE F77_DBLE &F77_RET_T const F77_DBLE F77_DBLE &F77_RET_T const F77_REAL F77_REAL &F77_RET_T const F77_DBLE * x