comparison src/samer/models/GaussianStats.java @ 0:bf79fb79ee13

Initial Mercurial check in.

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