DataMuseum.dk

top - download

    Length: 6373 (0x18e5)
    Types: TextFile
    Names: »pair.c«

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

TextFile

#include "unixstat.h"
PGM(pair,Paired Data Analysis and Plots,5.0,3/3/85)
/* Copyright (c) 1982 Gary Perlman (see Copyright file) */

#define	prettyprint(title, x, y, d) \
	printf ("%-16s %16.4f %16.4f %16.4f\n", title, x, y, d)

#ifdef SMALL_MEM
#define	MAXPAIRS         250     /* at most this many pairs */
#else
#define	MAXPAIRS        1000     /* at most this many pairs */
#endif

float	xy_data[2][MAXPAIRS];    /* data stored for plot in here */
char	Plotchar;                /* plotting character */
int 	Height = 20;             /* default height of plot */
int 	Width = 50;              /* default width of plot */ 
int 	Perline = 2;             /* number of points per line */
int 	Count = 0;               /* number of points counted in */
double	X, Y, D;                 /* X, Y, and difference */
double	sum_x, sum_y, sum_d;     /* sums of scores */
double	ss_x, ss_y, ss_d, sum_xy;/* sums of squares and cross product */
double	min_x, max_x, min_y;     /* minimums */
double	max_y, min_d, max_d;     /* maximums */
char	*name_x = "Column 1", *name_y = "Column 2";
int 	Storedata = 1;           /* assume data will be stored */
int 	Wantstats = 0;           /* are statistics wanted? */
int 	Wantplot = 0;            /* is a plot wanted? */

main (argc, argv) char **argv;
	{
	ARGV0;
	initial (argc, argv);
	readdata ();
	if (Wantstats) printstats ();
	if (Wantplot)
		scatterplot (xy_data[0], xy_data[1], Count, Plotchar, Height, Width, 2);
	exit (0);
	}

initial (argc, argv) char **argv;
	{
	extern	int 	optind;
	extern	char 	*optarg;
	int 	C;
	int 	opterr = 0;
	if (checkargv ("biplot", Argv0))
		{
		Wantstats = 0;
		Wantplot = 1;
		}
	else if (argc == 1)
		{
		Wantstats = 1;
		Wantplot = 0;
		}
#define	OPTSTRING "c:sbpw:h:x:y:"
	else while ((C = getopt (argc, argv, OPTSTRING)) != EOF)
		switch (C)
			{
			case 'c': Wantplot = 1; Plotchar = *optarg; break;
			case 's': Wantstats = 1; break;
			case 'b':
			case 'p': Wantplot = 1; break;
			case 'w':
				if (!INTEGER (optarg)) ERRNUM (optarg,plot width)
				Wantplot = 1;
				Width = atoi (optarg);
				if (Width < 1) ERRVAL (d,Width,plot width)
				break;
			case 'h':
			if (!INTEGER (optarg)) ERRNUM (optarg,plot height)
				Wantplot = 1;
				Height = atoi (optarg);
				if (Height < 1) ERRVAL (d,Height,plot height)
				break;
			case 'x': Wantstats = 1; name_x = optarg; break;
			case 'y': Wantstats = 1; name_y = optarg; break;
			default: opterr++;
			}
	if (opterr)
		USAGE ([-sbp] [-w width] [-h height] [-c char] [-x xname] [-y yname])
	ERROPT (optind)
	checkstdin (Argv0);
	if (!Wantplot)
		{
		Storedata = 0;
		Wantstats = 1;
		}
	}

