/*
 * Copyright 2009-2010 Yuichiro Moriguchi
 *
 * 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 net.morilib.awk.stat;

import java.math.BigInteger;
import java.util.Arrays;

import net.morilib.awk.stat.dist.GammaDistribution;

/**
 * 統計関数に関するユーティリティメソッドです。
 * 
 *
 * @author MORIGUCHI, Yuichiro 2012/02/17
 */
public final class StatisticsUtils {

	private static final BigInteger TWO = BigInteger.valueOf(2);

	//
	private StatisticsUtils() {}

	/**
	 * 和を計算します。
	 * 
	 * @param vals 値
	 * @return 和
	 */
	public static double sum(double... vals) {
		double r = 0.0;

		for(int i = 0; i < vals.length; i++) {
			r += vals[i];
		}
		return r;
	}

	/**
	 * 和を計算します。
	 * 
	 * @param vals 値
	 * @return 和
	 */
	public static long sum(int... vals) {
		long r = 0;

		for(int i = 0; i < vals.length; i++) {
			r += vals[i];
		}
		return r;
	}

	/**
	 * 算術平均を計算します。
	 * 
	 * @param vals 値
	 * @return 算術平均
	 */
	public static double mean(double... vals) {
		return sum(vals) / vals.length;
	}

	/**
	 * 幾何平均を計算します。
	 * 
	 * @param vals 値
	 * @return 幾何平均
	 */
	public static double geometricMean(double... vals) {
		double r = 0.0;

		for(int i = 0; i < vals.length; i++) {
			r += Math.log(vals[i]);
		}
		return Math.exp(r / vals.length);
	}

	/**
	 * 調和平均を計算します。
	 * 
	 * @param vals 値
	 * @return 調和平均
	 */
	public static double harmonicMean(double... vals) {
		double r = 0.0;

		for(int i = 0; i < vals.length; i++) {
			r += 1.0 / vals[i];
		}
		return vals.length / r;
	}

	/**
	 * メディアンを計算します。
	 * 
	 * @param vals 値
	 * @return メディアン
	 */
	public static double median(double... vals) {
		double[] ta = new double[vals.length];

		System.arraycopy(vals, 0, ta, 0, vals.length);
		Arrays.sort(ta);
		if(vals.length % 2 == 0) {
			return (ta[vals.length / 2 - 1] + ta[vals.length / 2]) / 2;
		} else {
			return ta[vals.length / 2];
		}
	}

	/**
	 * メディアンを計算します。
	 * 
	 * @param vals 値
	 * @return メディアン
	 */
	public static BigInteger median(BigInteger... vals) {
		BigInteger[] ta = new BigInteger[vals.length];

		System.arraycopy(vals, 0, ta, 0, vals.length);
		Arrays.sort(ta);
		if(vals.length % 2 == 0) {
			return (ta[vals.length / 2 - 1].add(
					ta[vals.length / 2])).divide(TWO);
		} else {
			return ta[vals.length / 2];
		}
	}

	/**
	 * パーセンタイルを計算します。
	 * 
	 * @param percent パーセント
	 * @param vals 値
	 * @return パーセンタイル
	 */
	public static double percentile(double percent, double... vals) {
		double[] ta = new double[vals.length];
		double p, q;
		int i;

		System.arraycopy(vals, 0, ta, 0, vals.length);
		Arrays.sort(ta);
		p = percent / 100 * (vals.length - 1);
		i = (int)p;
		q = p - i;
		if(percent < 0 || percent > 100) {
			throw new IllegalArgumentException();
		} else if(percent == 0) {
			return ta[0];
		} else if(percent == 100) {
			return ta[ta.length - 1];
		} else {
			return ta[i] * (1 - q) + ta[i + 1] * q;
		}
	}

	/**
	 * 最大値を計算します。
	 * 
	 * @param vals 値
	 * @return 最大値
	 */
	public static double max(double... vals) {
		double r = Double.MIN_VALUE;

		for(int i = 0; i < vals.length; i++) {
			if(vals[i] > r)  r = vals[i];
		}
		return r;
	}

	/**
	 * 最大値を計算します。
	 * 
	 * @param vals 値
	 * @return 最大値
	 */
	public static int max(int... vals) {
		int r = Integer.MIN_VALUE;

		for(int i = 0; i < vals.length; i++) {
			if(vals[i] > r)  r = vals[i];
		}
		return r;
	}

