/*
 *@(#)Fit.java 2.0 2005/04/30
 *
 *清華大學　精密儀器與機械學系
 *范燦升　fancansheng@163.com
 */

package algorithm;

import algorithm.LinearEquation;

/**
 *這個類是數據擬合所需用到的類，包括一元線性回歸、多項式回歸、多元線性回歸、三次樣條插值等。
 *@version	2.0, 2005/04/30
 *@author	范燦升
 *@see	Modeling
 *@see	input.FitInput
 */

public class Fit
{
	private double[][] xMatrix;
	private double[] yMatrix;
	private double[] coef;//回歸系數
	private double[] yExpected;//y的期望值
	private double yAverage;
	private int n;//自變量個數
	private int m;//數據組數
	private double u;//回歸平方和
	private double lyy;//總離差平方和
	
	/**
	 *構造擬合與插值類。
	 *@param	xMatrix		自變量數據對應的矩陣，這個矩陣在Fit類中產生克隆副本，Fit的成員方法不會對該矩陣生產影響。
	 *@param	yMatrix		因變量數據對應的矩陣，這個矩陣在Fit類中產生克隆副本，Fit的成員方法不會對該矩陣生產影響。
	 */
	public Fit(double[][] xMatrix,double[] yMatrix)
	{
		this.xMatrix=(double[][])xMatrix.clone();
		this.yMatrix=(double[])yMatrix.clone();
		m=yMatrix.length;
		n=xMatrix[0].length-1;//xMatrix的第一列為1
	}
	
	/**
	 *多元線性回歸。
	 *@return	double[]類的實例，它的成員域length為n+2，其中前n+1個元素為回歸方程中的系數A0、A1、A2……An。
	 *最后一個元素為相關系數，它越接近1，表示回歸效果越顯著。
	 *<p>如果數據錯誤以致生成的線性方程組無解則返回null。
	 */
	public double[] multiFit()
	{
		double[] result=new double[n+2];
		double[][] equationCoefficients;
		double[] equationConstants;
		int i;
		try
		{
			equationCoefficients=Matrix.product(Matrix.transpose(xMatrix),xMatrix);
			equationConstants=Matrix.product(Matrix.transpose(xMatrix),yMatrix);
			coef=new LinearEquation(equationCoefficients,equationConstants).solve();
			for(i=0;i<coef.length;i++)
				result[i]=coef[i];
			result[n+1]=rCoefficent();
			return result;
		}
		catch(NullPointerException e)
		{
			return null;
		}
		//矩陣變換和解方程時可能產生很多null
	}
	
	/**
	 *計算相關系數。
	 *@return	回歸分析的相關系數
	 */
	public double rCoefficent()
	{
		yAverage=0;
		int i,j;
		yExpected=new double[m];
		for(i=0;i<m;i++)
		{
			yAverage+=yMatrix[i];
			yExpected[i]=0;
		}
		yAverage/=m;
		
		for(i=0;i<m;i++)
		{
			for(j=0;j<n+1;j++)
				yExpected[i]+=
				(coef[j]*
				xMatrix[i][j]);
		}
		//計算y的期望值
		
		u=0;
		lyy=0;
		for(i=0;i<m;i++)
		{
			u+=((yExpected[i]-yAverage)*(yExpected[i]-yAverage));
			lyy+=((yMatrix[i]-yAverage)*(yMatrix[i]-yAverage));
		}
		
		double r=Math.sqrt(u/lyy);
		return r;
	}
}