# TITLE:  Inverse t and F probability distributions          
# AUTHOR: Roy F.A. Maclean 
# 	The beta and gamma subroutines are translated from 
#	"Numerical Recipes in Pascal" by Press, Flannery, Teukolsky and Vetterling.
# EMAIL:  rfamgm at gmail                 
# DATE:   10Jul2007
# MAKE:   CASIO
# MODEL:  9850      
# NOTES:  
#
#	For inverse t: Enter the area between t and infinity and the degree of freedom V.
#			The value of t giving rise to that area will be returned.
#			On a 9850 it takes about 30 seconds.
#
#	For inverse F: Enter the area between F and infinity and the degrees of freedom V1 and V2.
#			The value of F giving rise to that area will be returned.
#			On a 9850 it takes about 40 seconds.
#
#	\exp is the [EXP] button ( [*10^] button on french models )
#	\th is theta
#	\r is the polar r
#
#	The \ symbols at the end of the line in program GAMMALN mean don't enter a newline
#	just merge these lines together.

@@ Program "INVTFPR"

100->N
3\exp(-)3->H
5->Dim List 1
4->Dim List 2
3->Dim List 3
Fix 3
"1) T"
"2) F"?->Z
"AREA"?->\th
0->I
Z=2=>Goto 3
"V"?->V:31831->J
If \th>=.005
Then V=1=>64->J
V>=2=>10->J
V>=4=>5->J
IfEnd
Z=1=>Goto 4
Lbl 3
"V1"?->U:"V2"?->V:1\exp10->J
\th>=.01=>6367->J
\th>=.05=>255->J
Lbl 4
ClrText
Locate 1,2,"ITERATING..."
Lbl 2
(I+J)/2->\r
Z=1=>Prog "TPROB"
Z=2=>Prog "FPROB"
\th>Ans=>\r->J
\th<=Ans=>\r->I
Int 1000I=Int 1000J=>Int (.5+Frac (1000Frac I))=Int (.5+Frac (1000Frac J))=>Goto 1
Goto 2
Lbl 1
ClrText
"DONE":\r


@@ Program "TPROB"

\r->T
V/(V+T^2)->X
V/2->A:1/2->B
Prog "BETAI"
Ans/2


@@ Program "FPROB"

\r->F
V/(V+UF)->X
V/2->A:U/2->B
Prog "BETAI"


@@ Program "BETAI"

If (X<0) Or (X>1)
Then "BETAI-X MUST BE":"0<=X<=1"
Stop
IfEnd
If (X=0) Or (X=1)
Then 0->W
Else A+B->P
Prog "GAMMALN"
Ans->W:A->P
Prog "GAMMALN"
W-Ans->W:B->P
Prog "GAMMALN"
W-Ans->W
W+Aln X+Bln (1-X)->W
e^W->W
IfEnd
If X<(A+1)/(A+B+2)
Then Prog "BETACF"
WAns/A
Else 1-X->X
A->L:B->A:L->B
Prog "BETACF"
A->L:B->A:L->B:1-X->X
1-WAns/B
IfEnd


@@ Program "GAMMALN"

Rad:P->G:2.50662827465->S:G-1->K:K+5.5->Y
(K+.5)ln Y-Y->Y
1+76.18009173/(K+1)-86.50532033/(K+2)		\
+24.01409822/(K+3)-1.231739516/(K+4)		\
+(.120858003\exp(-)2)/(K+5)			\
-(.536382\exp(-)5)/(K+6)			
Y+ln SAns


@@ Program "BETACF"

1->List 1[1]
1->List 2[1]
1->List 1[4]
A-1->List 3[1]
A+1->List 3[2]
A+B->List 3[3]
1-XList 3[3]/List 3[2]->List 2[4]
For 1->M To N
M->E
E+E->L
EX(B-M)/((List 3[1]+L)(A+L))->D
List 1[4]+DList 1[1]->List 1[2]
List 2[4]+DList 2[1]->List 2[2]
(-)X(A+E)(List 3[3]+E)/((A+L)(List 3[2]+L))->D
List 1[2]+DList 1[4]->List 1[3]
List 2[2]+DList 2[4]->List 2[3]
List 1[4]->List 1[5]
List 1[2]/List 2[3]->List 1[1]
List 2[2]/List 2[3]->List 2[1]
List 1[3]/List 2[3]->List 1[4]
1->List 2[4]
Abs (List 1[4]-List 1[5])<HAbs List 1[4]=>Break
Next
If M=N+1
Then "A OR B TOO BIG":"OR N TOO SMALL"
Stop
IfEnd
Lbl 9
List 1[4]