code: package xyz.lejon.sampling;
  
  // This code is from the beast-mcmc package with the following copyright
  
  /*
   * GammaDistribution.java
   *
   * Copyright (C) 2002-2006 Alexei Drummond and Andrew Rambaut
   *
   * This file is part of BEAST.
   * See the NOTICE file distributed with this work for additional
   * information regarding copyright ownership and licensing.
   *
   * BEAST is free software; you can redistribute it and/or modify
   * it under the terms of the GNU Lesser General Public License as
   * published by the Free Software Foundation; either version 2
   * 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 Lesser General Public License for more details.
   *
   * You should have received a copy of the GNU Lesser General Public
   * License along with BEAST; if not, write to the
   * Free Software Foundation, Inc., 51 Franklin St, Fifth Floor,
   * Boston, MA  02110-1301  USA
   */
  
  import xyz.lejon.math.GammaFunction;
  
  /**
   * inverse gamma distribution.
   * <p/>
   * (Parameters: shape, scale; mean: ??; variance: ??)
   *
   * @author Joseph Heled
   * @version $Id$
   */
  public class InverseGammaDistribution {
  
      private double shape, scale;
  
      private final double factor;
      private final double logFacor;
  
      public InverseGammaDistribution(double shape, double scale) {
          this.shape = shape;
          this.scale = scale;
          this.factor = Math.pow(scale, shape) / Math.exp(GammaFunction.lnGamma(shape));
          this.logFacor = shape * Math.log(scale) - GammaFunction.lnGamma(shape);
      }
  
      public double getShape() {
          return shape;
      }
  
      public void setShape(double value) {
          shape = value;
      }
  
      public double getScale() {
          return scale;
      }
  
      public void setScale(double value) {
          scale = value;
      }
  
      public double pdf(double x) {
          return pdf(x, shape, scale, factor);
      }
  
      public double logPdf(double x) {
          return logPdf(x, shape, scale, logFacor);
      }
  
      public double cdf(double x) {
          return cdf(x, shape, scale);
      }
  
      public double quantile(double y) {
          return quantile(y, shape, scale);
      }
  
      public double mean() {
          return mean(shape, scale);
      }
  
      public double variance() {
          return variance(shape, scale);
      }
  
      public double nextInverseGamma() {
          return nextInverseGamma(shape, scale);
      }
  
  
      /**
       * probability density function of the Gamma distribution
       *
       * @param x     argument
       * @param shape shape parameter
       * @param scale scale parameter
       * @return pdf value
       */
      static double pdf(double x, double shape, double scale, double factor) {
          if (x <= 0)
              throw new IllegalArgumentException();
  
          final double a = Math.exp(-scale/x) / Math.pow(x, shape+1);
  
          return factor * a;
      }
  
      /**
       * the natural log of the probability density function of the distribution
       *
       * @param x     argument
       * @param shape shape parameter
       * @param scale scale parameter
       * @return log pdf value
       */
      static double logPdf(double x, double shape, double scale, double factor) {
  
          if (x <= 0)
              return Double.NEGATIVE_INFINITY;
  
          return  factor - (scale/x + (shape+1)*Math.log(x));
      }
  
      /**
       * cumulative density function of the Gamma distribution
       *
       * @param x     argument
       * @param shape shape parameter
       * @param scale scale parameter
       * @return cdf value
       */
      public static double cdf(double x, double shape, double scale) {        
          if (x == 0.0) {
              return 0.0;
          }
          return GammaFunction.incompleteGammaQ(shape, scale/x);
      }
  
      /**
       * quantile (inverse cumulative density function) of the Gamma distribution
       *
       * @param y     argument
       * @param shape shape parameter
       * @param scale scale parameter
       * @return icdf value
       */
      public static double quantile(double y, double shape, double scale) {
          // this is what R thinks
          //final GammaDistribution g = new GammaDistribution(shape, scale);
          //return 1/g.quantile(1-y);
          return 1/jdistlib.Gamma.quantile(1-y, shape, scale, true, false);
      }
  
      /**
       * mean of the Gamma distribution
       *
       * @param shape shape parameter
       * @param scale scale parameter
       * @return mean
       */
      public static double mean(double shape, double scale) {
          if( shape > 1 ) {
              return scale / (shape - 1);
          }
          return Double.POSITIVE_INFINITY;
      }
  
      /**
       * variance of the Gamma distribution.
       *
       * @param shape shape parameter
       * @param scale scale parameter
       * @return variance
       */
      public static double variance(double shape, double scale) {
           if( shape > 2 ) {
              return scale*scale / ((shape - 1)*(scale-1)*(scale-2));
          }
           return Double.POSITIVE_INFINITY;
      }
  
      /**
       * sample from the Gamma distribution. This could be calculated using
       * quantile, but current algorithm is faster.
       *
       * @param shape shape parameter
       * @param scale scale parameter
       * @return sample
       */
      public static double nextInverseGamma(double shape, double scale) {
          return 1.0 / new GammaDistribution(shape, 1/scale).nextGamma();
          //return 1.0 / Gamma.rgamma(shape,scale);
      }
      
  }