ZipfDistribution.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.  */

  18. package org.apache.commons.math3.distribution;

  20. import org.apache.commons.math3.exception.NotStrictlyPositiveException;
  21. import org.apache.commons.math3.exception.util.LocalizedFormats;
  22. import org.apache.commons.math3.random.RandomGenerator;
  23. import org.apache.commons.math3.random.Well19937c;
  24. import org.apache.commons.math3.util.FastMath;

  26. /**
  27.  * Implementation of the Zipf distribution.
  28.  *
  29.  * @see <a href="http://mathworld.wolfram.com/ZipfDistribution.html">Zipf distribution (MathWorld)</a>
  30.  */
  31. public class ZipfDistribution extends AbstractIntegerDistribution {
  32.     /** Serializable version identifier. */
  33.     private static final long serialVersionUID = -140627372283420404L;
  34.     /** Number of elements. */
  35.     private final int numberOfElements;
  36.     /** Exponent parameter of the distribution. */
  37.     private final double exponent;
  38.     /** Cached numerical mean */
  39.     private double numericalMean = Double.NaN;
  40.     /** Whether or not the numerical mean has been calculated */
  41.     private boolean numericalMeanIsCalculated = false;
  42.     /** Cached numerical variance */
  43.     private double numericalVariance = Double.NaN;
  44.     /** Whether or not the numerical variance has been calculated */
  45.     private boolean numericalVarianceIsCalculated = false;

  47.     /**
  48.      * Create a new Zipf distribution with the given number of elements and
  49.      * exponent.
  50.      * <p>
  51.      * <b>Note:</b> this constructor will implicitly create an instance of
  52.      * {@link Well19937c} as random generator to be used for sampling only (see
  53.      * {@link #sample()} and {@link #sample(int)}). In case no sampling is
  54.      * needed for the created distribution, it is advised to pass {@code null}
  55.      * as random generator via the appropriate constructors to avoid the
  56.      * additional initialisation overhead.
  57.      *
  58.      * @param numberOfElements Number of elements.
  59.      * @param exponent Exponent.
  60.      * @exception NotStrictlyPositiveException if {@code numberOfElements <= 0}
  61.      * or {@code exponent <= 0}.
  62.      */
  63.     public ZipfDistribution(final int numberOfElements, final double exponent) {
  64.         this(new Well19937c(), numberOfElements, exponent);
  65.     }

  67.     /**
  68.      * Creates a Zipf distribution.
  69.      *
  70.      * @param rng Random number generator.
  71.      * @param numberOfElements Number of elements.
  72.      * @param exponent Exponent.
  73.      * @exception NotStrictlyPositiveException if {@code numberOfElements <= 0}
  74.      * or {@code exponent <= 0}.
  75.      * @since 3.1
  76.      */
  77.     public ZipfDistribution(RandomGenerator rng,
  78.                             int numberOfElements,
  79.                             double exponent)
  80.         throws NotStrictlyPositiveException {
  81.         super(rng);

  83.         if (numberOfElements <= 0) {
  84.             throw new NotStrictlyPositiveException(LocalizedFormats.DIMENSION,
  85.                                                    numberOfElements);
  86.         }
  87.         if (exponent <= 0) {
  88.             throw new NotStrictlyPositiveException(LocalizedFormats.EXPONENT,
  89.                                                    exponent);
  90.         }

  92.         this.numberOfElements = numberOfElements;
  93.         this.exponent = exponent;
  94.     }

  96.     /**
  97.      * Get the number of elements (e.g. corpus size) for the distribution.
  98.      *
  99.      * @return the number of elements
  100.      */
  101.     public int getNumberOfElements() {
  102.         return numberOfElements;
  103.     }

  105.     /**
  106.      * Get the exponent characterizing the distribution.
  107.      *
  108.      * @return the exponent
  109.      */
  110.     public double getExponent() {
  111.         return exponent;
  112.     }

  114.     /** {@inheritDoc} */
  115.     public double probability(final int x) {
  116.         if (x <= 0 || x > numberOfElements) {
  117.             return 0.0;
  118.         }

  120.         return (1.0 / FastMath.pow(x, exponent)) / generalizedHarmonic(numberOfElements, exponent);
  121.     }

  123.     /** {@inheritDoc} */
  124.     @Override
  125.     public double logProbability(int x) {
  126.         if (x <= 0 || x > numberOfElements) {
  127.             return Double.NEGATIVE_INFINITY;
  128.         }

  130.         return -FastMath.log(x) * exponent - FastMath.log(generalizedHarmonic(numberOfElements, exponent));
  131.     }

  133.     /** {@inheritDoc} */
  134.     public double cumulativeProbability(final int x) {
  135.         if (x <= 0) {
  136.             return 0.0;
  137.         } else if (x >= numberOfElements) {
  138.             return 1.0;
  139.         }

  141.         return generalizedHarmonic(x, exponent) / generalizedHarmonic(numberOfElements, exponent);
  142.     }

  144.     /**
  145.      * {@inheritDoc}
  146.      *
  147.      * For number of elements {@code N} and exponent {@code s}, the mean is
  148.      * {@code Hs1 / Hs}, where
  149.      * <ul>
  150.      *  <li>{@code Hs1 = generalizedHarmonic(N, s - 1)},</li>
  151.      *  <li>{@code Hs = generalizedHarmonic(N, s)}.</li>
  152.      * </ul>
  153.      */
  154.     public double getNumericalMean() {
  155.         if (!numericalMeanIsCalculated) {
  156.             numericalMean = calculateNumericalMean();
  157.             numericalMeanIsCalculated = true;
  158.         }
  159.         return numericalMean;
  160.     }

  162.     /**
  163.      * Used by {@link #getNumericalMean()}.
  164.      *
  165.      * @return the mean of this distribution
  166.      */
  167.     protected double calculateNumericalMean() {
  168.         final int N = getNumberOfElements();
  169.         final double s = getExponent();

  171.         final double Hs1 = generalizedHarmonic(N, s - 1);
  172.         final double Hs = generalizedHarmonic(N, s);

  174.         return Hs1 / Hs;
  175.     }

  177.     /**
  178.      * {@inheritDoc}
  179.      *
  180.      * For number of elements {@code N} and exponent {@code s}, the mean is
  181.      * {@code (Hs2 / Hs) - (Hs1^2 / Hs^2)}, where
  182.      * <ul>
  183.      *  <li>{@code Hs2 = generalizedHarmonic(N, s - 2)},</li>
  184.      *  <li>{@code Hs1 = generalizedHarmonic(N, s - 1)},</li>
  185.      *  <li>{@code Hs = generalizedHarmonic(N, s)}.</li>
  186.      * </ul>
  187.      */
  188.     public double getNumericalVariance() {
  189.         if (!numericalVarianceIsCalculated) {
  190.             numericalVariance = calculateNumericalVariance();
  191.             numericalVarianceIsCalculated = true;
  192.         }
  193.         return numericalVariance;
  194.     }

  196.     /**
  197.      * Used by {@link #getNumericalVariance()}.
  198.      *
  199.      * @return the variance of this distribution
  200.      */
  201.     protected double calculateNumericalVariance() {
  202.         final int N = getNumberOfElements();
  203.         final double s = getExponent();

  205.         final double Hs2 = generalizedHarmonic(N, s - 2);
  206.         final double Hs1 = generalizedHarmonic(N, s - 1);
  207.         final double Hs = generalizedHarmonic(N, s);

  209.         return (Hs2 / Hs) - ((Hs1 * Hs1) / (Hs * Hs));
  210.     }

  212.     /**
  213.      * Calculates the Nth generalized harmonic number. See
  214.      * <a href="http://mathworld.wolfram.com/HarmonicSeries.html">Harmonic
  215.      * Series</a>.
  216.      *
  217.      * @param n Term in the series to calculate (must be larger than 1)
  218.      * @param m Exponent (special case {@code m = 1} is the harmonic series).
  219.      * @return the n<sup>th</sup> generalized harmonic number.
  220.      */
  221.     private double generalizedHarmonic(final int n, final double m) {
  222.         double value = 0;
  223.         for (int k = n; k > 0; --k) {
  224.             value += 1.0 / FastMath.pow(k, m);
  225.         }
  226.         return value;
  227.     }

  229.     /**
  230.      * {@inheritDoc}
  231.      *
  232.      * The lower bound of the support is always 1 no matter the parameters.
  233.      *
  234.      * @return lower bound of the support (always 1)
  235.      */
  236.     public int getSupportLowerBound() {
  237.         return 1;
  238.     }

  240.     /**
  241.      * {@inheritDoc}
  242.      *
  243.      * The upper bound of the support is the number of elements.
  244.      *
  245.      * @return upper bound of the support
  246.      */
  247.     public int getSupportUpperBound() {
  248.         return getNumberOfElements();
  249.     }

  251.     /**
  252.      * {@inheritDoc}
  253.      *
  254.      * The support of this distribution is connected.
  255.      *
  256.      * @return {@code true}
  257.      */
  258.     public boolean isSupportConnected() {
  259.         return true;
  260.     }
  261. }
Created with JaCoCo 0.7.2.201409121644