/********************************************************************/
/**                                                                **/
/**   chi2test.c                                                   **/
/**                                                                **/
/**   A universal generator for multivariate T-concave             **/
/**   distributions                                                **/
/**                                                                **/
/**   run a chi^2 test for normal distribution                     **/
/**                                                                **/
/**   domain:  IR^N                                                **/
/**   density: f(x_1,...x_N) = exp(-x_1^2-...-x_N^2)               **/
/**                                                                **/
/**                                                                **/
/**   authors: Josef Leydold and Wolfgang Hoermann                 **/
/**   last modified: Mon Mar  9 13:38:00 MET 1998                  **/
/**                                                                **/
/********************************************************************/              

/* the header file */
#include "tdrmv.h"

/*------------------------------------------------------------------*/
/* Compiling options in config.h                                    */
/*                                                                  */
/*    N          2, 3 or 4                                          */
/*    MODE       0 or 1                                             */
/*    RECTANGLE  0                                                  */
/*    GAMMA      1 or 2                                             */
/*                                                                  */
/*------------------------------------------------------------------*/

/* select distribution function for N(0,1)                          */
/*    1 ... use library error-function (not ANSI C)                 */
/*    0 ... use approximation (see below)                           */
#define USE_ERF 1

/* number of classes for each dimension for chi^2 test              */
#define N_CLASS    10
/* sample size for chi^2 test                                       */
#define N_SAMPLE   10000

/*------------------------------------------------------------------*/

/* functions                                                        */
/* distribution function of N(0,1)                                  */
static double dfnormal(double x); 

/* distribution function -                                          */
/*    CHI^2 distribution for n degrees of freedom.                  */
static double dfchia(double x, int n);

/* Chi^2 test for multivariate normal distribution                  */
static void chi2test(C_ROOT *cr);

/* Chi^2 test for uniform distribution                              */
static double chi2test_uniform(int *b, int classes, int size);

/* write error message into log file                                */
extern void warning(char *msg, char *filename, int line);
/* fatal error, exit                                                */
extern void fatal(char *msg, char *filename, int line);

/*------------------------------------------------------------------*/


/*******************************************************************

  Main

  ******************************************************************/

int main()
{
  C_ROOT *hat;                       /* pointer to hat */

  double normal();                   /* density function */
  void grad_normal();                /* gradient of density */

  /* compute hat function */
  hat = get_hat(get_constr_default(),normal,grad_normal);

  /* run chi^2 test */
  chi2test(hat);
      
#if DEBUG > 0
  db_close(hat);
#endif
  
  exit(0);
} /* end of main() */

/*-----------------------------------------------------------------*/

/*******************************************************************

  normal distribution

  density: f(x_1,...x_N) = exp(-x_1^2-...-x_N^2)              

  ******************************************************************/

/*-----------------------------------------------------------------*/

double normal(double *coord)
/*
  normal density function
*/
{
  int i;             /* aux counter */
  double f=0;       /* density function */

  /* gaussian density */
  for( i=0; i<N; i++ )
    f += coord[i] * coord[i];

  return( exp(-f) );
} /* end of normal() */

/*-----------------------------------------------------------------*/

void grad_normal(double *coord, double *grad)
/*
  gradient of normal density function
*/
{
  int i;              /* aux counter */
  double f=0;             /* density at point */

  /* d(e^{-\sum j*x_j^2})/d x_i = -2*i*x_i*e^{-\sum j*x_j^2} */

  for( i=0; i<N; i++ )
    f += coord[i] * coord[i];
  f = exp( -f );

  for( i=0; i<N; i++ )
    grad[i] = (-2) * coord[i] * f;

} /* end of grad_normal() */

/*-----------------------------------------------------------------*/

/*******************************************************************
  
  distribution functions

  ******************************************************************/

/*-----------------------------------------------------------------*/
#if USE_ERF == 1
/*-----------------------------------------------------------------*/

static double dfnormal(double x)
/*
  distribution function of N(0,1)
*/
{
   if (x>=0)
     return( (1+erf(x/1.414213562))/2 );
   else
     return(erfc(-x/1.414213562)/2);
} /* end of dfnormal() */

/*-----------------------------------------------------------------*/
#else
/*-----------------------------------------------------------------*/

