/*
 * Copyright (c) 2009 The openGion Project.
 *
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 *     http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND,
 * either express or implied. See the License for the specific language
 * governing permissions and limitations under the License.
 */
package org.opengion.penguin.math.statistics;

import java.util.Arrays;

import org.apache.commons.math3.stat.regression.OLSMultipleLinearRegression;

/**
 * apache.commons.mathを利用したOLS重回帰計算のクラスです。
 * y = c0 + x1c1 + x2c2 + x3c3 ...の係数を求めます。
 * c0の切片を考慮するかどうかはnoInterceptで決めます。
 *
 */
public class HybsMultiRegression implements HybsRegression {
	private double cnst[];			// 各係数(xの種類＋１になる？）
	private double rsquare;			// 決定係数
	private boolean noIntercept;	//切片を利用するかどうか

	/**
	 * コンストラクタ。
	 * 与えた二次元データを元に重回帰を計算します。
	 * xデータとして二次元配列を与えます。
	 * noInterceptで切片有り無しを選択します。
	 *
	 * @param in_x 説明変数
	 * @param in_y 目的変数
	 * @param noIntercept 切片利用有無(trueで利用しない)
	 */
	public HybsMultiRegression( final double[][] in_x, final double[] in_y, final boolean noIntercept ) {
		train( in_x, in_y, noIntercept );
	}

	/**
	 * 与えた二次元データを元に重回帰を計算します。
	 * xデータとして二次元配列を与えます。
	 * noInterceptで切片有り無しを選択します。
	 *
	 * @param in_x 説明変数
	 * @param in_y 目的変数
	 * @param noIntercept 切片利用有無(trueで利用しない)
	 */
	private void train( final double[][] in_x, final double[] in_y, final boolean noIntercept ) {
		this.noIntercept = noIntercept;

		// ここで重回帰計算
		final OLSMultipleLinearRegression regression = new OLSMultipleLinearRegression();
		regression.setNoIntercept(noIntercept);
        regression.newSampleData(in_y, in_x);

		cnst	= regression.estimateRegressionParameters();
		rsquare = regression.calculateRSquared();
	}

	/**
	 * 係数をセットした配列を返します。
	 *
	 * @return 係数の配列
	 */
	@Override
	public double[] getCoefficient() {
		return Arrays.copyOf( cnst,cnst.length );
	}

	/**
	 * 決定係数の取得。
	 * @return 決定係数
	 */
	@Override
	public double getRSquare() {
		return rsquare;
	}

	/**
	 * 計算( c0 + c1x1...)を行う。
	 * noInterceptによってc0の利用を決める。
	 * xの大きさが足りない場合は0を返す。
	 *
	 * @param in_x 必要な大きさの変数配列
	 * @return 計算結果
	 */
	@Override
	public double predict( final double... in_x ) {
		double rtn = 0;
		final int itr = noIntercept ? 0 : 1;
		if( in_x.length < cnst.length-itr ) {
			return rtn;
		}

		for( int i=0; i < in_x.length; i++ ) {
			rtn = rtn + in_x[i] * cnst[i+itr];
		}
		if( !noIntercept ) { rtn = rtn + cnst[0]; }

		return rtn;
	}

	//************** ここまでが本体 **************
	/**
	 * ここからテスト用mainメソッド 。
	 *
	 * @param args 引数
	 */
	public static void main( final String[] args ) {
		// データはhttp://mjin.doshisha.ac.jp/R/14.htmlより
		final double[] y = new double[] { 50, 60, 65, 65, 70, 75, 80, 85, 90, 95 };
		double[][] x = new double[10][];
		x[0] = new double[] { 165, 65 };
		x[1] = new double[] { 170, 68 };
		x[2] = new double[] { 172, 70 };
		x[3] = new double[] { 175, 65 };
		x[4] = new double[] { 170, 80 };
		x[5] = new double[] { 172, 85 };
		x[6] = new double[] { 183, 78 };
		x[7] = new double[] { 187, 79 };
		x[8] = new double[] { 180, 95 };
		x[9] = new double[] { 185, 97 };

		final HybsMultiRegression mr = new HybsMultiRegression(x,y,true);

		System.out.println( mr.getRSquare() );
		System.out.println( Arrays.toString( mr.getCoefficient()) );

		System.out.println( mr.predict( new double[] { 169,85 } ));
	}
}

