-- Incomplete beta function -- zlib License -- Copyright (c) 2016, 2017 Lewis Van Winkle -- https://codeplea.com/incomplete-beta-function-c incbeta(x,a,b) = test( x <= 0, 0, x >= 1, 1, x > (a + 1) / (a + b + 2), 1 - incbeta1(1 - x, b, a), incbeta1(x,a,b) ) incbeta1(x,a,b,c,d,f,i,m,n) = do( f = 1, c = 1, d = 0, for(i,0,200, m = floor(i / 2), test( i == 0, n = 1, mod(i,2) == 0, n = m (b - m) x / (a + 2 m - 1) / (a + 2 m), n = -(a + m) (a + b + m) x / (a + 2 m) / (a + 2 m + 1) ), d = 1 + n d, test(abs(d) < 10.0^(-30), d = 10.0^(-30)), d = 1 / d, c = 1 + n / c, test(abs(c) < 10.0^(-30), c = 10.0^(-30)), f = f c d, test(abs(1 - c d) < 10.0^(-8), break) ), m = lgamma(a) + lgamma(b) - lgamma(a + b), x = float(x), m = exp(log(x) a + log(1 - x) b - m) / a, test(abs(1 - c d) < 10.0^(-8), m (f - 1), tgamma(0)) ) fdist(t,df1,df2) = incbeta(t / (t + df2 / df1), 0.5 df1, 0.5 df2) tdist(t,df) = incbeta(0.5 (t+sqrt(t^2+df)) / sqrt(t^2+df), 0.5 df, 0.5 df) "Check fdist and tdist" epsilon = 10.0^(-6) check(incbeta(0,0.5,0.5) == 0) check(incbeta(1,0.5,0.5) == 1) check(fdist(0,1,1) == 0) check(abs(fdist(0.1,3,2) - 0.04710751) < epsilon) -- R: pf(0.1,3,2) check(abs(fdist(1,3,2) - 0.464758) < epsilon) check(abs(fdist(2,3,2) - 0.6495191) < epsilon) check(abs(fdist(3,3,2) - 0.7400733) < epsilon) check(abs(fdist(4,3,2) - 0.7935601) < epsilon) check(abs(tdist(-3,10) - 0.006671828) < epsilon) -- R: pt(-3,10) check(abs(tdist(-2,10) - 0.03669402) < epsilon) check(abs(tdist(-1,10) - 0.1704466) < epsilon) check(abs(tdist(0,10) - 0.5) < epsilon) check(abs(tdist(1,10) - 0.8295534) < epsilon) check(abs(tdist(2,10) - 0.963306) < epsilon) check(abs(tdist(3,10) - 0.9933282) < epsilon) "ok" "F distribution" xrange = (0,5) yrange = (0,1) draw(fdist(t,3,2),t) "t distribution" xrange = (-5,5) yrange = (0,1) draw(tdist(t,10),t)
Run