/*
 * 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;

import net.morilib.awk.stat.special.Beta;

/**
 * F分布です。
 *
 * @author MORIGUCHI, Yuichiro 2013/04/14
 */
public class FDistribution extends AbstractContinuousDistribution {

	private double d1, d2;

	/**
	 * 確率分布を生成します。
	 * 
	 * @param d1 自由度
	 * @param d2 自由度
	 */
	public FDistribution(int d1, int d2) {
		if(d1 <= 0)  throw new IllegalArgumentException();
		if(d2 <= 0)  throw new IllegalArgumentException();
		this.d1 = d1;
		this.d2 = d2;
	}

	/**
	 * 自由度d_1を得ます。
	 * 
	 * @return d_1
	 */
	public int getD1() {
		return (int)d1;
	}

	/**
	 * 自由度d_2を得ます。
	 * 
	 * @return d_2
	 */
	public int getD2() {
		return (int)d2;
	}

	public boolean isInSupport(double n) {
		return n >= 0;
	}

	public double supportInfimum() {
		return 0;
	}

	public double supportSupremum() {
//		return Double.POSITIVE_INFINITY;
		return 1.0e10;
	}

	private double _cdf(double x) {
		double k, r;

		k = (d1 * x) / (d1 * x + d2);
		r = Beta.I(k, d1 / 2.0, d2 / 2.0);
		return r;
	}

	public double cdf(double x) {
		double r, a, b;

		if(x <= 0) {
			return 0;
		} else {
			r = _cdf(x);
			a = _cdf(x - 0.0001);
			b = _cdf(x + 0.0001);
			if(a > r && r < b) {
				return (a + b) / 2;  // bug fix
			} else if(a < r && r > b) {
				return (a + b) / 2;  // bug fix
			} else if(Double.isNaN(r)) {
				return (a + b) / 2;  // bug fix
			} else if(r != 0.0) {
				return r;
			} else if(d1 > 2) {
				return x > mode() ? 1.0 : r;
			} else if(d2 > 2) {
				return x > expectedValue() ? 1.0 : r;
			} else {
				return (x > 5.0) ? 1.0 : r;
			}
		}
	}

	public double expectedValue() {
		if(d2 > 2) {
			return (double)d2 / (d2 - 2.0);
		} else {
			return Double.NaN;
		}
	}

	public double kurtosis() {
		double d12m2, d2m4, d2m2, d2m6, d2m8, r;

		if(d2 > 8) {
			d12m2 = d1 + d2 - 2.0;
			d2m2  = d2 - 2.0;
			d2m4  = d2 - 4.0;
			d2m6  = d2 - 6.0;
			d2m8  = d2 - 8.0;
			r  = 12.0;
			r *= (d1 * (5.0 * d2 - 22.0) * d12m2 + d2m4 * d2m2 + d2m2);
			r /= d1 * d2m6 * d2m8 * d12m2;
			return r;
		} else {
			return Double.NaN;
		}
	}

	public double mode() {
		if(d1 > 2) {
			return (d1 - 2.0) / d1 * d2 / (d2 + 2.0);
		} else {
			return Double.NaN;
		}
	}

	public double skewness() {
		double d2m4, d2m6, d12m2, d12m22, r;

		if(d2 > 6) {
			d12m2  = d1 + d2 - 2.0;
			d12m22 = 2 * d1 + d2 - 2.0;
			d2m4   = d2 - 4.0;
			d2m6   = d2 - 6.0;
			r  = d12m22 * Math.sqrt(8.0 * d2m4);
			r /= d2m6 * Math.sqrt(d1 * d12m2);
			return r;
		} else {
			return Double.NaN;
		}
	}

	public double variance() {
		double d12m2, d2m4, d2m2, r;

		if(d2 > 4) {
			d12m2 = d1 + d2 - 2.0;
			d2m2  = d2 - 2.0;
			d2m4  = d2 - 4.0;
			r  = 2.0 * d2 * d2 * d12m2;
			r /= d1 * d2m2 * d2m2 * d2m4;
			return r;
		} else {
			return Double.NaN;
		}
	}

	public double f(double x) {
		double l, d1x, d1xd2;

		if(x <= 0) {
			return 0;
		} else {
			d1x   = d1 * x;
			d1xd2 = d1 * x + d2;
			l  = d1 * Math.log(d1x) + d2 * Math.log(d2);
			l -= (d1 + d2) * Math.log(d1xd2);
			l /= 2.0;
			l -= Math.log(x) + Beta.logBeta(d1 / 2, d2 / 2);
			return Math.exp(l);
		}
	}

}
