#include <fstream.h>
#include <assert.h>
#include <string.h>
#include <math.h>
/*****************************************************************************
 Given an array of data in a file, this routine returns its mean, average 
 deviation adev, standard deviation sdev, skewness skew, and kurtosis, kurt.
 C.A. Bertulani  06/15/2000
******************************************************************************/

typedef double Number;

class Statistician
{
        private:
        int num;      // Number of data
        Number total;  // Sum of data
        Number high;	// Highest data value
        Number low;		// Lowest data value

        public:
        Statistician();
        void next_num(Number);   // Set number of data, total sum, lowest and largest
        int length();		// Get number of data
        Number min();       // Get minimun
        Number max();		// Get maximum
        Number mean();		// Get mean value
        Number adev(Number, Number []);  // Absolute deviation 
        Number sdev(Number, Number []);	// Standard deviation
        Number skew(Number, Number, Number []);	// Skewness
        Number kurt(Number, Number, Number []);	// Kurtosis
		Number corr(Number, Number, Number [], Number [], Number, Number); // Correlation
        void reset();

};

Statistician::Statistician() { num=0; total=0.;}

void Statistician::next_num(Number value)
{
        num++;
        total+=value;
        if(num==1) high=low=value;
        else if(value < low) low=value;
        else if(value > high) high=value;
}

int Statistician::length() { return num; }  

Number Statistician::min() { assert(num); return low; }

Number Statistician::max() { assert(num); return high; }

Number Statistician::mean() { assert(num); return total/num; }

Number Statistician::adev(Number mean, Number arr []) { 
	Number dev = 0.;
	for(int i=1;i<=num;i++)dev +=fabs(arr[i]-mean);
	return dev/num;
}

Number Statistician::sdev(Number mean, Number arr []) { 
	Number dev1 = 0.;
	Number dev2 = 0.;
	int i;
	for(i=1;i<=num;i++)dev1 +=(arr[i]-mean)*(arr[i]-mean);
	for(i=1;i<=num;i++)dev2 +=(arr[i]-mean);
	return sqrt(dev1/(num-1)-dev2*dev2/num/(num-1));
}

Number Statistician::skew(Number mean, Number s, Number arr []) { 
	Number dev = 0.;
	for(int i=1;i<=num;i++){
		dev +=(arr[i]-mean)*(arr[i]-mean)*(arr[i]-mean)/s/s/s;
	}
	return dev/num;
}

Number Statistician::kurt(Number mean, Number s, Number arr []) { 
	Number dev = 0.;
	for(int i=1;i<=num;i++){
		dev +=(arr[i]-mean)*(arr[i]-mean)*(arr[i]-mean)*(arr[i]-mean)/s/s/s/s;
	}
	return dev/num-3.;
}

Number Statistician::corr(Number mean1, Number mean2, 
						  Number arr1[], Number arr2[], Number s1, Number s2) {
	Number dev =0.;
	for(int i=1;i<=num;i++){
		dev +=(arr1[i]-mean1)*(arr2[i]-mean2)/s1/s2/(num-1);
	}
	return dev;
}

void Statistician::reset() { num = 0; total = 0; }

