GNU Octave  3.8.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-2013 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 // FIXME: These functions do not need to be dynamically loaded. They should
36 // be placed elsewhere in the Octave code hierarchy.
37 
38 DEFUN (betainc, args, ,
39  "-*- texinfo -*-\n\
40 @deftypefn {Mapping Function} {} betainc (@var{x}, @var{a}, @var{b})\n\
41 Return the regularized incomplete Beta function,\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, i.e., @var{x} such that\n\
336 \n\
337 @example\n\
338 @var{y} == betainc (@var{x}, @var{a}, @var{b}) \n\
339 @end example\n\
340 @seealso{betainc, beta, betaln}\n\
341 @end deftypefn")
342 {
343  octave_value retval;
344 
345  int nargin = args.length ();
346 
347  if (nargin == 3)
348  {
349  octave_value x_arg = args(0);
350  octave_value a_arg = args(1);
351  octave_value b_arg = args(2);
352 
353  if (x_arg.is_scalar_type ())
354  {
355  double x = x_arg.double_value ();
356 
357  if (a_arg.is_scalar_type ())
358  {
359  double a = a_arg.double_value ();
360 
361  if (! error_state)
362  {
363  if (b_arg.is_scalar_type ())
364  {
365  double b = b_arg.double_value ();
366 
367  if (! error_state)
368  retval = betaincinv (x, a, b);
369  }
370  else
371  {
372  Array<double> b = b_arg.array_value ();
373 
374  if (! error_state)
375  retval = betaincinv (x, a, b);
376  }
377  }
378  }
379  else
380  {
381  Array<double> a = a_arg.array_value ();
382 
383  if (! error_state)
384  {
385  if (b_arg.is_scalar_type ())
386  {
387  double b = b_arg.double_value ();
388 
389  if (! error_state)
390  retval = betaincinv (x, a, b);
391  }
392  else
393  {
394  Array<double> b = b_arg.array_value ();
395 
396  if (! error_state)
397  retval = betaincinv (x, a, b);
398  }
399  }
400  }
401  }
402  else
403  {
404  Array<double> x = x_arg.array_value ();
405 
406  if (a_arg.is_scalar_type ())
407  {
408  double a = a_arg.double_value ();
409 
410  if (! error_state)
411  {
412  if (b_arg.is_scalar_type ())
413  {
414  double b = b_arg.double_value ();
415 
416  if (! error_state)
417  retval = betaincinv (x, a, b);
418  }
419  else
420  {
421  Array<double> b = b_arg.array_value ();
422 
423  if (! error_state)
424  retval = betaincinv (x, a, b);
425  }
426  }
427  }
428  else
429  {
430  Array<double> a = a_arg.array_value ();
431 
432  if (! error_state)
433  {
434  if (b_arg.is_scalar_type ())
435  {
436  double b = b_arg.double_value ();
437 
438  if (! error_state)
439  retval = betaincinv (x, a, b);
440  }
441  else
442  {
443  Array<double> b = b_arg.array_value ();
444 
445  if (! error_state)
446  retval = betaincinv (x, a, b);
447  }
448  }
449  }
450  }
451 
452  // FIXME: It would be better to have an algorithm for betaincinv which
453  // accepted float inputs and returned float outputs. As it is, we do
454  // extra work to calculate betaincinv to double precision and then throw
455  // that precision away.
456  if (x_arg.is_single_type () || a_arg.is_single_type () ||
457  b_arg.is_single_type ())
458  {
459  retval = Array<float> (retval.array_value ());
460  }
461  }
462  else
463  print_usage ();
464 
465  return retval;
466 }
467 
468 /*
469 %!assert (betaincinv ([0.875 0.6875], [1 2], 3), [0.5 0.5], sqrt (eps))
470 %!assert (betaincinv (0.5, 3, 3), 0.5, sqrt (eps))
471 %!assert (betaincinv (0.34375, 4, 3), 0.5, sqrt (eps))
472 %!assert (betaincinv (0.2265625, 5, 3), 0.5, sqrt (eps))
473 %!assert (betaincinv (0.14453125, 6, 3), 0.5, sqrt (eps))
474 %!assert (betaincinv (0.08984375, 7, 3), 0.5, sqrt (eps))
475 %!assert (betaincinv (0.0546875, 8, 3), 0.5, sqrt (eps))
476 %!assert (betaincinv (0.03271484375, 9, 3), 0.5, sqrt (eps))
477 %!assert (betaincinv (0.019287109375, 10, 3), 0.5, sqrt (eps))
478 
479 ## Test class single as well
480 %!assert (betaincinv ([0.875 0.6875], [1 2], single (3)), [0.5 0.5], sqrt (eps ("single")))
481 %!assert (betaincinv (0.5, 3, single (3)), 0.5, sqrt (eps ("single")))
482 %!assert (betaincinv (0.34375, 4, single (3)), 0.5, sqrt (eps ("single")))
483 
484 ## Extreme values
485 %!assert (betaincinv (0, 42, 42), 0, sqrt (eps))
486 %!assert (betaincinv (1, 42, 42), 1, sqrt (eps))
487 
488 %!error betaincinv ()
489 %!error betaincinv (1, 2)
490 */
491 
Generated on Mon Dec 30 2013 03:04:22 for GNU Octave by   doxygen 1.8.1.2