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

/**
 * ベータ分布です。
 *
 * @author MORIGUCHI, Yuichiro 2013/04/14
 */
public class BetaDistribution extends AbstractContinuousDistribution {

	private double alpha, beta;

	/**
	 * 確率分布を生成します。
	 * 
	 * @param alpha 形状母数
	 * @param beta  形状母数
	 */
	public BetaDistribution(double alpha, double beta) {
		if(alpha <= 0)  throw new IllegalArgumentException();
		if(beta  <= 0)  throw new IllegalArgumentException();
		this.alpha = alpha;
		this.beta  = beta;
	}

	/**
	 * 形状母数αを得ます。
	 * 
	 * @return α
	 */
	public double getAlpha() {
		return alpha;
	}

	/**
	 * 形状母数βを得ます。
	 * 
	 * @return β
	 */
	public double getBeta() {
		return beta;
	}

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

	public double supportInfimum() {
		return 0;
	}

	public double supportSupremum() {
		return 1;
	}

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

		if(x <= 0) {
			return 0;
		} else if(x >= 1) {
			return 1;
		} else {
			r = Beta.I(x, alpha, beta);
			a = Beta.I(x - 0.000000000000001, alpha, beta);
			b = Beta.I(x + 0.000000000000001, alpha, beta);
			if(Double.isNaN(r) ||
					(a < r && r > b) ||
					(a > r && r < b)) {
				return (a + b) / 2;
			} else {
				return (r == 0.0 && x >= 0.9) ? 1.0 : r;
			}
		}
	}

	public double expectedValue() {
		return alpha / (alpha + beta);
	}

	public double kurtosis() {
		double ab1 = alpha + beta + 1;
		double ab2 = alpha + beta + 2;
		double ab3 = alpha + beta + 3;
		double abs = alpha - beta;
		double x;

		x  = 6 * (abs * abs * ab1 - alpha * beta * ab2);
		x /= alpha * beta * ab2 * ab3;
		return x;
	}

	public double mode() {
		return Double.NaN;
	}

	public double skewness() {
		double ab1 = alpha + beta + 1;
		double ab2 = alpha + beta + 2;
		double bas = beta  - alpha;
		double x;

		x  = 2 * bas * Math.sqrt(ab1);
		x /= ab2 * Math.sqrt(alpha * beta);
		return x;
	}

	public double variance() {
		double ab = alpha + beta;
		double ab1 = alpha + beta + 1;

		return alpha * beta / ab / ab / ab1;
	}

	public double f(double x) {
		double l, z, w;

		if(x <= 0 || x >= 1) {
			return 0;
		} else {
			z  = Math.log(x);  w = Math.log(1 - x);
			l  = (alpha - 1) * z + (beta - 1) * w;
			l -= Beta.logBeta(alpha, beta);
			return Math.exp(l);
		}
	}

}