	/**
	 * 最大値がどこにあるかを計算します。
	 * 
	 * @param vals 値
	 * @return 最大値の添字
	 */
	public static int whereIsMax(int... vals) {
		int r = Integer.MIN_VALUE, v = -1;

		for(int i = 0; i < vals.length; i++) {
			if(vals[i] > r) {
				r = vals[i];  v = i;
			}
		}
		return v;
	}

	/**
	 * 最小値を計算します。
	 * 
	 * @param vals 値
	 * @return 最大値
	 */
	public static double min(double... vals) {
		double r = Double.MAX_VALUE;

		for(int i = 0; i < vals.length; i++) {
			if(vals[i] < r)  r = vals[i];
		}
		return r;
	}

	/**
	 * 最小値を計算します。
	 * 
	 * @param vals 値
	 * @return 最大値
	 */
	public static int min(int... vals) {
		int r = Integer.MAX_VALUE;

		for(int i = 0; i < vals.length; i++) {
			if(vals[i] < r)  r = vals[i];
		}
		return r;
	}

	/**
	 * モードを計算します。
	 * 
	 * @param vals 値
	 * @return モード
	 */
	public static double mode(double... vals) {
		double max = max(vals), min = min(vals);
		int[] hist = toHistogram(vals);

		return mode(min, max, hist);
	}

	/**
	 * モードを計算します。
	 * 
	 * @param vals 値
	 * @return モード
	 */
	public static BigInteger mode(BigInteger... vals) {
		BigInteger[] ta = new BigInteger[vals.length];
		BigInteger p = null, m = p, v;
		int r = 0, c = 0;

		System.arraycopy(vals, 0, ta, 0, vals.length);
		Arrays.sort(ta);
		for(int i = 0; i < vals.length; i++) {
			v = ta[i];
			if(p == null || p != v) {
				p = v;
				if(c > r) {
					r = c;  m = p;
				}
				c = 1;
			} else {
				c++;
			}
		}
		return m;
	}

	/**
	 * 配列をヒストグラムに変換します。
	 * 
	 * @param vals 値
	 * @return ヒストグラム
	 */
	public static int[] toHistogram(double... vals) {
		int k = (int)Math.sqrt(vals.length);
		int[] r = new int[k + 2];
		double max = max(vals), min = min(vals);
		double h;

		if(vals.length == 0) {
			throw new IllegalArgumentException();
		}

		Arrays.fill(r, 0);
		h = (max - min) / k;
		for(int i = 0; i < vals.length; i++) {
			if(vals[i] == max) {
				r[k]++;
			} else if(vals[i] > max) {
				r[k + 1]++;
			} else if(vals[i] < min) {
				r[0]++;
			} else {
				r[(int)((vals[i] - min) / h) + 1]++;
			}
		}
		return r;
	}

	/**
	 * モードを計算します。
	 * 
	 * @param min モードを計算する下限
	 * @param max モードを計算する上限
	 * @param histogram ヒストグラム
	 * @return モード
	 */
	public static double mode(double min, double max,
			int... histogram) {
		int p = whereIsMax(histogram);
		double q = (double)(p - 0.5) / (histogram.length - 2);

		if(histogram.length < 2) {
			throw new IllegalArgumentException();
		} else if(p == 0) {
			return Double.NEGATIVE_INFINITY;
		} else if(p == histogram.length) {
			return Double.POSITIVE_INFINITY;
		} else {
			return (max - min) * q + min;
		}
	}

	/**
	 * 分散を計算します。
	 * 
	 * @param vals 値
	 * @return 分散
	 */
	public static double variance(double... vals) {
		double r = 0.0, m = mean(vals);

		for(int i = 0; i < vals.length; i++) {
			r += (vals[i] - m) * (vals[i] - m);
		}
		return r / vals.length;
	}

	/**
	 * 不偏分散を計算します。
	 * 
	 * @param vals 値
	 * @return 不偏分散
	 */
	public static double unbiasedVariance(double... vals) {
		double r = 0.0, m = mean(vals);

		for(int i = 0; i < vals.length; i++) {
			r += (vals[i] - m) * (vals[i] - m);
		}
		return r / (vals.length - 1);
	}