static double dfnormal(double x)
/*
  distribution function of N(0,1)
  Kennedy/Gentle "Statistical Computing" p 90 - 92.

  Implemented by Wolfgang Hoermann, April 1992
*/
{
  int  help;
  double xx,x2,x3,x4,x5,x6,x7,x8,zahler,nenner;

  static double
        ep0=242.66795523053175,ep1=21.979261618294152,
        ep2=6.9963834886191355,ep3= -3.5609843701815385e-2,
        eq0=215.0588758698612,eq1=91.164905404514901,
        eq2=15.082797630407787,/*eq3=1,*/
        zp0=300.4592610201616005,zq0=300.4592609569832933,
        zp1=451.9189537118729422,zq1=790.9509253278980272,
        zp2=339.3208167343436870,zq2=931.3540948506096211,
        zp3=152.9892850469404039,zq3=638.9802644656311665,
        zp4=43.16222722205673530,zq4=277.5854447439876434,
        zp5=7.211758250883093659,zq5=77.00015293522947295,
        zp6=0.5641955174789739711,zq6=12.78272731962942351,
        zp7= -1.368648573827167067e-7,/*zq7=1,*/
        dp0= -2.99610707703542174e-3,dq0=1.06209230528467918e-2,
        dp1= -4.94730910623250734e-2,dq1=1.91308926107829841e-1,
        dp2= -2.26956593539686930e-1,dq2=1.05167510706793207,
        dp3= -2.78661308609647788e-1,dq3=1.98733201817135256,
        dp4= -2.23192459734184686e-2/*,dq4=1*/;


  if(x<0)
    help=0;
  else
    help=1;
  x=fabs(x)/sqrt(2.);
  if(x<0.5) {
    x2=x*x;
    x4=x2*x2;
    x6=x4*x2;
    zahler=ep0+ep1*x2+ep2*x4+ep3*x6;
    nenner=eq0+eq1*x2+eq2*x4+x6;
    if(help)
      return(0.5*(1+x*zahler/nenner));
    else
      return(0.5*(1-x*zahler/nenner));
  }
  else if(x<4) {
    x2=x*x;
    x3=x2*x;
    x4=x3*x;
    x5=x4*x;
    x6=x5*x;
    x7=x6*x;
    zahler=zp0+zp1*x+zp2*x2+zp3*x3+zp4*x4+zp5*x5+zp6*x6+zp7*x7;
    nenner=zq0+zq1*x+zq2*x2+zq3*x3+zq4*x4+zq5*x5+zq6*x6+x7;

    if(help)
      return(0.5*(2-exp(-x2)*zahler/nenner));
    else
      return(0.5*exp(-x2)*zahler/nenner);
  }
  else if(x<50) {
    xx=x*x;
    x2=1/xx;
    x4=x2*x2;
    x6=x4*x2;
    x8=x6*x2;
    zahler=dp0+dp1*x2+dp2*x4+dp3*x6+dp4*x8;
    nenner=dq0+dq1*x2+dq2*x4+dq3*x6+x8;
    if(help)
      return(1-0.5*(exp(-xx)/x)*
             (1/sqrt(3.141592654)+zahler/(xx*nenner)));
    else
      return(0.5*(exp(-xx)/x)*(1/sqrt(3.141592654)+zahler/(xx*nenner)));
  }
  else if(help)
    return(1.);
  else
    return(0.);
}
/*-----------------------------------------------------------------*/
#endif
/*-----------------------------------------------------------------*/

static double dfchia(double x, int n)
/*
  cumulated distribution function - CHI^2 distribution
  for n degrees of freedom.
  Approximation from Wilson & Hilferty
  Kendall Stuart I, p 371
*/
{
  double chix,y;

  y=2./(9.*n);
  chix=(exp(log(x/n)/3)-1+y)/sqrt(y);

  return(dfnormal(chix));

} /* end of dfnormal() */

/*-----------------------------------------------------------------*/

/*******************************************************************
  
  chi2-tests

  ******************************************************************/

/*-----------------------------------------------------------------*/

/* constants                                                       */
#define SQRT2    1.4142135623731

/*-----------------------------------------------------------------*/

