/*
 *  Student's t-test, returns probability of x >= t
 *
 */

#include <u.h>
#include <libc.h>
#include <bio.h>
#include <stdio.h>
#include <ctype.h>
#include "config.h"
#include "gsl_types.h"
#include "gsl_rng.h"
#include "gsl_randist.h"
#include "gsl_cdf.h"

#pragma lib "libsys.a"
#pragma lib "libcheb.a"
#pragma lib "libspecfunc.a"
#pragma lib "librng.a"
#pragma lib "librandist.a"
#pragma lib "libcdf.a"




int df; 
double  t, probt, m1, m2, n1, n2, pv;
Biobuf bin;
int help=0;
int verbose=0;
int t_df=0;

void intro(void);
int readparams();
void t_test(void);
void results(void);


int
readparams(void)
// reads t, df for t-distribution
{
	char *buf, *av[256], first;
	int i, j, n;
	Binit(&bin, 0, OREAD);

	n=0;
	while(n==0){
		if((buf = Brdstr(&bin, '\n', 1)) == nil) return 0;			// nothing to read
		first=buf[0];
//				print("first==%c\n", first);
		n=tokenize(buf, av, 8);
		if(n!=0){
			// skip comment lines
			if((first=='#') || (first=='>')){
//				print("comments\n");
				n=0;
			}
		}
		else{
		// skip whitespace lines
//			print("whitespace\n");	
		}
	}	// ends up with n > 0
if (t_df) {
	t=atof(av[0]);
	df=atoi(av[1]);
}
else {
	m1=atof(av[0]);
	m2=atof(av[1]);
	n1=atoi(av[2]);
	n2=atoi(av[3]);
	pv=atof(av[4]);
}
	return 0;
} /*  readparams  */

void
calc_t(void) {
  /* calculate the t statistic */
  t = (m1 - m2) / (sqrt (pv * ((1.0 / n1) + (1.0 / n2)))); 
  df=n1+n2-2;
} /*  calc_t  */


void
t_test(void) {
  probt=gsl_ran_tdist_pdf(t, df);
} /*  t_test  */

void
intro(void) {
  print(" Student's t-test, \n");
  print("returns probability of x >= t\n");
if (t_df) 
    print("please, give me t, and df, ");
else
    print("please, give me mean1, mean2, n1, n2, and pooled variance, ");

  print("blank separated on a single line\n");
} /*  intro  */

void
results(void) {
  print("\n\nStudent's t-test\n\n");
  print("t = %f (has t-distribution with %d degrees of freedom)\n", t, df);
  print("Associated probability P(t, df) = %e \n", probt);

} /*  results  */



void
main(int argc, char *argv[])
{
	ARGBEGIN{
	case 'h':
		help = 1;
		break;
	case 'v':
		verbose = 1;
		break;
	case 't':
		t_df = 1;
		break;
	}ARGEND

	if(verbose) intro();
	readparams();
	if ( !t_df) calc_t();
	t_test();
	results();
} /*  main  */