Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- BEGIN {
- ITMAX = 100
- EPS7 = 3.0e-7
- }
- function fatalError(msg, errcode) {
- print msg > "/dev/null"
- exit(errcode+0)
- }
- function gammaln(xx,
- x,y,tmp,ser,cof,j) {
- cof[0] = 76.18009172947146
- cof[1] = -86.50532032941677
- cof[2] = 24.01409824083091
- cof[3] = -1.231739572450155
- cof[4] = 0.1208650973866179e-2
- cof[5] = -0.5395239384953e-5
- if(xx == 0.0)
- fatalError("gammaln: ERR_INVALID_PARAMETER_0, parameter was: " xx, 2)
- y = x = xx
- tmp = x + 5.5
- tmp -= (x+0.5) * log(tmp)
- ser = 1.000000000190015
- for (j=0;j<=5;j++) ser += cof[j]/++y
- return -tmp + log(2.5066282746310005*ser/x);
- }
- function fabs(x) {
- return x >= 0.0 ? x : -x
- }
- function betacf(a, b, x,
- qap,qam,qab,em,tem,d,
- bz,bm,bp,bpp,
- az,am,ap,app,aold,
- m)
- {
- bm = az = am = 1.0
- qab=a+b
- qap=a+1.0
- qam=a-1.0
- bz=1.0-qab*x/qap
- for (m=1;m<=ITMAX;m++)
- {
- em=m
- tem=em+em
- d=em*(b-em)*x/((qam+tem)*(a+tem))
- ap=az+d*am
- bp=bz+d*bm
- d = -(a+em)*(qab+em)*x/((qap+tem)*(a+tem))
- app=ap+d*az
- bpp=bp+d*bz
- aold=az
- am=ap/bpp
- bm=bp/bpp
- az=app/bpp
- bz=1.0
- if (fabs(az-aold) < (EPS7 * fabs(az))) return az
- }
- return(paramError = 1)
- }
- function betai(a, b, x,
- bt) {
- if (x < 0.0 || x > 1.0)
- {
- paramError = 1 # TODO: find out what this is for.
- return(0.0)
- }
- if (x == 0.0 || x == 1.0)
- bt = 0.0
- else
- bt = exp(gammaln(a+b)-gammaln(a)-gammaln(b)+a*log(x)+b*log(1.0-x))
- if (x < (a+1.0)/(a+b+2.0))
- return (bt * betacf(a,b,x) / a)
- else
- return (1.0 - bt * betacf(b,a,1.0-x) / b)
- }
- function probT(t, df,
- bta)
- {
- bta = betai(df/2.0, 0.5, 1.0/(1.0 + t*t/df));
- if(t > 0.0) return(1.0 - 0.5 * bta);
- else if(t < 0.0) return(0.5 * bta);
- return(0.5);
- }
- function critical_value(p, df1, df2, max_val, type,
- NEWTON_EPSILON, minval,
- maxval, critval, prob)
- {
- NEWTON_EPSILON = 0.000001 # Accuracy of Newton approximation
- minval = 0.0;
- maxval = max_val;
- critval;
- prob = 0.0;
- if (p <= 0.0) return 0.0;
- else if (p >= 1.0) return maxval;
- critval = (df1 + df2) / sqrt(p); # fair first value
- while ((maxval - minval) > NEWTON_EPSILON)
- {
- prob = probT(critval, df1);
- if (prob < p) minval = critval;
- else maxval = critval;
- critval = (maxval + minval) * 0.5;
- }
- return critval;
- }
- function criticalX(p, df1,
- df2, x)
- {
- df2 = 0;
- x = critical_value((1.0 - p), df1, df2, 99999.0, type);
- if(x < NEWTON_EPSILON) x = 0.0;
- return(x);
- }
- function criticalT(p, df) {
- return criticalX(p, df)
- }
- function confidenceT(alpha, stdev, samplesize) {
- return sqrt(stdev/samplesize) * criticalT(alpha/2, samplesize-1)
- }
Add Comment
Please, Sign In to add comment