runs_test.c

Go to the documentation of this file.

/*****************************************************************************/
#include "libtpcmodel.h"
/*****************************************************************************/
int runs_test(
  double *data1, 
  double *data2, 
  int N, 
  double alpha,
  double *p
) {
  int R=1, Nneg=0, Npos=0, NN;
  double expect, var, Z, pvalue=0.0, a, b;
  if(alpha<=0.0) alpha=5.0;
 
  residuals(data1, data2, N, &R, &Nneg, &Npos);
  //printf("* R=%d  Nneg=%d  Npos=%d\n", R, Nneg, Npos); 
 
  /* Expectation value for R. See Cobelli p. 262*/
  NN=2*Nneg*Npos;
  if(NN<1) {if(p!=NULL) *p=pvalue; return -1;}
  expect=(double)NN/(double)N + 1.0;
 
  /* Variance for R. See Cobelli Eq 8.6.4 p. 262*/
  a=(double)NN*(double)(NN-N); b=(double)(N-1)*(double)N*(double)N;
  var=a/b;
 
  /* R->Z that tends to _standardized_ normal distribution */
  Z=((double)R-expect)/sqrt(var); // printf("  Z=%g\n", Z);
  /* Get propability of getting value Z from normal distribution*/
  pvalue=ndtr(Z); if(p!=NULL) *p=pvalue;
 
  if (pvalue<0.01*alpha) return -1;
  else return 0;
}
/*****************************************************************************/
 
/*****************************************************************************/
int residuals(
  double *d1, 
  double *d2, 
  int Nr, 
  int *Rnr, 
  int *Nminus, 
  int *Nplus 
) {
  if(d1==NULL || d2==NULL || Nr<1 || Rnr==NULL || Nminus==NULL || Nplus==NULL) return(1);
 
  int i, *sign;
  double *res;
 
  /* Allocate memory for residuals */
  res=(double*)malloc((Nr)*sizeof(double));
  sign=(int*)malloc((Nr)*sizeof(int));
  if(res==NULL || sign==NULL) return(2);
 
  /* get residuals into res, nr of negative residuals in Nminus */
  /* nr of positive residuals into Nplus */
  /* Sign of each residual is saved into sign */
  /* and if sign is different from the last residual, increase runs nr */
  *Nminus=0; *Nplus=0; *Rnr=1;
  for(i=0; i<Nr; i++) {
    res[i]=d1[i]-d2[i];
    if(res[i]<0.0) {
      *Nminus+=1; sign[i]=-1;
    } else {
      *Nplus+=1; sign[i]=1;
    }
    if(i>0) if(sign[i]!=sign[i-1]) *Rnr+=1;
  }
  free(res); free(sign);
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/
int mrl_between_tacs(
  /* Pointer to the first data array */
  double y1[],
  /* Pointer to the second data array */
  double y2[],
  /* Length of data arrays */
  int n
) {
  int mrl=0, rl=0;
  char last_sign=0, sign;
 
  if(y1==NULL || y2==NULL) return(0);
  for(int i=0; i<n; i++) {
    if(isnan(y1[i]) || isnan(y2[i])) continue;
    if(y1[i]>y2[i]) sign=1; else if(y1[i]<y2[i]) sign=-1; else sign=0;
    if(sign!=last_sign) {
      rl=0; last_sign=sign;
    } else {
      if(sign!=0) {rl++; if(rl>mrl) mrl=rl;}
    }
  }
  return(mrl);
}
/*****************************************************************************/
 
/*****************************************************************************/