polyfit.java

一个用于排队系统仿真的开源软件,有非常形象的图象仿真过程!

/**    
  * Copyright (C) 2006, Laboratorio di Valutazione delle Prestazioni - Politecnico di Milano

  * This program is free software; you can redistribute it and/or modify
  * it under the terms of the GNU General Public License as published by
  * the Free Software Foundation; either version 2 of the License, or
  * (at your option) any later version.

  * This program is distributed in the hope that it will be useful,
  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  * GNU General Public License for more details.

  * You should have received a copy of the GNU General Public License
  * along with this program; if not, write to the Free Software
  * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
  */
  
/*
 * PolyFit.java
 *
 * Created on 2 novembre 2002, 20.53
 */

package jmt.engine.dataAnalysis;

import Jama.Matrix;
import Jama.SingularValueDecomposition;

/**
 *  This class performs polynomial fitting on given data: there are two
 *  alghorithms.
 *
 * @author  Federico Granata
 */
public class PolyFit {


	/**
     * fit:
	 * returns the coefficients of the polynomial which best approximates f
	 * (in the mean square error sense) for the given data.
	 *  y = f(x).
	 *  y' = a0 + a1 x + ..... + an x^(order)
	 *
	 *  @param  ndata   number of data;
	 *  @param  x       values of the indipendent variable
	 *  @param  y       values of the dipendent variable (data)
	 *  @param  order   the order of the polynomial
	 *
	 *  @return array (of order+1 double) of the coefficients of the polynomial
	 */
	public static double[] fit(int ndata, double[] x, double[] y, int order) {
		double[] a = new double[order + 1];//coefficient of polynomial
		double[][] alfa = new double[order + 1][order + 1];
		double[][] beta = new double[order + 1][1];
		double[] poli = new double[order + 1];//values of y'

		if ((y.length < ndata) || (x.length < ndata))
			return a;

		//create the matrix to solve to find the polynomial coefficient
		poli[0] = 1;
		for (int i = 0; i < ndata; i++) {
			for (int j = 1; j <= order; j++)
				poli[j] = Math.pow(x[i], j);//value of polynomial at x
			for (int k = 0; k <= order; k++) {
				beta[k][0] += y[i] * poli[k];
				for (int j = 0; j <= order; j++)
					alfa[k][j] += poli[k] * poli[j];
			}
		}

		Matrix mAlfa = new Matrix(alfa);
		Matrix mBeta = new Matrix(beta);
		Matrix mA = new Matrix(order + 1, 1);
		mA = mAlfa.solve(mBeta);
		beta = mA.getArrayCopy();
		for (int i = 0; i <= order; i++)
			a[i] = beta[i][0];
		return a;
	}


	/** fitSVD:
	 * returns the coefficients of the polynomial which best approximates f
	 * (in the mean square error sense) for the given data.
	 *  y = f(x).
	 *  y' = a0 + a1 x + ..... + an x^(order)
	 * it uses the Singular Value Decomposition tecnique, very stable.
	 * For a complete description see "Numerical recipes in C"
	 *
	 * @return array (of order+1 double) of the coefficients of the polynomial
	 * @param ndata number of data
	 * @param x values of the independent variable
	 * @param y values of the dependent variable (data)
	 * @param sd the standard deviation of dependent point
	 * @param order the order of the polynomial
	 */
	public static double[] fitSVD(int ndata, double[] x, double[] y, double[] sd, int order) {
		if ((y.length < ndata) || (x.length < ndata) || ((sd != null) && (sd.length < ndata)))
			throw new IllegalArgumentException("data paramenters incorrect too short vector");
		Matrix A = new Matrix(ndata, order + 1);
		double[] b = new double[ndata];
		SingularValueDecomposition SVD;
		double[] a = new double[order + 1];
		double xVal;
		double yVal;
		double small = 0.0000000001 * ndata;
		double[] coeff = new double[order + 1];
		double[][] U;
		double[][] V;
		double[] w;

		if (sd == null) {
			sd = new double[ndata];
			for (int i = 0; i < ndata; i++)
				sd[i] = 1;
		}

		for (int i = 0; i < ndata; i++)
			A.set(i, 0, 1 / sd[i]);
		for (int i = 0; i < ndata; i++) {
			xVal = 1;
			yVal = 1;
			for (int j = 1; j <= order; j++) {
				xVal *= x[i];
				yVal += xVal;
				A.set(i, j, xVal / sd[i]);
			}
			b[i] = yVal / sd[i];
		}
		SVD = A.svd();
		U = (SVD.getU()).getArray();
		V = (SVD.getV()).getArray();
		w = SVD.getSingularValues();

		for (int j = 0; j <= order; j++) {
			coeff[j] = 0;
			if (w[j] > small) {
				for (int i = 0; i < ndata; i++)
					coeff[j] += U[i][j] * b[i];
				coeff[j] /= w[j];
			}
		}

		for (int j = 0; j <= order; j++) {
			a[j] = 0;
			for (int k = 0; k <= order; k++)
				a[j] += coeff[k] * V[k][j];
		}

		return a;
	}

	

	/* writes some data in D:/java/ParameterAnalysis/provaPoly.dat
	 *  then print some polynomial fit coefficients, You have to test them
	 *  against matlab
	 *  @param  args the command line arguments
	 */

    /** tests the fitSVD method
	 */
	public static void testSVD() {
		//Log l = new Log("D:/java/ParameterAnalysis/provaPoly.dat");
		double[] y = new double[200];
		double[] x = new double[200];
		int order;

		for (int i = 0; i < 200; i++) {
			x[i] = i + 1;
			y[i] = Math.random();
			//System.out.println(y[i]+" ");
			//l.write(pippo);
		}
		System.out.println();
		for (order = 0; order < 10; order++) {
			double[] a;
			a = fitSVD(200, x, y, null, order);
			for (int i = 0; i <= order; i++)
				System.out.print(a[i] + "  ");
			System.out.println();
			a = fit(200, x, y, order);
			for (int i = 0; i <= order; i++)
				System.out.print(a[i] + "  ");
			System.out.println();
		}
		//l.close();
	}



}