/*
 *  test of equality of two binomial proportions
 *
 */
// written by ++pac following the text of http://www.ocair.org/files/KnowledgeBase/Statistics/Proportion.htm
// the following should hold: z^2 == t^2 == chisq

#include <u.h>
#include <libc.h>
#include <bio.h>
#include <ctype.h>

double n1, n2, total1, total2, p1, q1, p2, q2, fo1, fo2, fo3, fo4, fe1, fe2, fe3, fe4, diff, val2, val3, val, low, high, sigma, pu, x1, x2, x3, x4, t, z, chisq ;
Biobuf bin;

void intro(void);
int readparams(double *n1, double *total1, double *n2, double *total2);
void confid_int(void);
void t_test(void);
void z_test(void);
void chi_square(void);
void results(void);


int
readparams(double *n1, double *total1, double *n2, double *total2)
// reads n1, N1, n2, N2 for two binomial samples
{
	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, 4);
		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
	*n1=atof(av[0]);
	*total1=atof(av[1]);
	*n2=atof(av[2]);
	*total2=atof(av[3]);
	return 0;
} /*  readparams  */

void
confid_int(void) {
  p1 = n1/total1;
  p2 = n2/total2;
  q1 = 1-p1;
  q2 = 1-p2;
  fo1=n1;
  fo2=total1-n1;
  fo3=n2;
  fo4=total2-n2;

  fe1 = total1*(n1+n2)/(total1+total2);
  fe2 = total1*(total1-n1+total2-n2)/(total1+total2);
  fe3 = total2*(n1+n2)/(total1+total2);
  fe4 = total2*(total1-n1+total2-n2)/(total1+total2);
  diff = p1-p2;
  val = sqrt(p1*q1/total1+p2*q2/total2);
  low = diff-1.96*val;
  high = diff+1.96*val;
} /*  confid_int  */


void
t_test(void) {
  t=diff/val;
} /*  t_test  */


void
z_test(void) {
  pu = (n1+n2)/(total1+total2);
  val2=pu*(1-pu);
val3=((total1+total2))/(total1*total2);
  sigma=sqrt(val2)*sqrt(val3);
  z =  diff/sigma;
} /*  z_test  */


void
chi_square(void) {
  x1 = (fo1-fe1)*(fo1-fe1)/fe1;
  x2 = (fo2-fe2)*(fo2-fe2)/fe2;
  x3 = (fo3-fe3)*(fo3-fe3)/fe3;
  x4 = (fo4-fe4)*(fo4-fe4)/fe4;
  chisq = x1+x2+x3+x4;

} /*  chi_square  */

void
intro(void) {
  print("Binomial proportion difference test\n");
  print("please, give me n1, N1, n2, and N2,\n");
  print("whitespace separated\n");
} /*  intro  */

void
results(void) {
  print("\n\nBinomial proportion difference test\n");
  print("Confidence interval for p1-p2: (%f, %f)\n", low, high);
  print("t = %f (has t-distribution)\n", t);
  print("Z = %f (has Z-distribution)\n", z);
  print("chi-square = %f (has chi-squared-distribution)\n", chisq);
  print("the following should be approx. equal: \n");
  print("%f = %f = %f\n",  t*t, z*z, chisq);

} /*  results  */



void
main(void) {
//  intro();
  readparams(&n1, &total1, &n2, &total2);
  confid_int();
  t_test();
 z_test();
chi_square();
/*   */
  results();
exits("here");
} /*  main  */