Mercurial > octave

comparison scripts/specfun/betai.m @ 2303:5cffc4b8de57

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 1996-06-24 09:15:24 by jwe]
author jwe
date Mon, 24 Jun 1996 09:15:24 +0000
parents 5d29638dd524
children 2b5788792cad
comparison
equal deleted inserted replaced
2302:470c856bf55a 2303:5cffc4b8de57
1 # Copyright (C) 1996 John W. Eaton 1 ### Copyright (C) 1996 John W. Eaton
2 # 2 ###
3 # This file is part of Octave. 3 ### This file is part of Octave.
4 # 4 ###
5 # Octave is free software; you can redistribute it and/or modify it 5 ### Octave is free software; you can redistribute it and/or modify it
6 # under the terms of the GNU General Public License as published by the 6 ### under the terms of the GNU General Public License as published by
7 # Free Software Foundation; either version 2, or (at your option) any 7 ### the Free Software Foundation; either version 2, or (at your option)
8 # later version. 8 ### any later version.
9 # 9 ###
10 # Octave is distributed in the hope that it will be useful, but WITHOUT 10 ### Octave is distributed in the hope that it will be useful, but
11 # ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or 11 ### WITHOUT ANY WARRANTY; without even the implied warranty of
12 # FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License 12 ### MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
13 # for more details. 13 ### General Public License for more details.
14 # 14 ###
15 # You should have received a copy of the GNU General Public License 15 ### You should have received a copy of the GNU General Public License
16 # along with Octave; see the file COPYING. If not, write to the Free 16 ### along with Octave; see the file COPYING. If not, write to the Free
17 # Software Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. 17 ### Software Foundation, 59 Temple Place - Suite 330, Boston, MA
18 ### 02111-1307, USA.
18 19
19 function y = betai (a, b, x) 20 function y = betai (a, b, x)
20 21
21 # usage: betai (a, b, x) 22 ## usage: betai (a, b, x)
22 # 23 ##
23 # Returns the incomplete beta function 24 ## Returns the incomplete beta function
24 # betai (a, b, x) = BETA(a,b)^(-1) INT_0^x t^(a-1) (1-t)^(b-1) dt. 25 ## betai (a, b, x) = BETA(a,b)^(-1) INT_0^x t^(a-1) (1-t)^(b-1) dt.
25 # If x has more than one component, both a and b must be scalars. 26 ## If x has more than one component, both a and b must be scalars.
26 # If x is a scalar, a and b must be of compatible dimensions. 27 ## If x is a scalar, a and b must be of compatible dimensions.
27 28
28 # Written by KH (Kurt.Hornik@ci.tuwien.ac.at) on Aug 2, 1994. 29 ## Written by KH (Kurt.Hornik@ci.tuwien.ac.at) on Aug 2, 1994.
29 30
30 # Computation is based on the series expansion 31 ## Computation is based on the series expansion
31 # betai(a, b, x) 32 ## betai(a, b, x)
32 # = \frac{1}{B(a, b)} x^a 33 ## = \frac{1}{B(a, b)} x^a
33 # \sum_{k=0}^\infty \frac{(1-b)\cdots(k-b)}{a+k} \frac{x^k}{k!} 34 ## \sum_{k=0}^\infty \frac{(1-b)\cdots(k-b)}{a+k} \frac{x^k}{k!}
34 # for x <= 1/2. For x > 1/2, betai(a, b, x) = 1 - betai(b, a, 1-x). 35 ## for x <= 1/2. For x > 1/2, betai(a, b, x) = 1 - betai(b, a, 1-x).
35 36
36 if (nargin <> 3) 37 if (nargin <> 3)
37 usage (" betai (a, b, x)"); 38 usage (" betai (a, b, x)");
38 endif 39 endif
39 40
59 y = zeros (1, n); 60 y = zeros (1, n);
60 c = exp (lgamma (a+b) - lgamma (a) - lgamma (b)); 61 c = exp (lgamma (a+b) - lgamma (a) - lgamma (b));
61 62
62 y (find (x == 1)) = ones (1, sum (x == 1)); 63 y (find (x == 1)) = ones (1, sum (x == 1));
63 64
64 # Now do the series computation. The error when truncating at term K 65 ## Now do the series computation. The error when truncating at term K
65 # is always less than 2^(-K), hence the following choice of K. 66 ## is always less than 2^(-K), hence the following choice of K.
66 67
67 K = ceil (-log (eps) / log (2)); 68 K = ceil (-log (eps) / log (2));
68 k = (1:K)'; 69 k = (1:K)';
69 70
70 ind = find ((x > 0) & (x <= 1/2)); 71 ind = find ((x > 0) & (x <= 1/2));