SaddlePointExpansion.java

  1. /*
  2.  * Licensed to the Apache Software Foundation (ASF) under one or more
  3.  * contributor license agreements.  See the NOTICE file distributed with
  4.  * this work for additional information regarding copyright ownership.
  5.  * The ASF licenses this file to You under the Apache License, Version 2.0
  6.  * (the "License"); you may not use this file except in compliance with
  7.  * the License.  You may obtain a copy of the License at
  8.  *
  9.  *      http://www.apache.org/licenses/LICENSE-2.0
  10.  *
  11.  * Unless required by applicable law or agreed to in writing, software
  12.  * distributed under the License is distributed on an "AS IS" BASIS,
  13.  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
  14.  * See the License for the specific language governing permissions and
  15.  * limitations under the License.
  16.  */
  17. package org.apache.commons.math3.distribution;

  18. import org.apache.commons.math3.special.Gamma;
  19. import org.apache.commons.math3.util.FastMath;
  20. import org.apache.commons.math3.util.MathUtils;

  21. /**
  22.  * <p>
  23.  * Utility class used by various distributions to accurately compute their
  24.  * respective probability mass functions. The implementation for this class is
  25.  * based on the Catherine Loader's <a target="_blank"
  26.  * href="http://www.herine.net/stat/software/dbinom.html">dbinom</a> routines.
  27.  * </p>
  28.  * <p>
  29.  * This class is not intended to be called directly.
  30.  * </p>
  31.  * <p>
  32.  * References:
  33.  * <ol>
  34.  * <li>Catherine Loader (2000). "Fast and Accurate Computation of Binomial
  35.  * Probabilities.". <a target="_blank"
  36.  * href="http://www.herine.net/stat/papers/dbinom.pdf">
  37.  * http://www.herine.net/stat/papers/dbinom.pdf</a></li>
  38.  * </ol>
  39.  * </p>
  40.  *
  41.  * @since 2.1
  42.  */
  43. final class SaddlePointExpansion {

  44.     /** 1/2 * log(2 &#960;). */
  45.     private static final double HALF_LOG_2_PI = 0.5 * FastMath.log(MathUtils.TWO_PI);

  46.     /** exact Stirling expansion error for certain values. */
  47.     private static final double[] EXACT_STIRLING_ERRORS = { 0.0, /* 0.0 */
  48.     0.1534264097200273452913848, /* 0.5 */
  49.     0.0810614667953272582196702, /* 1.0 */
  50.     0.0548141210519176538961390, /* 1.5 */
  51.     0.0413406959554092940938221, /* 2.0 */
  52.     0.03316287351993628748511048, /* 2.5 */
  53.     0.02767792568499833914878929, /* 3.0 */
  54.     0.02374616365629749597132920, /* 3.5 */
  55.     0.02079067210376509311152277, /* 4.0 */
  56.     0.01848845053267318523077934, /* 4.5 */
  57.     0.01664469118982119216319487, /* 5.0 */
  58.     0.01513497322191737887351255, /* 5.5 */
  59.     0.01387612882307074799874573, /* 6.0 */
  60.     0.01281046524292022692424986, /* 6.5 */
  61.     0.01189670994589177009505572, /* 7.0 */
  62.     0.01110455975820691732662991, /* 7.5 */
  63.     0.010411265261972096497478567, /* 8.0 */
  64.     0.009799416126158803298389475, /* 8.5 */
  65.     0.009255462182712732917728637, /* 9.0 */
  66.     0.008768700134139385462952823, /* 9.5 */
  67.     0.008330563433362871256469318, /* 10.0 */
  68.     0.007934114564314020547248100, /* 10.5 */
  69.     0.007573675487951840794972024, /* 11.0 */
  70.     0.007244554301320383179543912, /* 11.5 */
  71.     0.006942840107209529865664152, /* 12.0 */
  72.     0.006665247032707682442354394, /* 12.5 */
  73.     0.006408994188004207068439631, /* 13.0 */
  74.     0.006171712263039457647532867, /* 13.5 */
  75.     0.005951370112758847735624416, /* 14.0 */
  76.     0.005746216513010115682023589, /* 14.5 */
  77.     0.005554733551962801371038690 /* 15.0 */
  78.     };

  79.     /**
  80.      * Default constructor.
  81.      */
  82.     private SaddlePointExpansion() {
  83.         super();
  84.     }

  85.     /**
  86.      * Compute the error of Stirling's series at the given value.
  87.      * <p>
  88.      * References:
  89.      * <ol>
  90.      * <li>Eric W. Weisstein. "Stirling's Series." From MathWorld--A Wolfram Web
  91.      * Resource. <a target="_blank"
  92.      * href="http://mathworld.wolfram.com/StirlingsSeries.html">
  93.      * http://mathworld.wolfram.com/StirlingsSeries.html</a></li>
  94.      * </ol>
  95.      * </p>
  96.      *
  97.      * @param z the value.
  98.      * @return the Striling's series error.
  99.      */
  100.     static double getStirlingError(double z) {
  101.         double ret;
  102.         if (z < 15.0) {
  103.             double z2 = 2.0 * z;
  104.             if (FastMath.floor(z2) == z2) {
  105.                 ret = EXACT_STIRLING_ERRORS[(int) z2];
  106.             } else {
  107.                 ret = Gamma.logGamma(z + 1.0) - (z + 0.5) * FastMath.log(z) +
  108.                       z - HALF_LOG_2_PI;
  109.             }
  110.         } else {
  111.             double z2 = z * z;
  112.             ret = (0.083333333333333333333 -
  113.                     (0.00277777777777777777778 -
  114.                             (0.00079365079365079365079365 -
  115.                                     (0.000595238095238095238095238 -
  116.                                             0.0008417508417508417508417508 /
  117.                                             z2) / z2) / z2) / z2) / z;
  118.         }
  119.         return ret;
  120.     }

  121.     /**
  122.      * A part of the deviance portion of the saddle point approximation.
  123.      * <p>
  124.      * References:
  125.      * <ol>
  126.      * <li>Catherine Loader (2000). "Fast and Accurate Computation of Binomial
  127.      * Probabilities.". <a target="_blank"
  128.      * href="http://www.herine.net/stat/papers/dbinom.pdf">
  129.      * http://www.herine.net/stat/papers/dbinom.pdf</a></li>
  130.      * </ol>
  131.      * </p>
  132.      *
  133.      * @param x the x value.
  134.      * @param mu the average.
  135.      * @return a part of the deviance.
  136.      */
  137.     static double getDeviancePart(double x, double mu) {
  138.         double ret;
  139.         if (FastMath.abs(x - mu) < 0.1 * (x + mu)) {
  140.             double d = x - mu;
  141.             double v = d / (x + mu);
  142.             double s1 = v * d;
  143.             double s = Double.NaN;
  144.             double ej = 2.0 * x * v;
  145.             v *= v;
  146.             int j = 1;
  147.             while (s1 != s) {
  148.                 s = s1;
  149.                 ej *= v;
  150.                 s1 = s + ej / ((j * 2) + 1);
  151.                 ++j;
  152.             }
  153.             ret = s1;
  154.         } else {
  155.             ret = x * FastMath.log(x / mu) + mu - x;
  156.         }
  157.         return ret;
  158.     }

  159.     /**
  160.      * Compute the logarithm of the PMF for a binomial distribution
  161.      * using the saddle point expansion.
  162.      *
  163.      * @param x the value at which the probability is evaluated.
  164.      * @param n the number of trials.
  165.      * @param p the probability of success.
  166.      * @param q the probability of failure (1 - p).
  167.      * @return log(p(x)).
  168.      */
  169.     static double logBinomialProbability(int x, int n, double p, double q) {
  170.         double ret;
  171.         if (x == 0) {
  172.             if (p < 0.1) {
  173.                 ret = -getDeviancePart(n, n * q) - n * p;
  174.             } else {
  175.                 ret = n * FastMath.log(q);
  176.             }
  177.         } else if (x == n) {
  178.             if (q < 0.1) {
  179.                 ret = -getDeviancePart(n, n * p) - n * q;
  180.             } else {
  181.                 ret = n * FastMath.log(p);
  182.             }
  183.         } else {
  184.             ret = getStirlingError(n) - getStirlingError(x) -
  185.                   getStirlingError(n - x) - getDeviancePart(x, n * p) -
  186.                   getDeviancePart(n - x, n * q);
  187.             double f = (MathUtils.TWO_PI * x * (n - x)) / n;
  188.             ret = -0.5 * FastMath.log(f) + ret;
  189.         }
  190.         return ret;
  191.     }
  192. }