Hi Andrew,
Andrew Douglas Pitonyak schrieb:
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"?
Thanks for looking at my problem. But no, I do not talk about the
cumulative case. I mean the density function, in terms with gamma():
gamma(a+b)/(gamma(a)*gamma(b)) * x^(a-1) * (1-a)^(b-1). But it makes no
difference, the critical term is in the cumulative case too.
This is the cdf_beta function in Maxima, which has GPL source code...
GPL would be unsuitable. But fortunately the normal algorithm, as you
have cited it, is well known mathematical public property.
Most people seem to quote Bosten and Battiste (1974) and their methods.
or use Didonato/Morris Algortihm 708 [1] as in Boost
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)
In normal cases continued fraction works fine. And I found no difference
whether you change to the complement at (a-1)/(a+b-2) or at mean
a/(a+b), what others do. But till now I found no good way to calculate
it in the near of mean. Both power series and continued fraction need a
large amount of iterations in that area. That is another problem to be
solved. My tests with a Basic macro of the solution BASYM in
DiDonato/Morris gave only 6 digit accuracy. Maybe the function ERFC,
which is needed in that solution, is not exact enough in OOo, I have not
investigate it nearer yet. If you have any idea for the case "x near
mean", I would be happy.
Now, take a look at the betai method from Numerical Recipes (you should
take a look at this book)
Yes, that book is a good start.
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.
That is suitable for normal cases, but has the same problem for x near 1
in b*ln(1.0-x).
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).
I fear, that there is no solution in principle. You lost the accuracy
while subtracting. OOo should have "warning" notes like Excel, where you
can bring such information to the user.
Kind regards
Regina
[1] A.R.DiDonato, A.H.Morris: Algorithm 708. Significant Digit
Computation if the Incomplete Beta Function Ratios
in ACM Transactions on Mathematical Software, Vol. 18, No 3, Sept. 1992,
Pages 360-373
Unfortunately no free access, you need a public library.
---------------------------------------------------------------------
To unsubscribe, e-mail: [EMAIL PROTECTED]
For additional commands, e-mail: [EMAIL PROTECTED]