readdata ()
	{
	double	atof ();
	char	line[BUFSIZ];
	char	*array[2];
	int 	fieldcount;
	while (fgets (line, BUFSIZ, stdin))
		{
		fieldcount = parseline (line, array, 3);
		if (fieldcount == 0) continue;
		if (Count == 0)
			{
			switch (Perline = fieldcount)
				{
				case 1: Perline = 1;
					min_x = 1.0;
					min_y = atof (array[0]);
					break;
				case 2: Perline = 2;
					min_x = atof (array[0]);
					min_y = atof (array[1]);
					break;
				default:
					ERRMANY (columns, 2)
				}
			max_x = min_x;
			max_y = min_y;
			min_d = max_d = min_x - min_y;
			}
		else if (fieldcount != Perline)
			ERRMSG1 (Must have 1 or 2 numbers per line (see line %d), Count+1)
		if (!number (array[0]) || ((Perline == 2) && !number (array[1])))
			ERRMSG1 (Non-numerical input at line %d, Count+1)
		if (Perline == 2)
			{
			X	= atof (array[0]);
			Y	= atof (array[1]);
			}
		else
			{
			X	= (double) (Count+1);
			Y	= atof (array[0]);
			}
		D	= X - Y;
		sum_d	+= D;
		ss_d	+= D*D;
		sum_x	+= X;
		ss_x	+= X*X;
		sum_xy	+= X*Y;
		sum_y	+= Y;
		ss_y	+= Y*Y;
		if (X > max_x) max_x = X;
		else if (X < min_x) min_x = X;
		if (Y > max_y) max_y = Y;
		else if (Y < min_y) min_y = Y;
		if (D > max_d) max_d = D;
		else if (D < min_d) min_d = D;
		if (Storedata)
		if (Count < MAXPAIRS)
			{
			xy_data[0][Count] = X;
			xy_data[1][Count] = Y;
			}
		else
			{
			WARNING (too much data for a plot)
			Storedata = 0;
			Wantstats = 1;
			}
		Count++;
		}
	}

printstats ()
	{
	double	mean_x, mean_y, mean_d;         /* means stored here */
	double	sd_x, sd_y, sd_d, standev ();   /* standard deviations */
	double	t_x, t_y, t_d;                  /* t test values */
	double	p_x, p_y, p_d, pof ();          /* probability levels */
	char	tmpstr[20];                     /* for internal formatting */
	double	a, b;                           /* intercept and slope */
	double	r, t_r, p_r;                    /* correlation, t val & prob */
	double	fcount = (double) Count;        /* double Count */
	if (Count <= 1)
		ERRDATA
	printf ("%16s %16.16s %16.16s %16s\n", "", name_x, name_y, "Difference");
	prettyprint ("Minimums", min_x, min_y, min_d);
	prettyprint ("Maximums", max_x, max_y, max_d);
	prettyprint ("Sums", sum_x, sum_y, sum_d);
	prettyprint ("SumSquares", ss_x, ss_y, ss_d);
	mean_x = sum_x / fcount;
	mean_y = sum_y / fcount;
	mean_d = sum_d / fcount;
	prettyprint ("Means", mean_x, mean_y, mean_d);
	sd_x = standev (sum_x, ss_x, Count);
	sd_y = standev (sum_y, ss_y, Count);
	sd_d = standev (sum_d, ss_d, Count);
	prettyprint ("SDs", sd_x, sd_y, sd_d);
	if (sd_x)
		{
		t_x = mean_x*sqrt(fcount)/sd_x;
		p_x = pof (t_x*t_x, 1, Count-1);
		}
	else	p_x = t_x = 0.0;
	if (sd_y)
		{
		t_y = mean_y*sqrt(fcount)/sd_y;
		p_y = pof (t_y*t_y, 1, Count-1);
		}
	else	p_y = t_y = 0.0;
	if (sd_d)
		{
		t_d = mean_d*sqrt(fcount)/sd_d;
		p_d = pof (t_d*t_d, 1, Count-1);
		}
	else	p_d = t_d = 0.0;
	sprintf (tmpstr, "t(%d)", Count-1);
	prettyprint (tmpstr, t_x, t_y, t_d);
	prettyprint ("p", p_x, p_y, p_d);
	if (sd_x > 0.0)
		{
		b = (sum_xy - sum_x*sum_y/fcount) / (ss_x - sum_x*mean_x);
		a = mean_y - b * mean_x;
		if (sd_y) r = b * sd_x / sd_y;
		else r = 0.0;
		}
	else a = b = r = 0.0;
	if (Count > 2 && r != 1.0 && r != -1.0)
		{
		t_r = r / sqrt ((1.0 - r*r) / (fcount - 2.0));
		p_r = pof (t_r*t_r, 1, Count-2);
		}
	else	t_r = p_r = 0.0;
	putchar ('\n');
	sprintf (tmpstr, "t(%d)", Count-2);
	printf ("%16s %16s %16s %16s\n", "Correlation", "r-squared", tmpstr, "p");
	printf ("%16.4f %16.4f %16.4f %16.4f\n", r, r*r, t_r, p_r);
	printf ("%16s %16s\n", "Intercept", "Slope");
	printf ("%16.4f %16.4f\n", a, b);
	}

double
standev (sum, ss, count) double sum, ss;
	{
	if (count <= 1) return (0.0);
	return (sqrt ((ss-sum*sum/count)/(count-1)));
	}