DataMuseum.dk

top - download

    Length: 5636 (0x1604)
    Types: TextFile
    Names: »ts.c«

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

TextFile

#include "unixstat.h"
PGM(ts,Time Series Analysis and Plots,5.0,3/2/85)
/* Copyright (c) 1982 Gary Perlman (see Copyright file) */

#define	WIDTH 70           /* default plot width */

/* OPTIONS */
double	Base = 1.0;        /* base value of the plot interval */
double	Interval = 1.0;    /* plot interval size for label */
int 	Axes = 1;          /* type of axes */
int 	Echo;              /* print out (possibly transformed) ts */
int 	Width = WIDTH;     /* width of plot */
int 	Length;            /* rescaled length of ts */
int 	Plot;              /* request a plot */
int 	Plotstyle;         /* specify plot style */
int 	Stats;             /* request summary statistics */
int 	Onumber;           /* number output */
int 	Autocor;           /* maximum lag for autocorrelations */

#ifdef SMALL_MEM
#define	MAXTSIZE    250       /* max size for raw time series */
#else
#define	MAXTSIZE   1000       /* max size for raw time series */
#endif

typedef struct
	{
	unsigned	int	size;
	float	*data;
	} TS;

#define	check_ts(ts) if(!(ts || ts->size))return(NULL)

TS *
new_ts (size)
	{
	TS	*tss;
	if ((tss = (TS *) malloc (sizeof (TS))) == NULL) return (NULL);
	if ((tss->data = (float *) malloc ((unsigned)size * sizeof (*tss->data))) == NULL)
		return (NULL);
	tss->size = 0;
	return (tss);
	}

TS *
make_ts (tss, n)
TS	*tss;
	{
	TS  	*newts = new_ts (n);
	int 	i;
	int 	tsn = tss->size;
	if (newts)
		{
		newts->size = n;
		for (i = 0; i < n; i++)
			newts->data[i] = 0.0;
		for (i = 0; i < n*tsn; i++)
			newts->data[i/tsn] += tss->data[i/n];
		for (i = 0; i < n; i++)
			newts->data[i] /= tsn;
		}
	free (tss);
	return (newts);
	}

read_ts (tss)
TS	*tss;
	{
	char	s[BUFSIZ];
	check_ts (tss);
	while (getword (s, stdin))
		{
		if (!number (s))
			ERRNUM (s,input value)
		if (tss->size >= MAXTSIZE)
			ERRSPACE (time series)
		tss->data[tss->size++] = atof (s);
		}
	}

print_num (n)
	{
	if (!Onumber) return;
	if ((n >= 0) && ((n%Onumber) == 0))
		printf ("%8g	", n*Interval + Base);
	else
		printf ("%8s	", "");
	}

print_ts (tss)
TS	*tss;
	{
	int 	ii;
	check_ts (tss);
	for (ii = 0; ii < tss->size; ii++)
		{
		print_num (ii);
		printf ("%.5g\n", tss->data[ii]);
		}
	}

#define	tabprint(label,format,value) printf ("%s	= %format\n", "label", (value))

printstats (vector, n)
float	*vector;
unsigned n;
	{
	double	sum = 0.0;
	double	ss  = 0.0;
	double	minx = vector[0];
	double	maxx = vector[0];
	int 	count = n;
	while (count-- > 0)
		{
		sum += vector[0];
		ss  += vector[0] * vector[0];
		if (vector[0] > maxx)
			maxx = vector[0];
		else if (vector[0] < minx)
			minx = vector[0];
		vector++;
		}
	tabprint (n,d,n);
	tabprint (sum,g,sum);
	tabprint (ss,g,ss);
	tabprint (min,g,minx);
	tabprint (max,g,maxx);
	tabprint (range,g,maxx-minx);
	tabprint (midpt,g,(maxx+minx)/2.0);
	tabprint (mean,g,sum/n);
	tabprint (sd,g,n <= 1 ? 0.0 : sqrt ((ss-sum*sum/n)/(n-1)));
	}

