GNU Octave  4.2.1
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
betainc.cc
Go to the documentation of this file.
1 /*
2 
3 Copyright (C) 1997-2017 John W. Eaton
4 
5 This file is part of Octave.
6 
7 Octave is free software; you can redistribute it and/or modify it
8 under the terms of the GNU General Public License as published by the
9 Free Software Foundation; either version 3 of the License, or (at your
10 option) any later version.
11 
12 Octave is distributed in the hope that it will be useful, but WITHOUT
13 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
14 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
15 for more details.
16 
17 You should have received a copy of the GNU General Public License
18 along with Octave; see the file COPYING. If not, see
19 <http://www.gnu.org/licenses/>.
20 
21 */
22 
23 #if defined (HAVE_CONFIG_H)
24 # include "config.h"
25 #endif
26 
27 #include "lo-specfun.h"
28 
29 #include "defun.h"
30 #include "error.h"
31 #include "errwarn.h"
32 #include "ovl.h"
33 #include "utils.h"
34 
35 
36 DEFUN (betainc, args, ,
37  doc: /* -*- texinfo -*-
38 @deftypefn {} {} betainc (@var{x}, @var{a}, @var{b})
39 Compute the regularized incomplete Beta function.
40 
41 The regularized incomplete Beta function is defined by
42 @tex
43 $$
44  I (x, a, b) = {1 \over {B (a, b)}} \int_0^x t^{(a-z)} (1-t)^{(b-1)} dt.
45 $$
46 @end tex
47 @ifnottex
48 @c Set example in small font to prevent overfull line
49 
50 @smallexample
51 @group
52  x
53  1 /
54 betainc (x, a, b) = ----------- | t^(a-1) (1-t)^(b-1) dt.
55  beta (a, b) /
56  t=0
57 @end group
58 @end smallexample
59 
60 @end ifnottex
61 
62 If @var{x} has more than one component, both @var{a} and @var{b} must be
63 scalars. If @var{x} is a scalar, @var{a} and @var{b} must be of
64 compatible dimensions.
65 @seealso{betaincinv, beta, betaln}
66 @end deftypefn */)
67 {
68  if (args.length () != 3)
69  print_usage ();
70 
72 
73  octave_value x_arg = args(0);
74  octave_value a_arg = args(1);
75  octave_value b_arg = args(2);
76 
77  // FIXME: Can we make a template version of the duplicated code below
78  if (x_arg.is_single_type () || a_arg.is_single_type ()
79  || b_arg.is_single_type ())
80  {
81  if (x_arg.is_scalar_type ())
82  {
83  float x = x_arg.float_value ();
84 
85  if (a_arg.is_scalar_type ())
86  {
87  float a = a_arg.float_value ();
88 
89  if (b_arg.is_scalar_type ())
90  {
91  float b = b_arg.float_value ();
92 
93  retval = octave::math::betainc (x, a, b);
94  }
95  else
96  {
97  Array<float> b = b_arg.float_array_value ();
98 
99  retval = octave::math::betainc (x, a, b);
100  }
101  }
102  else
103  {
104  Array<float> a = a_arg.float_array_value ();
105 
106  if (b_arg.is_scalar_type ())
107  {
108  float b = b_arg.float_value ();
109 
110  retval = octave::math::betainc (x, a, b);
111  }
112  else
113  {
114  Array<float> b = b_arg.float_array_value ();
115 
116  retval = octave::math::betainc (x, a, b);
117  }
118  }
119  }
120  else
121  {
122  Array<float> x = x_arg.float_array_value ();
123 
124  if (a_arg.is_scalar_type ())
125  {
126  float a = a_arg.float_value ();
127 
128  if (b_arg.is_scalar_type ())
129  {
130  float b = b_arg.float_value ();
131 
132  retval = octave::math::betainc (x, a, b);
133  }
134  else
135  {
136  Array<float> b = b_arg.float_array_value ();
137 
138  retval = octave::math::betainc (x, a, b);
139  }
140  }
141  else
142  {
143  Array<float> a = a_arg.float_array_value ();
144 
145  if (b_arg.is_scalar_type ())
146  {
147  float b = b_arg.float_value ();
148 
149  retval = octave::math::betainc (x, a, b);
150  }
151  else
152  {
153  Array<float> b = b_arg.float_array_value ();
154 
155  retval = octave::math::betainc (x, a, b);
156  }
157  }
158  }
159  }
160  else
161  {
162  if (x_arg.is_scalar_type ())
163  {
164  double x = x_arg.double_value ();
165 
166  if (a_arg.is_scalar_type ())
167  {
168  double a = a_arg.double_value ();
169 
170  if (b_arg.is_scalar_type ())
171  {
172  double b = b_arg.double_value ();
173 
174  retval = octave::math::betainc (x, a, b);
175  }
176  else
177  {
178  Array<double> b = b_arg.array_value ();
179 
180  retval = octave::math::betainc (x, a, b);
181  }
182  }
183  else
184  {
185  Array<double> a = a_arg.array_value ();
186 
187  if (b_arg.is_scalar_type ())
188  {
189  double b = b_arg.double_value ();
190 
191  retval = octave::math::betainc (x, a, b);
192  }
193  else
194  {
195  Array<double> b = b_arg.array_value ();
196 
197  retval = octave::math::betainc (x, a, b);
198  }
199  }
200  }
201  else
202  {
203  Array<double> x = x_arg.array_value ();
204 
205  if (a_arg.is_scalar_type ())
206  {
207  double a = a_arg.double_value ();
208 
209  if (b_arg.is_scalar_type ())
210  {
211  double b = b_arg.double_value ();
212 
213  retval = octave::math::betainc (x, a, b);
214  }
215  else
216  {
217  Array<double> b = b_arg.array_value ();
218 
219  retval = octave::math::betainc (x, a, b);
220  }
221  }
222  else
223  {
224  Array<double> a = a_arg.array_value ();
225 
226  if (b_arg.is_scalar_type ())
227  {
228  double b = b_arg.double_value ();
229 
230  retval = octave::math::betainc (x, a, b);
231  }
232  else
233  {
234  Array<double> b = b_arg.array_value ();
235 
236  retval = octave::math::betainc (x, a, b);
237  }
238  }
239  }
240  }
241 
242  return retval;
243 }
244 
245 /*
246 ## Double precision
247 %!test
248 %! a = [1, 1.5, 2, 3];
249 %! b = [4, 3, 2, 1];
250 %! v1 = betainc (1,a,b);
251 %! v2 = [1,1,1,1];
252 %! x = [.2, .4, .6, .8];
253 %! v3 = betainc (x, a, b);
254 %! v4 = 1 - betainc (1.-x, b, a);
255 %! assert (v1, v2, sqrt (eps));
256 %! assert (v3, v4, sqrt (eps));
257 
258 ## Single precision
259 %!test
260 %! a = single ([1, 1.5, 2, 3]);
261 %! b = single ([4, 3, 2, 1]);
262 %! v1 = betainc (1,a,b);
263 %! v2 = single ([1,1,1,1]);
264 %! x = single ([.2, .4, .6, .8]);
265 %! v3 = betainc (x, a, b);
266 %! v4 = 1 - betainc (1.-x, b, a);
267 %! assert (v1, v2, sqrt (eps ("single")));
268 %! assert (v3, v4, sqrt (eps ("single")));
269 
270 ## Mixed double/single precision
271 %!test
272 %! a = single ([1, 1.5, 2, 3]);
273 %! b = [4, 3, 2, 1];
274 %! v1 = betainc (1,a,b);
275 %! v2 = single ([1,1,1,1]);
276 %! x = [.2, .4, .6, .8];
277 %! v3 = betainc (x, a, b);
278 %! v4 = 1-betainc (1.-x, b, a);
279 %! assert (v1, v2, sqrt (eps ("single")));
280 %! assert (v3, v4, sqrt (eps ("single")));
281 
282 %!error betainc ()
283 %!error betainc (1)
284 %!error betainc (1,2)
285 %!error betainc (1,2,3,4)
286 */
287 
288 DEFUN (betaincinv, args, ,
289  doc: /* -*- texinfo -*-
290 @deftypefn {} {} betaincinv (@var{y}, @var{a}, @var{b})
291 Compute the inverse of the incomplete Beta function.
292 
293 The inverse is the value @var{x} such that
294 
295 @example
296 @var{y} == betainc (@var{x}, @var{a}, @var{b})
297 @end example
298 @seealso{betainc, beta, betaln}
299 @end deftypefn */)
300 {
301  if (args.length () != 3)
302  print_usage ();
303 
305 
306  octave_value x_arg = args(0);
307  octave_value a_arg = args(1);
308  octave_value b_arg = args(2);
309 
310  if (x_arg.is_scalar_type ())
311  {
312  double x = x_arg.double_value ();
313 
314  if (a_arg.is_scalar_type ())
315  {
316  double a = a_arg.double_value ();
317 
318  if (b_arg.is_scalar_type ())
319  {
320  double b = b_arg.double_value ();
321 
322  retval = octave::math::betaincinv (x, a, b);
323  }
324  else
325  {
326  Array<double> b = b_arg.array_value ();
327 
328  retval = octave::math::betaincinv (x, a, b);
329  }
330  }
331  else
332  {
333  Array<double> a = a_arg.array_value ();
334 
335  if (b_arg.is_scalar_type ())
336  {
337  double b = b_arg.double_value ();
338 
339  retval = octave::math::betaincinv (x, a, b);
340  }
341  else
342  {
343  Array<double> b = b_arg.array_value ();
344 
345  retval = octave::math::betaincinv (x, a, b);
346  }
347  }
348  }
349  else
350  {
351  Array<double> x = x_arg.array_value ();
352 
353  if (a_arg.is_scalar_type ())
354  {
355  double a = a_arg.double_value ();
356 
357  if (b_arg.is_scalar_type ())
358  {
359  double b = b_arg.double_value ();
360 
361  retval = octave::math::betaincinv (x, a, b);
362  }
363  else
364  {
365  Array<double> b = b_arg.array_value ();
366 
367  retval = octave::math::betaincinv (x, a, b);
368  }
369  }
370  else
371  {
372  Array<double> a = a_arg.array_value ();
373 
374  if (b_arg.is_scalar_type ())
375  {
376  double b = b_arg.double_value ();
377 
378  retval = octave::math::betaincinv (x, a, b);
379  }
380  else
381  {
382  Array<double> b = b_arg.array_value ();
383 
384  retval = octave::math::betaincinv (x, a, b);
385  }
386  }
387  }
388 
389  // FIXME: It would be better to have an algorithm for betaincinv which
390  // accepted float inputs and returned float outputs. As it is, we do
391  // extra work to calculate betaincinv to double precision and then throw
392  // that precision away.
393  if (x_arg.is_single_type () || a_arg.is_single_type ()
394  || b_arg.is_single_type ())
395  {
396  retval = Array<float> (retval.array_value ());
397  }
398 
399  return retval;
400 }
401 
402 /*
403 %!assert (betaincinv ([0.875 0.6875], [1 2], 3), [0.5 0.5], sqrt (eps))
404 %!assert (betaincinv (0.5, 3, 3), 0.5, sqrt (eps))
405 %!assert (betaincinv (0.34375, 4, 3), 0.5, sqrt (eps))
406 %!assert (betaincinv (0.2265625, 5, 3), 0.5, sqrt (eps))
407 %!assert (betaincinv (0.14453125, 6, 3), 0.5, sqrt (eps))
408 %!assert (betaincinv (0.08984375, 7, 3), 0.5, sqrt (eps))
409 %!assert (betaincinv (0.0546875, 8, 3), 0.5, sqrt (eps))
410 %!assert (betaincinv (0.03271484375, 9, 3), 0.5, sqrt (eps))
411 %!assert (betaincinv (0.019287109375, 10, 3), 0.5, sqrt (eps))
412 
413 ## Test class single as well
414 %!assert (betaincinv ([0.875 0.6875], [1 2], single (3)), [0.5 0.5], sqrt (eps ("single")))
415 %!assert (betaincinv (0.5, 3, single (3)), 0.5, sqrt (eps ("single")))
416 %!assert (betaincinv (0.34375, 4, single (3)), 0.5, sqrt (eps ("single")))
417 
418 ## Extreme values
419 %!assert (betaincinv (0, 42, 42), 0, sqrt (eps))
420 %!assert (betaincinv (1, 42, 42), 1, sqrt (eps))
421 
422 %!error betaincinv ()
423 %!error betaincinv (1, 2)
424 */
OCTINTERP_API void print_usage(void)
Definition: defun.cc:52
double betaincinv(double y, double p, double q)
Definition: lo-specfun.cc:3193
bool is_scalar_type(void) const
Definition: ov.h:673
#define DEFUN(name, args_name, nargout_name, doc)
Definition: defun.h:46
float float_value(bool frc_str_conv=false) const
Definition: ov.h:778
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:398
JNIEnv void * args
Definition: ov-java.cc:67
FloatNDArray float_array_value(bool frc_str_conv=false) const
Definition: ov.h:796
octave_value retval
Definition: data.cc:6294
NDArray array_value(bool frc_str_conv=false) const
Definition: ov.h:793
double betainc(double x, double a, double b)
Definition: lo-specfun.cc:2265
b
Definition: cellfun.cc:398
bool is_single_type(void) const
Definition: ov.h:627
double double_value(bool frc_str_conv=false) const
Definition: ov.h:775
F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &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 F77_REAL &F77_RET_T F77_DBLE &F77_RET_T F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T const F77_DBLE F77_DBLE &F77_RET_T const F77_REAL F77_REAL &F77_RET_T F77_REAL F77_REAL &F77_RET_T F77_DBLE F77_DBLE &F77_RET_T const F77_DBLE * x