static void chi2test(C_ROOT *cr)
/* 
   Chi2 test for uniform distribution
   returns approx. p-value
*/
{
#if N == 2
  double z[N];
  int i,j,idx[N];
  int b[N][N_CLASS];
  int bg[N_CLASS][N_CLASS];

  /* clear arrays */
  for( i=0; i<N; i++ )
    for( j=0; j<N_CLASS; j++ )
      b[i][j] = 0;
  for( idx[0]=0; idx[0]<N_CLASS; idx[0]++ )
    for( idx[1]=0; idx[1]<N_CLASS; idx[1]++ )
      bg[idx[0]][idx[1]] = 0;

  /* sample */
  for( i=0; i<N_SAMPLE; i++) {
    sample_point(cr,z);
    for( j=0; j<N; j++ ) {
      idx[j] = (int)( dfnormal(z[j] * SQRT2) * N_CLASS );
      ++b[j][idx[j]];
    }

    ++bg[idx[0]][idx[1]];
  }

  /* make test */
  for( j=0; j<N; j++ ) {
    printf("Marginal x[%d]:\n",j);
    chi2test_uniform(b[j],N_CLASS,N_SAMPLE);
  }
  printf("%d-dimension:\n",N);
  chi2test_uniform(bg[0],N_CLASS*N_CLASS,N_SAMPLE);
/*.................................................................*/
#elif N == 3
  double z[N];
  int i,j,idx[N];
  int b[N][N_CLASS];
  int bg[N_CLASS][N_CLASS][N_CLASS];

  /* clear arrays */
  for( i=0; i<N; i++ )
    for( j=0; j<N_CLASS; j++ )
      b[i][j] = 0;
  for( idx[0]=0; idx[0]<N_CLASS; idx[0]++ )
    for( idx[1]=0; idx[1]<N_CLASS; idx[1]++ )
      for( idx[2]=0; idx[2]<N_CLASS; idx[2]++ )
	bg[idx[0]][idx[1]][idx[2]] = 0;

  /* sample */
  for( i=0; i<N_SAMPLE; i++) {
    sample_point(cr,z);

    for( j=0; j<N; j++ ) {
      idx[j] = (int)( dfnormal(z[j] * SQRT2) * N_CLASS );
      ++b[j][idx[j]];
    }
    ++bg[idx[0]][idx[1]][idx[2]];
  }

  /* make test */
  for( j=0; j<N; j++ ) {
    printf("Marginal x[%d]:\n",j);
    chi2test_uniform(b[j],N_CLASS,N_SAMPLE);
  }
  printf("%d-dimension:\n",N);
  chi2test_uniform(bg[0][0],N_CLASS*N_CLASS*N_CLASS,N_SAMPLE);
/*.................................................................*/
#elif N == 4
  double z[N];
  int i,j,idx[N];
  int b[N][N_CLASS];
  int bg[N_CLASS][N_CLASS][N_CLASS][N_CLASS];

  /* clear arrays */
  for( i=0; i<N; i++ )
    for( j=0; j<N_CLASS; j++ )
      b[i][j] = 0;
  for( idx[0]=0; idx[0]<N_CLASS; idx[0]++ )
    for( idx[1]=0; idx[1]<N_CLASS; idx[1]++ )
      for( idx[2]=0; idx[2]<N_CLASS; idx[2]++ )
	for( idx[3]=0; idx[3]<N_CLASS; idx[3]++ )
	  bg[idx[0]][idx[1]][idx[2]][idx[3]] = 0;

  /* sample */
  for( i=0; i<N_SAMPLE; i++) {
    sample_point(cr,z);

    for( j=0; j<N; j++ ) {
      idx[j] = (int)( dfnormal(z[j] * SQRT2) * N_CLASS );
      ++b[j][idx[j]];
    }

    ++bg[idx[0]][idx[1]][idx[2]][idx[3]];
  }

  /* make test */
  for( j=0; j<N; j++ ) {
    printf("Marginal x[%d]:\n",j);
    chi2test_uniform(b[j],N_CLASS,N_SAMPLE);
  }
  printf("%d-dimension:\n",N);
  chi2test_uniform(bg[0][0][0],N_CLASS*N_CLASS*N_CLASS*N_CLASS,N_SAMPLE);
#else
  FATAL("N must be 2, 3 or 4");
#endif
} /* end of chi2test() */


/*-----------------------------------------------------------------*/

static double chi2test_uniform(int *b, int classes, int size)
/* 
   chi2 test for uniform distribution
   returns p-value (value of power function of test)
*/
{
  double chi2 = 0.;    /* chi2 value */
  double pval;         /* p value */
  double expected;     /* expected number of observations */
  int i;

  /* compute number of expected observations */
  expected = (double)size/classes;
  if ( expected < 5. )
    WARNING("chi2test: Expected number of observations per class < 5");

  /* compute chi2 value */
  for( i=0; i<classes; i++ )
    chi2 += (b[i]-expected)*(b[i]-expected);
  chi2 /= expected;

  /* compute p-value */
  pval=1.-dfchia(chi2,classes-1);

  /* print result */
  printf("\tsample  = %d \t\t chi2-value      = %f\n",size,chi2);
  printf("\tclasses = %d     \t\t approx. P-value = %f\n\n",classes,pval);

  return pval;
} /* end of chi2test_uniform() */

/*-----------------------------------------------------------------*/