samer@0: /*
samer@0:  *	Copyright (c) 2000, Samer Abdallah, King's College London.
samer@0:  *	All rights reserved.
samer@0:  *
samer@0:  *	This software is provided AS iS and WITHOUT ANY WARRANTY; 
samer@0:  *	without even the implied warranty of MERCHANTABILITY or 
samer@0:  *	FITNESS FOR A PARTICULAR PURPOSE.
samer@0:  */
samer@0: 
samer@0: package samer.models;
samer@0: import  samer.core.*;
samer@0: import  samer.maths.*;
samer@0: import  samer.tools.*;
samer@0: import  java.util.*;
samer@0: 
samer@0: /**
samer@0: 	Collect statistics about a vector: accumulates
samer@0: 	sums and sums of products, then computes mean
samer@0: 	and covariance.
samer@0: 
samer@0: 	This is NOT a model! Computing energies and probabilities
samer@0: 	directly from the first and second moments is horrible, since
samer@0: 	you need the inverse of the covariance matrix. You could,
samer@0: 	however, use this to train a Gaussian model.
samer@0:  */
samer@0: 
samer@0: public class GaussianStats extends NullTask
samer@0: {
samer@0: 	VVector	vector;
samer@0: 	VVector	sumx;
samer@0: 	Matrix		sumxx;
samer@0: 	int			n, count;
samer@0: 	double	[] _sumx, x;
samer@0: 	double	[][] _sumxx;
samer@0: 	double	[] sig;
samer@0: 
samer@0: 	public GaussianStats(VVector vec) throws Exception
samer@0: 	{
samer@0: 		n=vec.size();
samer@0: 		sumx = new VVector("t1",n);	// first moment: sum
samer@0: 		sumxx = new Matrix("t2",n,n);	// second moment: sum products
samer@0: 
samer@0: 		_sumx=sumx.array();
samer@0: 		_sumxx=sumxx.getArray();
samer@0: 		x=vec.array();
samer@0: 		vector=vec;
samer@0: 		sig= new double[n];
samer@0: 		
samer@0: 	}
samer@0: 
samer@0: 	public void dispose() {
samer@0: 		sumx.dispose();
samer@0: 		sumxx.dispose();
samer@0: 	}
samer@0: 
samer@0: 	public void getMean(VVector mu) { getMean(mu.array()); mu.changed(); }
samer@0: 	public void getMean(double [] mu) { Mathx.mul(mu, _sumx, 1.0/count);	}
samer@0: 	public void getCovariance(Matrix C) {
samer@0: 		// compute mean products
samer@0: 		C.assign(sumxx);
samer@0: 		C.timesEquals(1.0/count);
samer@0: 
samer@0: 		// subtract square mean from mean square
samer@0: 		double [][] _C=C.getArray();
samer@0: 
samer@0: 		for (int i=0; i<n; i++) {
samer@0: 			double tmp=_sumx[i]/count; // diagonal
samer@0: 			for (int j=0; j<i; j++) { // off diagonal
samer@0: 				_C[i][j] -= tmp*(_sumx[j]/count);
samer@0: 				_C[j][i] = _C[i][j];
samer@0: 			}
samer@0: 			_C[i][i] -= tmp*tmp;
samer@0: 		}
samer@0: 		C.changed();
samer@0: 	}
samer@0: 
samer@0: 	public void getCorrelation(Matrix R) {
samer@0: 		// compute mean products
samer@0: 		R.assign(sumxx);
samer@0: 		R.timesEquals(1.0/count);
samer@0: 
samer@0: 		double [][] _R=R.getArray();
samer@0: 
samer@0: 		// first get std devs and put 1s down diagonal
samer@0: 		for (int i=0; i<n; i++) {
samer@0: 			double tmp=_sumx[i]/count; // diagonal
samer@0: 			sig[i] = Math.sqrt(_R[i][i] - tmp*tmp);
samer@0: 			_R[i][i]=1;
samer@0: 		}
samer@0: 
samer@0: 		for (int i=0; i<n; i++) {
samer@0: 			double tmp=_sumx[i]/count;
samer@0: 			for (int j=0; j<i; j++) { // off diagonal
samer@0: 				_R[i][j] -= tmp*(_sumx[j]/count);
samer@0: 				_R[i][j] /= sig[i]*sig[j];
samer@0: 				_R[j][i] = _R[i][j];
samer@0: 			}
samer@0: 		}
samer@0: 		R.changed();
samer@0: 	}
samer@0: 
samer@0: 	public void reset() {
samer@0: 		Mathx.zero(_sumx);
samer@0: 		sumxx.zero(); count=0;
samer@0: 	}
samer@0: 
samer@0: 	public void starting() {}
samer@0: 	public void stopping() {}
samer@0: 	public void run()
samer@0: 	{
samer@0: 		Mathx.add(_sumx,x);
samer@0: 
samer@0: 		for (int i=0; i<n; i++) {
samer@0: 			double  [] a=_sumxx[i];
samer@0: 			double  k=x[i];
samer@0: 			for (int j=0; j<=i; j++) {
samer@0: 				a[j] += k*x[j];
samer@0: 			}
samer@0: 		}
samer@0: 		sumx.changed();
samer@0: 		sumxx.changed();
samer@0: 		count++;
samer@0: 	}
samer@0: }
samer@0: