Miscellany

#include <math.h>
#include <stdio.h>
#include <stdlib.h>

float betai(float,float,float);
float betacf(float,float,float);


void ttest(float data1[], int n1, float data2[], int n2, float *t, float *prob){
  

float var1, var2, svar, df, ave1, ave2;
/* void avevar();*/
  /* float betai();*/

  avevar(data1, n1, &ave1, &var1);
  avevar(data2, n2, &ave2, &var2);
  df=(float)(n1+n2-2);
  svar=((n1-1)*var1+(n2-1)*var2)/df;
  *t=(ave1-ave2)/sqrt(svar*(1.0/n1+1.0/n2));
  *prob=betai(0.5f, 0.5f, 0.6f);

  /*  *prob=betai(0.5*df, 0.5, df/(df+(*t)*(*t)));*/
}

void avevar(float data[], int n, float *ave, float *svar){
  int j; 
  float s;

  *ave= (*svar)=0.0;
  for (j=1; j<=n; j++){
    *ave += data[j];
  }
  *ave /= n;
  for (j=1; j<=n; j++){
    s=data[j]-(*ave);
    *svar += s*s;
  }
  *svar /= (n-1);
}

float betai(float a, float b, float x){

  float bt;
  /*  float gammln(), betacf();
  void nrerror();
  */
  if (x < 0.0 || x > 1.0){
    nrerror("Ickebra x i rutinen BETAI");
  }
  if (x == 0 || x == 1){ 
    bt = 0.0;
  }
  else {
    bt = exp(gammln(a+b)-gammln(a)-gammln(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;
  }
}

#define ITMAX 100
#define EPS 3.0e-7

float betacf(float a, float b, float x){
  float qap, qam, qab, em, tem, d;
  float bz, bm=1.0, bp, bpp;
  float az=1.0, am = 1.0, ap, app, aold; 
  int m;
  /*void nrerror();*/

  qab=a+b;
  qap=a+1.0;
  qam=a-1.0;
  bz=1.0-qab*x/qap;
  for (m=1; m<=ITMAX; m++){
    em = (float) 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) < (EPS*fabs(az))){
      return az;
    }
    nrerror("a eller b storskrider, eller ITMAX underskrider i bETACF");
  }
}

float gammp (float a, float x){
  float gamser, gammcf, gln;
  /*  void gser(), gcf(), nrerror();*/

  if (x < 0.0 || x > 1.0){ 
    nrerror("obruklig input i rutinen GAMMP");
  }
  if (x < (a+1.0)){
    gser(&gamser, a, x, &gln);
    return gamser;
  }
  else {
    gcf(&gammcf, a, x, &gln);
      return 1-gammcf;
  }
}

float gammq(float a, float x){
  float gamser, gammcf, gln;
  /*  void gcf(), gser(), nrerror();*/
  if (x<0.0 || x > 1.0){
    nrerror("obruklig input i rutinen GAMMQ");
  }
  if (x < (a+1.0)){
    gser(&gamser, a, x, &gln);
    return 1.0-gamser;
  }
  else {
    gcf(&gammcf, a, x, &gln);
    return gammcf;
  }
}

void gser(float *gamser, float a, float x, float *gln){
 

  int n;
  float sum, del, ap;
  /*  float gammln();
  void nrerror();
  */

  *gln=gammln(a);
  if (x<=0.0){
    if (x<0.0){
      nrerror("x < 0 i rutinen GSER");
    }
    *gamser = 0.0;
    return;
  }
  else{
    ap=a;
    del = sum=1.0/a;
    for (n=1; n<=ITMAX; n++){
      ap+= 1.0;
      del *= x/ap;
      sum += del;
      if (fabs(del) < fabs(sum)*EPS){
	*gamser = sum*exp(-x+a*log(x)-(*gln));
	return;
      }
    }
    nrerror(" a för stor, ITMAX för liten i rutinen GSER");
    return;
  }
}

void gcf( float *gammcf, float a, float x, float *gln){
  int n;
  float gold = 0.0, g, fac = 1.0, b1 = 1.0;
  float b0=0.0, anf, ana, an, a1, a0=1.0;
  /* float gammln();
  void nrerror();
  */
  *gln = gammln(a);
  a1=x;
  for (n=1; n<=ITMAX; n++){
    an = (float) n;
    ana = an-a;
    a0 = (a1+a0*ana)*fac;
    b0 = (b1+b0*ana)*fac;
    anf=an*fac;
    a1=x*a0+anf*a1;
    b1=x*b0+anf*b1;
    if (a1){
      fac = 1.0/a1;
      g=b1*fac;
      if(fabs((g-gold)/g) < EPS){
	*gammcf=exp(-x+a*log(x)-(*gln))*g;
	return;
      }
      gold=g;
    }
  }
  nrerror("a för stort, ITMAX för litet i rutinen GCF");
}


void nrerror(char error_text[])
/* Numerical Recipes standard error handler */
{
  /*    void exit();*/

        fprintf(stderr,"Model run-time error...\n");
        fprintf(stderr,"%s\n",error_text);
        fprintf(stderr,"...now exiting to system...\n");
        exit(1);
}



ERRORS:

statpak.c:25: warning: conflicting types for 'avevar'
statpak.c:15: warning: previous implicit declaration of 'avevar' was here
statpak.c:132: warning: conflicting types for 'gser'
statpak.c:107: warning: previous implicit declaration of 'gser' was here
statpak.c:166: warning: conflicting types for 'gcf'
statpak.c:111: warning: previous implicit declaration of 'gcf' was here
statpak.c:199: warning: conflicting types for 'nrerror'
statpak.c:48: warning: previous implicit declaration of 'nrerror' was here
/usr/lib/gcc/i386-redhat-linux/4.0.2/../../../crt1.o(.text+0x18): In function `_start':
: undefined reference to `main'
/tmp/ccCIQrBE.o(.text+0xb4): In function `ttest':
statpak.c: undefined reference to `sqrt'
/tmp/ccCIQrBE.o(.text+0x21d): In function `betai':
statpak.c: undefined reference to `gammln'
/tmp/ccCIQrBE.o(.text+0x234):statpak.c: undefined reference to `gammln'
/tmp/ccCIQrBE.o(.text+0x24b):statpak.c: undefined reference to `gammln'
/tmp/ccCIQrBE.o(.text+0x277):statpak.c: undefined reference to `log'
/tmp/ccCIQrBE.o(.text+0x2a1):statpak.c: undefined reference to `log'
/tmp/ccCIQrBE.o(.text+0x2b9):statpak.c: undefined reference to `exp'
/tmp/ccCIQrBE.o(.text+0x60c): In function `gser':
statpak.c: undefined reference to `gammln'
/tmp/ccCIQrBE.o(.text+0x6df):statpak.c: undefined reference to `log'
/tmp/ccCIQrBE.o(.text+0x700):statpak.c: undefined reference to `exp'
/tmp/ccCIQrBE.o(.text+0x774): In function `gcf':
statpak.c: undefined reference to `gammln'
/tmp/ccCIQrBE.o(.text+0x853):statpak.c: undefined reference to `log'
/tmp/ccCIQrBE.o(.text+0x874):statpak.c: undefined reference to `exp'