cor_ts (tss)
TS	*tss;
	{
	int 	lag;
	double	cor ();
	double	r, r2;
	double	Fr;
	double	p, pof ();
	int 	df;
	printf ("%8s %6s %6s %5s %12s %5s %6s\n",
		"Lag", "r", "r^2", "n'", "F", "df", "p");
	for (lag = 0; lag <= Autocor; lag++)
	    {
	    df = tss->size - lag - 2;
		if (lag == 0)
			{
			float	*seq;
			int 	i;
			if (seq = (float *) malloc (tss->size))
				{
				for (i = 0; i < tss->size; i++)
					seq[i] = (float) i;
				r = cor (seq, tss->data, tss->size);
				free (seq);
				}
			else r = 0.0;
			}
		else
			r = cor (tss->data, tss->data+lag, tss->size-lag);
	    r2 = r*r;
	    if (r2 == 1.0)
			{
			Fr = 0.0;
			p = 0.0;
			}
	    else
			{
			Fr = (r2 * df) / (1.0 - r2);
			p = pof (Fr, 1, df);
			}
	    printf ("%8d %6.3f %6.3f %5d %12.3f %5d %6.3f\n",
		         lag,    r,   r2, tss->size-lag, Fr, df, p);
	    }
	}

initial (argc, argv) char **argv;
	{
	extern	char	*optarg;
	extern	int 	optind;
	int 	errflg = 0;
	int 	C;
#define	OPTSTRING "ei:l:b:pP:sn:ac:w:"
	while ((C = getopt (argc, argv, OPTSTRING)) != EOF)
		switch (C)
			{
			case 'e': Echo = 1; break;
			case 'i':
				if (!number (optarg)) ERRNUM (optarg,interval size)
				Interval = atof (optarg);
				if (Interval <= FZERO)
					ERRVAL (g,Interval,interval size)
				break;
			case 'b':
				if (!number (optarg)) ERRNUM (optarg,base)
				Base = atof (optarg);
				break;
			case 'l':
				if (!INTEGER (optarg)) ERRNUM (optarg,series length)
				Length = atoi (optarg);
				if (Length <= 0) ERRVAL (d,Length,series length)
				break;
			case 'P': Plotstyle = atoi (optarg);
			case 'p': Plot = 1; break;
			case 's': Stats = 1; break;
			case 'n':
				if (!INTEGER (optarg)) ERRNUM (optarg,number option)
				Onumber = atoi (optarg);
				if (Onumber <= 0) ERRVAL (d,Onumber,number option)
				break;
			case 'a': Plot = (Axes = !Axes); break;
			case 'c':
				if (!INTEGER (optarg)) ERRNUM (optarg,autocorrelation lag)
				Autocor = atoi (optarg);
				if (Autocor <= 0) ERRVAL (d,Autocor,autocorrelation lag)
				break;
			case 'w':
				if (!INTEGER (optarg)) ERRNUM (optarg,plot width)
				Width = atoi (optarg);
				if (Width <= 1) ERRVAL (d,Width,plot width)
				break;
			default: errflg++; break;
			}
	if (errflg)
		USAGE ([-aes] [-i interval] [-l length] [-n number] [-c autocor] [-w width])
	if (!Plot && !Autocor && !Echo) Stats = 1;
	ERROPT (optind)
	}

main (argc, argv) char**argv;
	{
	TS	*tss = new_ts (MAXTSIZE);
	ARGV0;
	initial (argc, argv);
	checkstdin (Argv0);
	read_ts (tss);
	if (Length)
		tss = make_ts (tss, Length);
	if (Echo)
		print_ts (tss);
	if (Stats)
		printstats (tss->data, tss->size);
	if (Autocor)
		cor_ts (tss);
	if (Plot)
		barplot (tss->data, tss->size, Plotstyle, Axes, Onumber, Width, Base, Interval);
	free (tss);
	exit (0);
	}