/*
 * Copyright 2013 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.dist;


/**
 * 超幾何分布です。
 * 
 *
 * @author MORIGUCHI, Yuichiro 2013/04/20
 */
public class HypergeometricDistribution
extends AbstractDiscreteDistribution {

	private int num, N, K;

	/**
	 * 確率分布を生成します。
	 * 
	 * @param all   母集団の数
	 * @param part  部分集団の数
	 * @param trial 試行回数
	 */
	public HypergeometricDistribution(int all, int part, int trial) {
		this.N   = all;
		this.K   = part;
		this.num = trial;
		if(K   >= N)  throw new IllegalArgumentException();
		if(num >= N)  throw new IllegalArgumentException();
	}

	/**
	 * 母集団の数を得ます。
	 * 
	 * @return 母集団の数
	 */
	public int getAll() {
		return N;
	}

	/**
	 * 部分集団の数を得ます。
	 * 
	 * @return 部分集団
	 */
	public int getPart() {
		return K;
	}

	/**
	 * 試行回数を得ます。
	 * 
	 * @return 試行回数
	 */
	public int getTrial() {
		return num;
	}

	public double f(int n) {
		int k = num, m = 0;
		double r = 1, c = 1;

		if(!isInSupport(n))  return 0;
		for(; k > n; k--, m++) {
			r *= (double)(N   - K - m) / (double)(N - m);
			c *= (double)(num - n - m) / (double)k;
		}

		for(; k > 0; k--, m++) {
			r *= (double)(K + k - n) / (double)(N - m);
		}
		return r / c;
	}

	public boolean isInSupport(int n) {
		return supportMinimum() <= n && n <= supportMaximum();
	}

	public int supportMinimum() {
		return Math.max(0, num + K - N);
	}

	public int supportMaximum() {
		return Math.min(K, num);
	}

	public double expectedValue() {
		return num * K / N;
	}

	public double kurtosis() {
		double NK, N1, N2, N3, Ni, Nn, xd, xn, n = num;

		NK  = N - K;
		N1  = N - 1;  N2 = N - 2;  N3 = N - 3;
		Ni  = N + 1;  Nn = N - n;
		xd  = n * K * Nn * N2 * N3;
		xn  = N1 * N * N * (N * Ni - 6.0 * K * NK - 6.0 * n * Nn);
		xn += 6.0 * n * K * NK * Nn * (5.0 * N - 6.0);
		return xn / xd;
	}

	public double mode() {
		return Math.floor((num + 1.0) * (K + 1.0) / (N + 2.0));
	}

	public double skewness() {
		double NK, N1, N2, Nn, N2n, N2K, xd, xn, n = num;

		NK  = N - K;
		N1  = N - 1;  N2 = N - 2;
		Nn  = N - n;
		N2n = N - 2 * n;
		N2K = N - 2 * K;
		xn  = N2K * Math.sqrt(N1) * N2n;
		xd  = Math.sqrt(n * K * NK * Nn) * N2;
		return xn / xd;
	}

	public double variance() {
		double NK, N1, Nn, n = num;

		NK  = N - K;
		N1  = N - 1;
		Nn  = N - num;
		return n * K / N * NK / N * Nn / N1;
	}

}
