/*
* 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 */
|