DataMuseum.dk

top - download

    Length: 5178 (0x143a)
    Types: TextFile
    Names: »oneway.c«

└─⟦87ddcff64⟧ Bits:30001253 CPHDIST85 Tape, 1985 Autumn Conference Copenhagen
    └─ ⟦this⟧ »cph85dist/stat/src/oneway.c« 

TextFile

#include "unixstat.h"
PGM(oneway,one-way between groups t-test/anova,1.0,3/9/85)

main (argc, argv) char **argv;
	{
	ARGV0;
	initial (argc, argv);
	checkstdin (Argv0);
	readdata ();
	printstats ();
	oneway ();
	exit (0);
	}

#define	MAXGROUP    20
char	*Name[MAXGROUP];     /* names of the different groups */
int 	Ngroups = 0;         /* number of groups of data */
int 	Count[MAXGROUP];     /* number of cells in each group */
double	Sum[MAXGROUP];       /* sum of all values in the group */
double	Sumsq[MAXGROUP];     /* sums of squares */
double	Min[MAXGROUP];       /* min value in the group */
double	Max[MAXGROUP];       /* max value in the group */

/* OPTIONS */
double	Splitter = -1.0;    /* used to indicate group differences */
int 	Ttest = 0;          /* Ttest format for two group case */

#define	OPTSTRING "ts:"
initial (argc, argv) char **argv;
	{
	extern	char	*optarg;
	extern	int 	optind;
	int 	errflg = 0;
	int 	C;
	int 	group;
	while ((C = getopt (argc, argv, OPTSTRING)) != EOF)
		switch (C)
			{
			case 's':
				if (!number (optarg)) ERRNUM (optarg,group splitter)
				Splitter = atof (optarg);
				break;
			case 't': Ttest = 1; break;
			default: errflg++; break;
			}
	if (errflg)
		USAGE ([-t] [-s splitter] [names])
	for (group = 0; optind + group < argc; group++)
		Name[group] = argv[optind+group];
	}

readdata ()
	{
	char	word[BUFSIZ];       /* string input data read in here */
	double	datum;              /* each data value */
	int 	grindex = 0;        /* group index == Ngroups - 1 */
	while (getword (word, stdin))
		{
		if (!number (word)) ERRNUM (word,input data)
		datum = atof (word);
		if (datum == Splitter) /* new group */
			{
			if (Count[grindex] != 0) /* until some data read for this group */
				grindex = Ngroups++;
			}
		else /* real datum */
			{
			if (grindex >= MAXGROUP)
				ERRMANY (groups,MAXGROUP)
			if (Count[grindex] == 0)
				{
				if (Ngroups == 0) Ngroups = 1;
				Min[grindex] = Max[grindex] = datum;
				}
			else
				{
				if (datum < Min[grindex]) Min[grindex] = datum;
				else if (datum > Max[grindex]) Max[grindex] = datum;
				}
			Sum[grindex] += datum;
			Sumsq[grindex] += datum * datum;
			Count[grindex]++;
			}
		}
	if (Count[grindex] == 0) Ngroups--; /* ignore empty last group */
	if (Ngroups > MAXGROUP) ERRMANY (groups,MAXGROUP)
	if (Ngroups == 0) ERRDATA
	}

printstats ()
	{
	int 	group;
	if (Ngroups)
		printf ("%-9.9s %5s %8s %8s %8s %8s\n",
			"Name", "N", "Mean", "SD", "Min", "Max");
	for (group = 0; group < Ngroups; group++)
		{
		if (Name[group] && *Name[group])
			printf ("%-9.9s ", Name[group]);
		else printf ("Group-%-3d ", group+1);
		printf ("%5d ", Count[group]);
		printf ("%8.3f ", Sum[group]/Count[group]);
#define	var(count,sum,ss) (count<=1) ? 0.0 : ((ss-sum*sum/count)/(count-1))
#define	sd(count,sum,ss) sqrt (var (count,sum,ss))
		printf ("%8.3f ", sd (Count[group], Sum[group], Sumsq[group]));
		printf ("%8.3f ", Min[group]);
		printf ("%8.3f ", Max[group]);
		putchar ('\n');
		}
	}

/* compute one-way anova based on Guilford & Fruchter (5th Ed) p 239 */
oneway ()
	{
	double	Gsum;          /* grand sum */
	double	Gsumsq;        /* grand sum of squares (SStotal) */
	int 	Gcount;        /* N total */
	double	SSn;           /* will be sigma(i) Sum[i]*Sum[i]/Count[i] */
	int 	dfbetween;     /* degrees of freedom numerator */
	int 	dfwithin;      /* degrees of freedom denominator */
	double	SSbetween;     /* sum of squares between groups */
	double	SSwithin;      /* sum of squares within groups */
	double	MSbetween;     /* mean square between groups */
	double	MSwithin;      /* mean square within groups */
	double	F;             /* F ratio */
	double	p, pof ();     /* probability of F */
	int 	group;         /* loop index */
	Gsum = Gsumsq = 0.0;
	Gcount = 0;
	for (group = 0; group < Ngroups; group++)
		{
		Gsum += Sum[group];
		Gsumsq += Sumsq[group];
		Gcount += Count[group];
		}
	if (Gcount <= Ngroups) 
		ERRMSG0 (You need more than one datum per group for comparison)
	if (Ngroups < 2)
		ERRMSG0 (You need at least two groups of data for a comparison)
	SSn = 0.0;
	for (group = 0; group < Ngroups; group++)
		SSn += (Sum[group] * Sum[group]) / Count[group];
	SSbetween = SSn - Gsum*Gsum/Gcount;
	SSwithin = Gsumsq - SSn;
	dfbetween = Ngroups - 1;
	dfwithin = Gcount - Ngroups;
	MSbetween = SSbetween / dfbetween;
	MSwithin = SSwithin / dfwithin;
	if (fzero (MSwithin))
		{
		F = 999.999;
		p = 0.0;
		}
	else
		{
		F = MSbetween/MSwithin;
		p = pof (F, dfbetween, dfwithin);
		}
	if (Ngroups == 2 && (Ttest || checkargv ("ttest", Argv0)))
		printf ("t(%d) = %.3f   p = %.3f\n", dfwithin, sqrt (F), p);
	else
		{
	putchar ('\n');
		printf ("%-8s %10s %5s %10s %8s %5s\n",
			"Source", "SS", "df", "MS", "F", "p");
		printf ("%-8s ", "Between");
		printf ("%10.3f ", SSbetween);
		printf ("%5d ", dfbetween);
		printf ("%10.3f ", MSbetween);
		printf ("%8.3f ", F);
		printf ("%5.3f ", p);
		if (p < .001) putchar ('*');
		if (p < .01) putchar ('*');
		if (p < .05) putchar ('*');
		putchar ('\n');
		printf ("%-8s ", "Within");
		printf ("%10.3f ", SSwithin);
		printf ("%5d ", dfwithin);
		printf ("%10.3f ", MSwithin);
		putchar ('\n');
		}
	}