/*        Main program                        */
void main()
{
		void moment(Number [], int, Number *, Number *, Number *,
	                           Number *, Number *, Number *);

        Number data, *arr, *arr2, mean, sdev;

		// Prompt input
        char a[21], Q[]="Q";
		ifstream infile;

		// Generate our own set of data 
		ofstream outfile1;
		char  InputFileName1[] = "in1.dat";
		ofstream outfile2;
		char  InputFileName2[] = "in2.dat";

		outfile1.open (InputFileName1, ios::out);
        for(Number x=0.;x<3.14159;x+=3.14159/100)
			outfile1 << sin(x) << endl;  // Generate a Sine data
        outfile1.close();

		outfile2.open (InputFileName2, ios::out);
        for(x=0.;x<3.14159;x+=3.14159/100)
			outfile2 << cos(x) << endl;  // Generate a Cosine data
        outfile2.close();

        Statistician S, S2;

        cout    << "\n\t\t*** STATISTICAL PROGRAM ***"
                << "\nEnter the name of the data file or Q to quit: ";


        cin.getline(a,20);
		while(strcmp(a, Q))
        {
		   // First set of data 
		   infile.open(a);    
           if(!infile)
                cout << "\nFile doesn't exist";
           else {
			   // ind out data size, min, max, and mean
                while(infile >> data)
                {
                        S.next_num(data);
                }
                if(!S.length())
                        cout << "\nThere was no data in that file.";    
                else{
                        cout    << "\n\t\t*** Summary of Data1 ***";
                        cout    << "\n\tNumber of values: " << S.length()
                                << "\n\tSmallest value: " << S.min()
                                << "\tLargest value: " << S.max()
                                << "\n\tMean value: " << S.mean();
                        infile.close();
                }
				// 
				mean = S.mean();
				Number min = S.min();
				int n = S.length();
				Number max = S.mean();
				arr = new Number[n+1];
				infile.open(a);
				n = 1;
				while(infile >> data)
				{
					arr[n] = data;
					n +=1;
				}
				n = n-1;
				Number adev = S.adev(mean, arr);
				sdev = S.sdev(mean, arr);
				Number skew = S.skew(mean,sdev,arr);
				Number kurt = S.kurt(mean,sdev,arr);
				cout    << "\n\tAbsolute deviation: " << adev
                     << "\tStandard deviation: " << sdev
                     << "\n\tSkewness: " << skew
                         << "\tKurtosis: " << kurt;
				cout << endl;
				infile.close();
				// Compare with Numerical Recipes routine
				Number var;
				moment(arr, n, &mean, &adev, &sdev, &var, &skew, &kurt);
				cout    << "\n\t\t*** Compare with Numerical Recipes ***";
				cout    << "\n\tMean value: " << mean
					<< "\n\tAbsolute deviation: " << adev           
                     << "\tStandard deviation: " << sdev
                     << "\n\tSkewness: " << skew
                         << "\tKurtosis: " << kurt;
		   }

		   // Second set of data 
           cout << "\n\nEnter the name of the next data file or Q to quit: ";
           cin.getline(a,20);   		
		   infile.open(a);
           if(!infile)
                cout << "\nFile doesn't exist";
           else {
			   // ind out data size, min, max, and mean
                while(infile >> data)
                {
                        S2.next_num(data);
                }
                if(!S2.length())
                        cout << "\nThere was no data in that file.";    
                else{
                        cout    << "\n\t\t*** Summary of Data2 ***";
                        cout    << "\n\tNumber of values: " << S2.length()
                                << "\n\tSmallest value: " << S2.min()
                                << "\tLargest value: " << S2.max()
                                << "\n\tMean value: " << S2.mean();
                        infile.close();
                }
				// 
				Number mean2 = S2.mean();
				Number min2 = S2.min();
				int n2 = S2.length();
				Number max2 = S2.mean();
				arr2 = new Number[n2+1];
				infile.open(a);
				n2 = 1;
				while(infile >> data)
				{
					arr2[n2] = data;
					n2 +=1;
				}
				n2 = n2-1;
				Number adev2 = S2.adev(mean2, arr2);
				Number sdev2 = S2.sdev(mean2, arr2);
				Number skew2 = S2.skew(mean2,sdev2,arr2);
				Number kurt2 = S2.kurt(mean2,sdev2,arr2);
				cout    << "\n\tAbsolute deviation: " << adev2
                     << "\tStandard deviation: " << sdev2
                     << "\n\tSkewness: " << skew2
                         << "\tKurtosis: " << kurt2;

				cout << "\n\n\tLinear correlation (Pearson's r): "
						<< S2.corr(mean,mean2,arr,arr2,sdev,sdev2);

				cout << endl;
				infile.close();
		   }
      
//		   S.reset(); S2.reset;
           cout << "\nEnter the name of the next data file or Q to quit: ";
           cin.getline(a,20);   		
		}
}
/*****************************************************************************
 Given an array of data[1..n], this routine returns its mean ave, average 
 deviation adev, standard deviation sdev, variance var, skewness skew, and 
 kurtosis curt.
******************************************************************************/
void moment(Number data[], int n, Number *ave, Number *adev, Number *sdev,
	Number *var, Number *skew, Number *curt)
{
	int j;
	Number ep=0.0,s,p;

	if (n <= 1) cerr << "n must be at least 2 in moment";
	s=0.0;					/* First pass to get the mean. */
	for (j=1;j<=n;j++) s += data[j];
	*ave=s/n;
	*adev=(*var)=(*skew)=(*curt)=0.0; // Second pass to get the first (absolute),
	for (j=1;j<=n;j++) {			// second, third and fourth moments of the	
		*adev += fabs(s=data[j]-(*ave));   // deviation from the mean. 
		*var += (p=s*s);
		*skew += (p *= s);
		*curt += (p *= s);
	}
	*adev /= n;
	*var=(*var-ep*ep/n)/(n-1);		// Corrected two-pass formula.
	*sdev=sqrt(*var);		// Put the pieces together according to the
	if (*var) {				// conventional definitions.
		*skew /= (n*(*var)*(*sdev));
		*curt=(*curt)/(n*(*var)*(*var))-3.0;
	} else cerr << "No skew/kurtosis when variance = 0 (in moment)";
}