Regina Henschel wrote:
Hello,
I'm working on the beta distribution for Calc. In the case
BETADIST(x,1,b,0,1,FALSE()), which is the beta distribution density
function, the term b*(1-x)^(b-1) has to be calculated. The domain for
x is 0<=x<=1. There is no problem for x<0.99. But you get large
relative errors for x near 1. Is there any trick for calculating such
a term in this case?
Using the power series -(x+1/2x^2+1/3x^3+...) for log(1-x) is no
solution because of its very very bad convergence for x near 1.
kind regard
Regina
Do you mean the "cumulative beta probability density function"?
This is the cdf_beta function in Maxima, which has GPL source code...
Most people seem to quote Bosten and Battiste (1974) and their methods.
For example, see here:
http://books.google.com/books?id=-i3bCxwg7kUC&pg=PA105&lpg=PA105&dq=Bosten+and+Battiste+1974&source=web&ots=3M8J2xGMCh&sig=Iwt_mlsrJ-LbqnJ3Vt-3vYwlme4&hl=en&sa=X&oi=book_result&resnum=2&ct=result#PPA105,M1
That would be Statistical Computing around page 105. Or, further down,
near page 108, the note that continued fractions yields good results as
long as
x < (a-1)/(a+b-2)
In this case, they note that I_x(a,b) = 1 - I_{1-x}(b, a)
Now, take a look at the betai method from Numerical Recipes (you should
take a look at this book)
FUNCTION betai(a,b,x: real): real;
VAR
bt: real;
BEGIN
IF ((x < 0.0) OR (x > 1.0)) THEN BEGIN
writeln('pause in routine BETAI'); readln
END;
IF ((x = 0.0) OR (x = 1.0)) THEN bt := 0.0
ELSE bt := exp(gammln(a+b)-gammln(a)-gammln(b)
+a*ln(x)+b*ln(1.0-x));
IF (x < ((a+1.0)/(a+b+2.0))) THEN
betai := bt*betacf(a,b,x)/a
ELSE betai := 1.0-bt*betacf(b,a,1.0-x)/b
END;
Do you recognize that little trick in there? They then call betacf,
which performs the continued fractions.
Something similar from the Maxima source code:
;; Incomplete beta.
;; Reference:
;; Numerical Recipes in C: The Art of Scientific Computing
;; Comments: Translated from C. In R, pbeta(x,a,b)
;; Conditions: 0<=x<=1; a, b>0
(defun ibeta (x a b)
(declare (type flonum x a b))
(cond ((or (= x 0) (= x 1)) 0.0)
(t (let ((bt (exp (+ (- (lnbeta a b))
(* a (log x))
(* b (log (- 1 x)))))))
(cond ((< x (/ (+ 1.0 a) (+ a b 2.0)))
(/ (* bt (betacf x a b)) a))
(t (- 1.0 (/ (* bt (betacf (- 1.0 x) b a)) b))))))))
Then again, perhaps none of this will help you, because I see that they
are taking the log(1-x).
--
Andrew Pitonyak
My Macro Document: http://www.pitonyak.org/AndrewMacro.odt
My Book: http://www.hentzenwerke.com/catalog/oome.htm
Info: http://www.pitonyak.org/oo.php
See Also: http://documentation.openoffice.org/HOW_TO/index.html
---------------------------------------------------------------------
To unsubscribe, e-mail: [EMAIL PROTECTED]
For additional commands, e-mail: [EMAIL PROTECTED]
