betainc.cc

Go to the documentation of this file.
00001 /*
00002 
00003 Copyright (C) 1997-2012 John W. Eaton
00004 
00005 This file is part of Octave.
00006 
00007 Octave is free software; you can redistribute it and/or modify it
00008 under the terms of the GNU General Public License as published by the
00009 Free Software Foundation; either version 3 of the License, or (at your
00010 option) any later version.
00011 
00012 Octave is distributed in the hope that it will be useful, but WITHOUT
00013 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
00014 FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
00015 for more details.
00016 
00017 You should have received a copy of the GNU General Public License
00018 along with Octave; see the file COPYING.  If not, see
00019 <http://www.gnu.org/licenses/>.
00020 
00021 */
00022 
00023 #ifdef HAVE_CONFIG_H
00024 #include <config.h>
00025 #endif
00026 
00027 #include "lo-specfun.h"
00028 
00029 #include "defun-dld.h"
00030 #include "error.h"
00031 #include "gripes.h"
00032 #include "oct-obj.h"
00033 #include "utils.h"
00034 
00035 DEFUN_DLD (betainc, args, ,
00036   "-*- texinfo -*-\n\
00037 @deftypefn {Mapping Function} {} betainc (@var{x}, @var{a}, @var{b})\n\
00038 Return the regularized incomplete Beta function,\n\
00039 @tex\n\
00040 $$\n\
00041  I (x, a, b) = {1 \\over {B (a, b)}} \\int_0^x t^{(a-z)} (1-t)^{(b-1)} dt.\n\
00042 $$\n\
00043 @end tex\n\
00044 @ifnottex\n\
00045 @c Set example in small font to prevent overfull line\n\
00046 \n\
00047 @smallexample\n\
00048 @group\n\
00049 @c spacing appears odd here, but is correct after Makeinfo\n\
00050                                      x\n\
00051                           1         /\n\
00052 betainc (x, a, b) = -----------    | t^(a-1) (1-t)^(b-1) dt.\n\
00053                      beta (a, b)    /\n\
00054                                  t=0\n\
00055 @end group\n\
00056 @end smallexample\n\
00057 \n\
00058 @end ifnottex\n\
00059 \n\
00060 If @var{x} has more than one component, both @var{a} and @var{b} must be\n\
00061 scalars.  If @var{x} is a scalar, @var{a} and @var{b} must be of\n\
00062 compatible dimensions.\n\
00063 @end deftypefn")
00064 {
00065   octave_value retval;
00066 
00067   int nargin = args.length ();
00068 
00069   if (nargin == 3)
00070     {
00071       octave_value x_arg = args(0);
00072       octave_value a_arg = args(1);
00073       octave_value b_arg = args(2);
00074 
00075       // FIXME Can we make a template version of the duplicated code below
00076       if (x_arg.is_single_type () || a_arg.is_single_type () ||
00077           b_arg.is_single_type ())
00078         {
00079           if (x_arg.is_scalar_type ())
00080             {
00081               float x = x_arg.float_value ();
00082 
00083               if (a_arg.is_scalar_type ())
00084                 {
00085                   float a = a_arg.float_value ();
00086 
00087                   if (! error_state)
00088                     {
00089                       if (b_arg.is_scalar_type ())
00090                         {
00091                           float b = b_arg.float_value ();
00092 
00093                           if (! error_state)
00094                             retval = betainc (x, a, b);
00095                         }
00096                       else
00097                         {
00098                           FloatNDArray b = b_arg.float_array_value ();
00099 
00100                           if (! error_state)
00101                             retval = betainc (x, a, b);
00102                         }
00103                     }
00104                 }
00105               else
00106                 {
00107                   FloatNDArray a = a_arg.float_array_value ();
00108 
00109                   if (! error_state)
00110                     {
00111                       if (b_arg.is_scalar_type ())
00112                         {
00113                           float b = b_arg.float_value ();
00114 
00115                           if (! error_state)
00116                             retval = betainc (x, a, b);
00117                         }
00118                       else
00119                         {
00120                           FloatNDArray b = b_arg.float_array_value ();
00121 
00122                           if (! error_state)
00123                             retval = betainc (x, a, b);
00124                         }
00125                     }
00126                 }
00127             }
00128           else
00129             {
00130               FloatNDArray x = x_arg.float_array_value ();
00131 
00132               if (a_arg.is_scalar_type ())
00133                 {
00134                   float a = a_arg.float_value ();
00135 
00136                   if (! error_state)
00137                     {
00138                       if (b_arg.is_scalar_type ())
00139                         {
00140                           float b = b_arg.float_value ();
00141 
00142                           if (! error_state)
00143                             retval = betainc (x, a, b);
00144                         }
00145                       else
00146                         {
00147                           FloatNDArray b = b_arg.float_array_value ();
00148 
00149                           if (! error_state)
00150                             retval = betainc (x, a, b);
00151                         }
00152                     }
00153                 }
00154               else
00155                 {
00156                   FloatNDArray a = a_arg.float_array_value ();
00157 
00158                   if (! error_state)
00159                     {
00160                       if (b_arg.is_scalar_type ())
00161                         {
00162                           float b = b_arg.float_value ();
00163 
00164                           if (! error_state)
00165                             retval = betainc (x, a, b);
00166                         }
00167                       else
00168                         {
00169                           FloatNDArray b = b_arg.float_array_value ();
00170 
00171                           if (! error_state)
00172                             retval = betainc (x, a, b);
00173                         }
00174                     }
00175                 }
00176             }
00177         }
00178       else
00179         {
00180           if (x_arg.is_scalar_type ())
00181             {
00182               double x = x_arg.double_value ();
00183 
00184               if (a_arg.is_scalar_type ())
00185                 {
00186                   double a = a_arg.double_value ();
00187 
00188                   if (! error_state)
00189                     {
00190                       if (b_arg.is_scalar_type ())
00191                         {
00192                           double b = b_arg.double_value ();
00193 
00194                           if (! error_state)
00195                             retval = betainc (x, a, b);
00196                         }
00197                       else
00198                         {
00199                           NDArray b = b_arg.array_value ();
00200 
00201                           if (! error_state)
00202                             retval = betainc (x, a, b);
00203                         }
00204                     }
00205                 }
00206               else
00207                 {
00208                   NDArray a = a_arg.array_value ();
00209 
00210                   if (! error_state)
00211                     {
00212                       if (b_arg.is_scalar_type ())
00213                         {
00214                           double b = b_arg.double_value ();
00215 
00216                           if (! error_state)
00217                             retval = betainc (x, a, b);
00218                         }
00219                       else
00220                         {
00221                           NDArray b = b_arg.array_value ();
00222 
00223                           if (! error_state)
00224                             retval = betainc (x, a, b);
00225                         }
00226                     }
00227                 }
00228             }
00229           else
00230             {
00231               NDArray x = x_arg.array_value ();
00232 
00233               if (a_arg.is_scalar_type ())
00234                 {
00235                   double a = a_arg.double_value ();
00236 
00237                   if (! error_state)
00238                     {
00239                       if (b_arg.is_scalar_type ())
00240                         {
00241                           double b = b_arg.double_value ();
00242 
00243                           if (! error_state)
00244                             retval = betainc (x, a, b);
00245                         }
00246                       else
00247                         {
00248                           NDArray b = b_arg.array_value ();
00249 
00250                           if (! error_state)
00251                             retval = betainc (x, a, b);
00252                         }
00253                     }
00254                 }
00255               else
00256                 {
00257                   NDArray a = a_arg.array_value ();
00258 
00259                   if (! error_state)
00260                     {
00261                       if (b_arg.is_scalar_type ())
00262                         {
00263                           double b = b_arg.double_value ();
00264 
00265                           if (! error_state)
00266                             retval = betainc (x, a, b);
00267                         }
00268                       else
00269                         {
00270                           NDArray b = b_arg.array_value ();
00271 
00272                           if (! error_state)
00273                             retval = betainc (x, a, b);
00274                         }
00275                     }
00276                 }
00277             }
00278         }
00279     }
00280   else
00281     print_usage ();
00282 
00283   return retval;
00284 }
00285 
00286 /*
00287 
00288 %% test/octave.test/arith/betainc-1.m
00289 %!test
00290 %! a=[1, 1.5, 2, 3];
00291 %! b=[4, 3, 2, 1];
00292 %! v1=betainc(1,a,b);
00293 %! v2=[1,1,1,1];
00294 %! x = [.2, .4, .6, .8];
00295 %! v3=betainc(x, a, b);
00296 %! v4 = 1-betainc(1.-x, b, a);
00297 %! assert(v1, v2, sqrt(eps));
00298 %! assert(v3, v4, sqrt(eps));
00299 
00300 %% Single precision
00301 %!test
00302 %! a=single ([1, 1.5, 2, 3]);
00303 %! b=single ([4, 3, 2, 1]);
00304 %! v1=betainc(1,a,b);
00305 %! v2=single ([1,1,1,1]);
00306 %! x = single ([.2, .4, .6, .8]);
00307 %! v3=betainc(x, a, b);
00308 %! v4 = 1-betainc(1.-x, b, a);
00309 %! assert(v1, v2, sqrt(eps ('single')));
00310 %! assert(v3, v4, sqrt(eps ('single')));
00311 
00312 %% Mixed double/single precision
00313 %!test
00314 %! a=single ([1, 1.5, 2, 3]);
00315 %! b=[4, 3, 2, 1];
00316 %! v1=betainc(1,a,b);
00317 %! v2=single ([1,1,1,1]);
00318 %! x = [.2, .4, .6, .8];
00319 %! v3=betainc(x, a, b);
00320 %! v4 = 1-betainc(1.-x, b, a);
00321 %! assert(v1, v2, sqrt(eps ('single')));
00322 %! assert(v3, v4, sqrt(eps ('single')));
00323 
00324 %% test/octave.test/arith/betainc-2.m
00325 %!error <Invalid call to betainc> betainc();
00326 
00327 %% test/octave.test/arith/betainc-3.m
00328 %!error <Invalid call to betainc> betainc(1);
00329 
00330 %% test/octave.test/arith/betainc-4.m
00331 %!error <Invalid call to betainc> betainc(1,2);
00332 
00333 */
 All Classes Files Functions Variables Typedefs Enumerations Enumerator Friends Defines