GNU Octave  4.0.0
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-2015 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 #ifdef 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 "gripes.h"
32 #include "oct-obj.h"
33 #include "utils.h"
34 
35 
36 DEFUN (betainc, args, ,
37  "-*- texinfo -*-\n\
38 @deftypefn {Mapping Function} {} betainc (@var{x}, @var{a}, @var{b})\n\
39 Compute the regularized incomplete Beta function.\n\
40 \n\
41 The regularized incomplete Beta function is defined by\n\
42 @tex\n\
43 $$\n\
44  I (x, a, b) = {1 \\over {B (a, b)}} \\int_0^x t^{(a-z)} (1-t)^{(b-1)} dt.\n\
45 $$\n\
46 @end tex\n\
47 @ifnottex\n\
48 @c Set example in small font to prevent overfull line\n\
49 \n\
50 @smallexample\n\
51 @group\n\
52  x\n\
53  1 /\n\
54 betainc (x, a, b) = ----------- | t^(a-1) (1-t)^(b-1) dt.\n\
55  beta (a, b) /\n\
56  t=0\n\
57 @end group\n\
58 @end smallexample\n\
59 \n\
60 @end ifnottex\n\
61 \n\
62 If @var{x} has more than one component, both @var{a} and @var{b} must be\n\
63 scalars. If @var{x} is a scalar, @var{a} and @var{b} must be of\n\
64 compatible dimensions.\n\
65 @seealso{betaincinv, beta, betaln}\n\
66 @end deftypefn")
67 {
68  octave_value retval;
69 
70  int nargin = args.length ();
71 
72  if (nargin == 3)
73  {
74  octave_value x_arg = args(0);
75  octave_value a_arg = args(1);
76  octave_value b_arg = args(2);
77 
78  // FIXME: Can we make a template version of the duplicated code below
79  if (x_arg.is_single_type () || a_arg.is_single_type ()
80  || b_arg.is_single_type ())
81  {
82  if (x_arg.is_scalar_type ())
83  {
84  float x = x_arg.float_value ();
85 
86  if (a_arg.is_scalar_type ())
87  {
88  float a = a_arg.float_value ();
89 
90  if (! error_state)
91  {
92  if (b_arg.is_scalar_type ())
93  {
94  float b = b_arg.float_value ();
95 
96  if (! error_state)
97  retval = betainc (x, a, b);
98  }
99  else
100  {
101  Array<float> b = b_arg.float_array_value ();
102 
103  if (! error_state)
104  retval = betainc (x, a, b);
105  }
106  }
107  }
108  else
109  {
110  Array<float> a = a_arg.float_array_value ();
111 
112  if (! error_state)
113  {
114  if (b_arg.is_scalar_type ())
115  {
116  float b = b_arg.float_value ();
117 
118  if (! error_state)
119  retval = betainc (x, a, b);
120  }
121  else
122  {
123  Array<float> b = b_arg.float_array_value ();
124 
125  if (! error_state)
126  retval = betainc (x, a, b);
127  }
128  }
129  }
130  }
131  else
132  {
133  Array<float> x = x_arg.float_array_value ();
134 
135  if (a_arg.is_scalar_type ())
136  {
137  float a = a_arg.float_value ();
138 
139  if (! error_state)
140  {
141  if (b_arg.is_scalar_type ())
142  {
143  float b = b_arg.float_value ();
144 
145  if (! error_state)
146  retval = betainc (x, a, b);
147  }
148  else
149  {
150  Array<float> b = b_arg.float_array_value ();
151 
152  if (! error_state)
153  retval = betainc (x, a, b);
154  }
155  }
156  }
157  else
158  {
159  Array<float> a = a_arg.float_array_value ();
160 
161  if (! error_state)
162  {
163  if (b_arg.is_scalar_type ())
164  {
165  float b = b_arg.float_value ();
166 
167  if (! error_state)
168  retval = betainc (x, a, b);
169  }
170  else
171  {
172  Array<float> b = b_arg.float_array_value ();
173 
174  if (! error_state)
175  retval = betainc (x, a, b);
176  }
177  }
178  }
179  }
180  }
181  else
182  {
183  if (x_arg.is_scalar_type ())
184  {
185  double x = x_arg.double_value ();
186 
187  if (a_arg.is_scalar_type ())
188  {
189  double a = a_arg.double_value ();
190 
191  if (! error_state)
192  {
193  if (b_arg.is_scalar_type ())
194  {
195  double b = b_arg.double_value ();
196 
197  if (! error_state)
198  retval = betainc (x, a, b);
199  }
200  else
201  {
202  Array<double> b = b_arg.array_value ();
203 
204  if (! error_state)
205  retval = betainc (x, a, b);
206  }
207  }
208  }
209  else
210  {
211  Array<double> a = a_arg.array_value ();
212 
213  if (! error_state)
214  {
215  if (b_arg.is_scalar_type ())
216  {
217  double b = b_arg.double_value ();
218 
219  if (! error_state)
220  retval = betainc (x, a, b);
221  }
222  else
223  {
224  Array<double> b = b_arg.array_value ();
225 
226  if (! error_state)
227  retval = betainc (x, a, b);
228  }
229  }
230  }
231  }
232  else
233  {
234  Array<double> x = x_arg.array_value ();
235 
236  if (a_arg.is_scalar_type ())
237  {
238  double a = a_arg.double_value ();
239 
240  if (! error_state)
241  {
242  if (b_arg.is_scalar_type ())
243  {
244  double b = b_arg.double_value ();
245 
246  if (! error_state)
247  retval = betainc (x, a, b);
248  }
249  else
250  {
251  Array<double> b = b_arg.array_value ();
252 
253  if (! error_state)
254  retval = betainc (x, a, b);
255  }
256  }
257  }
258  else
259  {
260  Array<double> a = a_arg.array_value ();
261 
262  if (! error_state)
263  {
264  if (b_arg.is_scalar_type ())
265  {
266  double b = b_arg.double_value ();
267 
268  if (! error_state)
269  retval = betainc (x, a, b);
270  }
271  else
272  {
273  Array<double> b = b_arg.array_value ();
274 
275  if (! error_state)
276  retval = betainc (x, a, b);
277  }
278  }
279  }
280  }
281  }
282  }
283  else
284  print_usage ();
285 
286  return retval;
287 }
288 
289 /*
290 ## Double precision
291 %!test
292 %! a = [1, 1.5, 2, 3];
293 %! b = [4, 3, 2, 1];
294 %! v1 = betainc (1,a,b);
295 %! v2 = [1,1,1,1];
296 %! x = [.2, .4, .6, .8];
297 %! v3 = betainc (x, a, b);
298 %! v4 = 1 - betainc (1.-x, b, a);
299 %! assert (v1, v2, sqrt (eps));
300 %! assert (v3, v4, sqrt (eps));
301 
302 ## Single precision
303 %!test
304 %! a = single ([1, 1.5, 2, 3]);
305 %! b = single ([4, 3, 2, 1]);
306 %! v1 = betainc (1,a,b);
307 %! v2 = single ([1,1,1,1]);
308 %! x = single ([.2, .4, .6, .8]);
309 %! v3 = betainc (x, a, b);
310 %! v4 = 1 - betainc (1.-x, b, a);
311 %! assert (v1, v2, sqrt (eps ("single")));
312 %! assert (v3, v4, sqrt (eps ("single")));
313 
314 ## Mixed double/single precision
315 %!test
316 %! a = single ([1, 1.5, 2, 3]);
317 %! b = [4, 3, 2, 1];
318 %! v1 = betainc (1,a,b);
319 %! v2 = single ([1,1,1,1]);
320 %! x = [.2, .4, .6, .8];
321 %! v3 = betainc (x, a, b);
322 %! v4 = 1-betainc (1.-x, b, a);
323 %! assert (v1, v2, sqrt (eps ("single")));
324 %! assert (v3, v4, sqrt (eps ("single")));
325 
326 %!error betainc ()
327 %!error betainc (1)
328 %!error betainc (1,2)
329 %!error betainc (1,2,3,4)
330 */
331 
332 DEFUN (betaincinv, args, ,
333  "-*- texinfo -*-\n\
334 @deftypefn {Mapping Function} {} betaincinv (@var{y}, @var{a}, @var{b})\n\
335 Compute the inverse of the incomplete Beta function.\n\
336 \n\
337 The inverse is the value @var{x} such that\n\
338 \n\
339 @example\n\
340 @var{y} == betainc (@var{x}, @var{a}, @var{b})\n\
341 @end example\n\
342 @seealso{betainc, beta, betaln}\n\
343 @end deftypefn")
344 {
345  octave_value retval;
346 
347  int nargin = args.length ();
348 
349  if (nargin == 3)
350  {
351  octave_value x_arg = args(0);
352  octave_value a_arg = args(1);
353  octave_value b_arg = args(2);
354 
355  if (x_arg.is_scalar_type ())
356  {
357  double x = x_arg.double_value ();
358 
359  if (a_arg.is_scalar_type ())
360  {
361  double a = a_arg.double_value ();
362 
363  if (! error_state)
364  {
365  if (b_arg.is_scalar_type ())
366  {
367  double b = b_arg.double_value ();
368 
369  if (! error_state)
370  retval = betaincinv (x, a, b);
371  }
372  else
373  {
374  Array<double> b = b_arg.array_value ();
375 
376  if (! error_state)
377  retval = betaincinv (x, a, b);
378  }
379  }
380  }
381  else
382  {
383  Array<double> a = a_arg.array_value ();
384 
385  if (! error_state)
386  {
387  if (b_arg.is_scalar_type ())
388  {
389  double b = b_arg.double_value ();
390 
391  if (! error_state)
392  retval = betaincinv (x, a, b);
393  }
394  else
395  {
396  Array<double> b = b_arg.array_value ();
397 
398  if (! error_state)
399  retval = betaincinv (x, a, b);
400  }
401  }
402  }
403  }
404  else
405  {
406  Array<double> x = x_arg.array_value ();
407 
408  if (a_arg.is_scalar_type ())
409  {
410  double a = a_arg.double_value ();
411 
412  if (! error_state)
413  {
414  if (b_arg.is_scalar_type ())
415  {
416  double b = b_arg.double_value ();
417 
418  if (! error_state)
419  retval = betaincinv (x, a, b);
420  }
421  else
422  {
423  Array<double> b = b_arg.array_value ();
424 
425  if (! error_state)
426  retval = betaincinv (x, a, b);
427  }
428  }
429  }
430  else
431  {
432  Array<double> a = a_arg.array_value ();
433 
434  if (! error_state)
435  {
436  if (b_arg.is_scalar_type ())
437  {
438  double b = b_arg.double_value ();
439 
440  if (! error_state)
441  retval = betaincinv (x, a, b);
442  }
443  else
444  {
445  Array<double> b = b_arg.array_value ();
446 
447  if (! error_state)
448  retval = betaincinv (x, a, b);
449  }
450  }
451  }
452  }
453 
454  // FIXME: It would be better to have an algorithm for betaincinv which
455  // accepted float inputs and returned float outputs. As it is, we do
456  // extra work to calculate betaincinv to double precision and then throw
457  // that precision away.
458  if (x_arg.is_single_type () || a_arg.is_single_type ()
459  || b_arg.is_single_type ())
460  {
461  retval = Array<float> (retval.array_value ());
462  }
463  }
464  else
465  print_usage ();
466 
467  return retval;
468 }
469 
470 /*
471 %!assert (betaincinv ([0.875 0.6875], [1 2], 3), [0.5 0.5], sqrt (eps))
472 %!assert (betaincinv (0.5, 3, 3), 0.5, sqrt (eps))
473 %!assert (betaincinv (0.34375, 4, 3), 0.5, sqrt (eps))
474 %!assert (betaincinv (0.2265625, 5, 3), 0.5, sqrt (eps))
475 %!assert (betaincinv (0.14453125, 6, 3), 0.5, sqrt (eps))
476 %!assert (betaincinv (0.08984375, 7, 3), 0.5, sqrt (eps))
477 %!assert (betaincinv (0.0546875, 8, 3), 0.5, sqrt (eps))
478 %!assert (betaincinv (0.03271484375, 9, 3), 0.5, sqrt (eps))
479 %!assert (betaincinv (0.019287109375, 10, 3), 0.5, sqrt (eps))
480 
481 ## Test class single as well
482 %!assert (betaincinv ([0.875 0.6875], [1 2], single (3)), [0.5 0.5], sqrt (eps ("single")))
483 %!assert (betaincinv (0.5, 3, single (3)), 0.5, sqrt (eps ("single")))
484 %!assert (betaincinv (0.34375, 4, single (3)), 0.5, sqrt (eps ("single")))
485 
486 ## Extreme values
487 %!assert (betaincinv (0, 42, 42), 0, sqrt (eps))
488 %!assert (betaincinv (1, 42, 42), 1, sqrt (eps))
489 
490 %!error betaincinv ()
491 %!error betaincinv (1, 2)
492 */
493 
octave_value
Definition: ov.h:68
print_usage
OCTINTERP_API void print_usage(void)
Definition: defun.cc:51
octave_value::is_scalar_type
bool is_scalar_type(void) const
Definition: ov.h:657
DEFUN
#define DEFUN(name, args_name, nargout_name, doc)
Definition: defun.h:44
octave_value::float_value
float float_value(bool frc_str_conv=false) const
Definition: ov.h:762
betaincinv
double betaincinv(double y, double p, double q)
Definition: lo-specfun.cc:3295
octave_value::float_array_value
FloatNDArray float_array_value(bool frc_str_conv=false) const
Definition: ov.h:782
error_state
int error_state
Definition: error.cc:101
octave_value::length
octave_idx_type length(void) const
Definition: ov.cc:1525
betainc
double betainc(double x, double a, double b)
Definition: lo-specfun.cc:2277
octave_value::array_value
NDArray array_value(bool frc_str_conv=false) const
Definition: ov.h:779
octave_value::is_single_type
bool is_single_type(void) const
Definition: ov.h:611
octave_value::double_value
double double_value(bool frc_str_conv=false) const
Definition: ov.h:759
x
F77_RET_T const double * x
Definition: __pchip_deriv__.cc:38
Generated on Thu Jun 4 2015 23:30:12 for GNU Octave by   doxygen 1.8.8