/* stats.h -- header file including some basic statistical functions.
   Written by Allin Cottrell, October 1995. */

#include <stdlib.h>
#include <math.h>
#if !defined(PI)
#define PI 3.14159265358979  
#endif


double mean (double a[], double w[], int n)
	{ double wsum = 0, sum = 0; int i;
	  for (i = 0; i < n; i++) wsum += w[i];
	  for (i = 0; i < n; i++) sum += a[i] * w[i]/wsum;
	  return sum;
	}

double variance (double a[], double w[], int n)
        { double xbar = 0, ss = 0, wsum = 0; int i;
	  for (i = 0; i < n; i++) wsum += w[i];  
	  xbar = mean(a, w, n);
	  for (i = 0; i < n; i++) 
             ss+= (a[i] - xbar) * (a[i] - xbar) * w[i]/wsum;
	  return ss;
	}

double sd (double a[], double w[], int n)
	{ double v = 0, s = 0;
	  v = variance(a, w, n);
	  s = sqrt(v);
	  return s;
	}

double cv (double a[], double w[], int n)  
	{ double coeff = sd(a, w, n) / mean(a, w, n);
	  return coeff;
	}

double covariance (double a[], double b[], double w[], int n)  
	{ double xbar, ybar, ss = 0, wsum = 0 ; int i;
	  xbar = mean(a, w, n); ybar = mean(b, w, n);
	  for (i = 0; i < n; i ++) wsum += w[i];
	  for (i = 0; i < n; i++) 
             ss += (a[i] - xbar) * (b[i] - ybar) * w[i]/wsum;
	  return ss;
	}

double slope (double a[], double b[], double w[], int n)
        { return covariance(a, b, w, n) / variance(a, w, n);
	}

double intercept (double a[], double b[], double w[], int n)
        { return mean(b, w, n) - slope(a, b, w, n) * mean(a, w, n);
	}

double correlation (double a[], double b[], double w[], int n)  
	{ double xbar, ybar, ss = 0, wsum = 0; int i;
	  xbar = mean(a, w, n); ybar = mean(b, w, n);
	  for (i = 0; i < n; i ++) wsum += w[i];
	  for (i = 0; i < n; i++) 
             ss += (a[i] - xbar) * (b[i] - ybar) * w[i]/wsum;
	  return ss / (sd(a, w, n) * sd(b, w, n));
	}

double normpdf(double mu, double sigma, double x)
	{ double fraction, d;
	  d = sigma * sqrt(2.0 * PI);
	  fraction = -(x - mu) * (x - mu) / (2 * sigma * sigma);
	  return (1.0 / d) * exp(fraction);
	}

void conv(double *obs, double *w, int n, double l, 
          double h, double s, int steps, char *plotfile)
	{ double x, *fx, area, wsum, step = (h - l) / steps;
	  int i, j, m = 2 * steps;
	  FILE* cnvfile;
	  fx = malloc(m * sizeof(double));
	  x = l;
	  for (i = 0; i < n; i++) wsum += w[i];
	  for (i = 0; i <= steps; i++)
	    { x = l + i * step;
	      fx[i] = 0;
	      for (j = 0; j < n; j++) 
                 fx[i] += normpdf(obs[j], s, x) * w[j]/wsum;
	      area += fx[i];
	    }
	  for (i = 0; i < steps; i++) fx[i] *= 1.0 / area;
	  for (i = steps; i < m; i++) fx[i] = l + (i - steps) * step;
	  cnvfile = fopen(plotfile, "w");
          fprintf(cnvfile, "#  x\t\tpdf(x)\n");
          for (i = 0; i < steps; i++)
            { fprintf(cnvfile, "%f\t%f\n", fx[i + steps], fx[i]);
            }
          fclose(cnvfile);
	}

void basicstats(double a[], double w[], int n)
        { printf("\nWeighted mean = %f\n", mean(a, w, n));
	  printf("Standard Deviation = %f\n", sd(a, w, n));
	  printf("Coefficient of variation = %f\n\n", cv(a, w, n));
	}

void printcorr(double a[], double b[], double w[], int n)
        { printf("\nCorrelation coefficient = %f\n\n", 
             correlation(a, b, w, n));
        }