Java tutorial
/* * BetaDistribution.java * * BEAST: Bayesian Evolutionary Analysis by Sampling Trees * Copyright (C) 2014 BEAST Developers * * BEAST is free software: you can redistribute it and/or modify it * under the terms of the GNU General Public License as published by * the Free Software Foundation, either version 3 of the License, * or (at your option) any later version. * * BEAST is distributed in the hope that it will be useful, but * WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU General Public License for more details. * * You should have received a copy of the GNU General Public License * along with BEAST. If not, see <http://www.gnu.org/licenses/>. */ package beast.math.distributions; import beast.math.UnivariateFunction; import org.apache.commons.math3.special.Beta; import org.apache.commons.math3.special.Gamma; /** * User: dkuh004 * Date: Mar 25, 2011 * Time: 11:32:25 AM */ public class BetaDistribution implements Distribution { // Default inverse cumulative probability accurac public static final double DEFAULT_INVERSE_ABSOLUTE_ACCURACY = 1e-9; // first shape parameter private double alpha; // second shape parameter private double beta; // Normalizing factor used in density computations. updated whenever alpha or beta are changed. private double z; // Inverse cumulative probability accuracy private final double solverAbsoluteAccuracy; /** * This general constructor creates a new beta distribution with a * specified mean and scale * * @param alpha shape parameter * @param beta shape parameter */ public BetaDistribution(double alpha, double beta) { this.alpha = alpha; this.beta = beta; z = Double.NaN; solverAbsoluteAccuracy = DEFAULT_INVERSE_ABSOLUTE_ACCURACY; } public double getAlpha() { return alpha; } public double getBeta() { return beta; } // Recompute the normalization factor. private void recomputeZ() { if (Double.isNaN(z)) { z = Gamma.logGamma(alpha) + Gamma.logGamma(beta) - Gamma.logGamma(alpha + beta); } } /** * probability density function of the distribution * * @param x argument * @return pdf value */ public double pdf(double x) { recomputeZ(); if (x < 0 || x > 1) { return 0; } else if (x == 0) { if (alpha < 1) { // AR - throwing exceptions deep in numerical code causes trouble. Catching runtime // exceptions is bad. Better to return NaN and let the calling code deal with it. return Double.NaN; // throw MathRuntimeException.createIllegalArgumentException( // "Cannot compute beta density at 0 when alpha = {0,number}", alpha); } return 0; } else if (x == 1) { if (beta < 1) { // AR - throwing exceptions deep in numerical code causes trouble. Catching runtime // exceptions is bad. Better to return NaN and let the calling code deal with it. return Double.NaN; // throw MathRuntimeException.createIllegalArgumentException( // "Cannot compute beta density at 1 when beta = %.3g", beta); } return 0; } else { double logX = Math.log(x); double log1mX = Math.log1p(-x); return Math.exp((alpha - 1) * logX + (beta - 1) * log1mX - z); } } /** * the natural log of the probability density function of the distribution * * @param x argument * @return log pdf value */ public double logPdf(double x) { recomputeZ(); if (x < 0 || x > 1) { return 0; } else if (x == 0) { if (alpha < 1) { // AR - throwing exceptions deep in numerical code causes trouble. Catching runtime // exceptions is bad. Better to return NaN and let the calling code deal with it. return Double.NaN; // throw MathRuntimeException.createIllegalArgumentException( // "Cannot compute beta density at 0 when alpha = {0,number}", alpha); } return 0; } else if (x == 1) { if (beta < 1) { // AR - throwing exceptions deep in numerical code causes trouble. Catching runtime // exceptions is bad. Better to return NaN and let the calling code deal with it. return Double.NaN; // throw MathRuntimeException.createIllegalArgumentException( // "Cannot compute beta density at 1 when beta = %.3g", beta); } return 0; } else { double logX = Math.log(x); double log1mX = Math.log1p(-x); return (alpha - 1) * logX + (beta - 1) * log1mX - z; } } /** * quantile (inverse cumulative density function) of the distribution * * @param y argument * @return icdf value */ public double quantile(double y) { if (y == 0) { return 0; } else if (y == 1) { return 1; } else { throw new RuntimeException(); // try{ // return super.inverseCumulativeProbability(y); // } catch (MathException e) { //// throw MathRuntimeException.createIllegalArgumentException( // AR - throwing exceptions deep in numerical code causes trouble. Catching runtime // // exceptions is bad. Better to return NaN and let the calling code deal with it. // return Double.NaN; // //// "Couldn't calculate beta quantile for alpha = " + alpha + ", beta = " + beta + ": " +e.getMessage()); // } } } protected double getInitialDomain(double p) { return p; } protected double getDomainLowerBound(double p) { return 0; } protected double getDomainUpperBound(double p) { return 1; } public double cdf(double x) { if (x <= 0) { return 0; } else if (x >= 1) { return 1; } else { return Beta.regularizedBeta(x, alpha, beta); } } public double cumulativeProbability(double x) { if (x <= 0) { return 0; } else if (x >= 1) { return 1; } else { return Beta.regularizedBeta(x, alpha, beta); } } public double cumulativeProbability(double x0, double x1) { return cumulativeProbability(x1) - cumulativeProbability(x0); } //Return the absolute accuracy setting of the solver used to estimate inverse cumulative probabilities. protected double getSolverAbsoluteAccuracy() { return solverAbsoluteAccuracy; } /** * cumulative density function of the distribution * * @param x argument * @return cdf value */ // public double cdf(double x){ // throw new UnsupportedOperationException(); // } /** * mean of the distribution * * @return mean */ public double mean() { return (alpha / (alpha + beta)); } /** * variance of the distribution * * @return variance */ public double variance() { return (alpha * beta) / ((alpha + beta) * (alpha + beta) * (alpha + beta + 1)); } /** * @return a probability density function representing this distribution */ public final UnivariateFunction getProbabilityDensityFunction() { return pdfFunction; } private final UnivariateFunction pdfFunction = new UnivariateFunction() { public final double evaluate(double x) { return pdf(x); } public final double getLowerBound() { return 0.0; } public final double getUpperBound() { return 1.0; } }; }