/***********************************************************************
 * Program to compute single-factor analysis of variance
 * with unequal replication.
 *
 * Author: M.V.Wickerhauser
 * Date: 22 February 2006
 *
 */
#include <assert.h>
#include <stdio.h>
#include <stdlib.h>

void 
usage(char *name) 
{
  fprintf(stderr, "Usage: %s k n1 n2 ... nk < datafile\n", name);
  fprintf(stderr, " where    k  = number of factor groups (1<k<100)\n");
  fprintf(stderr, "          n1 = number of observations in group 1\n");
  fprintf(stderr, "          n2 = number of observations in group 2\n");
  fprintf(stderr, "          nk = number of observations in group k\n");
  fprintf(stderr, "    datafile = text file of observations in the following order:\n");
  fprintf(stderr, "       x11  x12  ... x1n1 \n");
  fprintf(stderr, "       x21  x22 ... x2n2  \n");
  fprintf(stderr, "       ...  ...  ...  \n");
  fprintf(stderr, "       xk1  xk2   ...  xknk \n");
  exit(1);
}

int
main(int argc, char **argv) {
  int k, n[100], i, l, N;
  
  float **x;			/* Observations x[i][l] go here */
  float mX, *mXi;		/* Total and group means */
  float totalSS, groupSS, errorSS; /* Sums of squares */
  float totalMS, groupMS, errorMS; /* Mean squares */
  int totalDF, groupDF, errorDF; /* Degrees of freedom */

  /* CHECK COMMAND LINE ARGUMENTS: */
  if(argc<2) usage(argv[0]);
  k=atoi(argv[1]);		/* number of groups */
  if(k<2 || k>99) usage(argv[0]); /* accept only `1<k<100' groups */
  printf("Number of groups: k=%d\n", k);
  if(argc!= 1+1+k) usage(argv[0]);
  /* Make space for the data */
  x = (float **)malloc(k*sizeof(float *)); assert(x);

  /* READ THE DATA as `k' rows from `stdin' */
  for(i=0; i<k; i++) {
    n[i]=atoi(argv[i+2]);	/* `n[i]' replications in group `i' */
    x[i] = (float *)malloc(n[i]*sizeof(float)); assert(x[i]);
    for(l=0; l<n[i]; l++) scanf("%f", x[i]+l);
  }
  /* Count the total degrees of freedom */
  N=0;
  for(i=0; i<k; i++) N += n[i];
  printf("Total number of samples: N=%d\n", N);


  /* ******************** MEANS ******************** */
  /* TOTAL MEAN mX */
  mX = 0;
  for(i=0; i<k; i++) 
    for(l=0; l<n[i]; l++) mX += x[i][l];
  mX /= N;
  printf("Total mean: %10.3f\n", mX);

  /* Make space for the group means */
  mXi = (float *)malloc(k*sizeof(float)); assert(mXi);
  /* GROUP MEANS mXi */
  for(i=0; i<k; i++) {
    mXi[i] = 0;
    for(l=0; l<n[i]; l++)
      mXi[i] += x[i][l];
    mXi[i] /= n[i];
    printf("Group %d mean (%d replications): %10.3f\n", i+1, n[i], mXi[i]);
  }


  /* ******************** SUMS OF SQUARES ******************** */
  /* Total sum of squares */
  totalSS = 0;
  for(i=0; i<k; i++)
    for(l=0; l<n[i]; l++)
      totalSS += (x[i][l] - mX)*(x[i][l] - mX);

  /* Groups sum of squares */
  groupSS = 0;
  for(i=0; i<k; i++)
    groupSS += (mXi[i] - mX)*(mXi[i] - mX) * n[i];

  /* Error sum of squares */
  errorSS = 0;
  for(i=0; i<k; i++)
    for(l=0; l<n[i]; l++)
	errorSS += (x[i][l] - mXi[i])*(x[i][l] - mXi[i]);

  /* ******************** DEGREES OF FREEDOM ******************** */
  totalDF = N-1;
  groupDF = k-1;
  errorDF = N-k;

  /* ******************** MEAN SQUARE DEVIATION ******************** */

  totalMS = totalSS/totalDF;
  errorMS = errorSS/errorDF;
  groupMS = groupSS/groupDF;


  /* ******************** ANOVA OUTPUT ******************** */

  printf("----------------------------------------------\n");
  printf("|           SINGLE-FACTOR ANOVA              |\n");
  printf("----------------------------------------------\n");
  printf("Source             SS      DF               MS\n");
  printf("----------------------------------------------\n");

  printf("Total    %12.3f     %3d     %12.3f\n", totalSS, totalDF, totalMS );
  printf("Group    %12.3f     %3d     %12.3f\n", groupSS, groupDF, groupMS );
  printf("Error    %12.3f     %3d     %12.3f\n", errorSS, errorDF, errorMS );

  printf("-------------------------------------------------------------\n");
  printf("Source of F value            num DF    den DF            F   \n");
  printf("-------------------------------------------------------------\n");
  printf("(Group MS) / (Error MS)     %3d       %3d     %12.3f\n",
	 groupDF, errorDF, groupMS/errorMS );
  printf("-------------------------------------------------------------\n");

  return 0;
}