	/**
	 * 歪度を計算します。
	 * 
	 * @param vals 値
	 * @return 歪度
	 */
	public static double skewness(double... vals) {
		double r = 0.0, m = mean(vals);
		double x, s = Math.sqrt(variance(vals));

		for(int i = 0; i < vals.length; i++) {
			x  = (vals[i] - m) / s;
			r += x * x * x;
		}
		return r / vals.length;
	}

	/**
	 * 尖度を計算します。
	 * 
	 * @param vals 値
	 * @return 尖度
	 */
	public static double kurtosis(double... vals) {
		double r = 0.0, m = mean(vals), s = variance(vals);

		for(int i = 0; i < vals.length; i++) {
			r += ((vals[i] - m) * (vals[i] - m) *
					(vals[i] - m) * (vals[i] - m));
		}
		return r / vals.length / s / s - 3.0;
	}

	/**
	 * 共分散を計算します。
	 * 
	 * @param x 値1
	 * @param y 値2
	 * @return 共分散
	 */
	public static double covariance(double[] x, double[] y) {
		double r = 0.0, mx = mean(x), my = mean(y);

		if(x == null || y == null) {
			throw new NullPointerException();
		} else if(x.length != y.length) {
			throw new IllegalArgumentException();
		} else if(x.length == 0) {
			throw new IllegalArgumentException();
		}

		for(int i = 0; i < x.length; i++) {
			r += (x[i] - mx) * (y[i] - my);
		}
		return r / x.length;
	}

	/**
	 * 相関係数を計算します。
	 * 
	 * @param x 値1
	 * @param y 値2
	 * @return 相関係数
	 */
	public static double correlation(double[] x, double[] y) {
		double r = 0.0, xs = 0.0, ys = 0.0;
		double mx = mean(x), my = mean(y);

		if(x == null || y == null) {
			throw new NullPointerException();
		} else if(x.length != y.length) {
			throw new IllegalArgumentException();
		} else if(x.length == 0) {
			throw new IllegalArgumentException();
		}

		for(int i = 0; i < x.length; i++) {
			r  += (x[i] - mx) * (y[i] - my);
			xs += (x[i] - mx) * (x[i] - mx);
			ys += (y[i] - my) * (y[i] - my);
		}
		return r / Math.sqrt(xs) / Math.sqrt(ys);
	}

	/**
	 * 2つのデータについて線形回帰をします。<br />
	 * 結果はr[1]*x + r[0]の形で返されます。
	 * 
	 * @param x 値1
	 * @param y 値2
	 * @return 
	 */
	public static double[] regression(double[] x, double[] y) {
		double x11 = 0.0, x12 = 0.0, x21 = 0.0, x22 = 0.0;
		double y1  = 0.0, y2  = 0.0;
		double bt0, bt1, det;

		if(x == null || y == null) {
			throw new NullPointerException();
		} else if(x.length != y.length) {
			throw new IllegalArgumentException();
		} else if(x.length == 0) {
			throw new IllegalArgumentException();
		}

		for(int i = 0; i < x.length; i++) {
			x11 += 1;
			x12 += x[i];
			x21 += x[i];
			x22 += x[i] * x[i];
			y1  += y[i];
			y2  += x[i] * y[i];
		}
		det  = x11 * x22 - x12 * x21;
		bt0  = (x11 *  y2 - x21 *  y1) / det;
		bt1  = (y1  * x22 - y2  * x12) / det;
		return new double[] { bt0, bt1 };
	}

	/**
	 * カイ二乗検定をします。<br />
	 * 帰無仮説が棄却されたときfalseを得ます。
	 * 
	 * @param obs 観測値
	 * @param exp 理論値
	 * @param alpha 有意水準
	 * @return 帰無仮説が棄却されたときfalse
	 */
	public static boolean testChiSquared(double[] obs, double[] exp,
			double alpha) {
		GammaDistribution dis;
		double x, chi = 0;
		int i = 0;

		for(; i < obs.length && i < exp.length; i++) {
			x    = obs[i] - exp[i];
			chi += x * x / exp[i];
		}
		dis = GammaDistribution.chiSquared(i - 1);
		return chi < dis.invCdf(1 - alpha);
	}

}
