001/*
002 * Licensed to the Apache Software Foundation (ASF) under one or more
003 * contributor license agreements.  See the NOTICE file distributed with
004 * this work for additional information regarding copyright ownership.
005 * The ASF licenses this file to You under the Apache License, Version 2.0
006 * (the "License"); you may not use this file except in compliance with
007 * the License.  You may obtain a copy of the License at
008 *
009 *      http://www.apache.org/licenses/LICENSE-2.0
010 *
011 * Unless required by applicable law or agreed to in writing, software
012 * distributed under the License is distributed on an "AS IS" BASIS,
013 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
014 * See the License for the specific language governing permissions and
015 * limitations under the License.
016 */
017
018package org.apache.commons.rng.sampling.distribution;
019
020import org.apache.commons.rng.UniformRandomProvider;
021
022/**
023 * Implementation of the <a href="https://en.wikipedia.org/wiki/Zipf's_law">Zipf distribution</a>.
024 *
025 * @since 1.0
026 */
027public class RejectionInversionZipfSampler
028    extends SamplerBase
029    implements DiscreteSampler {
030    /** Number of elements. */
031    private final int numberOfElements;
032    /** Exponent parameter of the distribution. */
033    private final double exponent;
034    /** {@code hIntegral(1.5) - 1}. */
035    private final double hIntegralX1;
036    /** {@code hIntegral(numberOfElements + 0.5)}. */
037    private final double hIntegralNumberOfElements;
038    /** {@code 2 - hIntegralInverse(hIntegral(2.5) - h(2)}. */
039    private final double s;
040
041    /**
042     * @param rng Generator of uniformly distributed random numbers.
043     * @param numberOfElements Number of elements.
044     * @param exponent Exponent.
045     * @throws IllegalArgumentException if {@code numberOfElements <= 0}
046     * or {@code exponent <= 0}.
047     */
048    public RejectionInversionZipfSampler(UniformRandomProvider rng,
049                                         int numberOfElements,
050                                         double exponent) {
051        super(rng);
052        if (numberOfElements <= 0) {
053            throw new IllegalArgumentException("number of elements is not strictly positive: " + numberOfElements);
054        }
055        if (exponent <= 0) {
056            throw new IllegalArgumentException("exponent is not strictly positive: " + exponent);
057        }
058
059        this.numberOfElements = numberOfElements;
060        this.exponent = exponent;
061        this.hIntegralX1 = hIntegral(1.5) - 1;
062        this.hIntegralNumberOfElements = hIntegral(numberOfElements + 0.5);
063        this.s = 2 - hIntegralInverse(hIntegral(2.5) - h(2));
064    }
065
066    /**
067     * Rejection inversion sampling method for a discrete, bounded Zipf
068     * distribution that is based on the method described in
069     * <blockquote>
070     *   Wolfgang Hörmann and Gerhard Derflinger.
071     *   <i>"Rejection-inversion to generate variates from monotone discrete
072     *    distributions",</i><br>
073     *   <strong>ACM Transactions on Modeling and Computer Simulation</strong> (TOMACS) 6.3 (1996): 169-184.
074     * </blockquote>
075     */
076    @Override
077    public int sample() {
078        // The paper describes an algorithm for exponents larger than 1
079        // (Algorithm ZRI).
080        // The original method uses
081        //   H(x) = (v + x)^(1 - q) / (1 - q)
082        // as the integral of the hat function.
083        // This function is undefined for q = 1, which is the reason for
084        // the limitation of the exponent.
085        // If instead the integral function
086        //   H(x) = ((v + x)^(1 - q) - 1) / (1 - q)
087        // is used,
088        // for which a meaningful limit exists for q = 1, the method works
089        // for all positive exponents.
090        // The following implementation uses v = 0 and generates integral
091        // number in the range [1, numberOfElements].
092        // This is different to the original method where v is defined to
093        // be positive and numbers are taken from [0, i_max].
094        // This explains why the implementation looks slightly different.
095
096        while(true) {
097            final double u = hIntegralNumberOfElements + nextDouble() * (hIntegralX1 - hIntegralNumberOfElements);
098            // u is uniformly distributed in (hIntegralX1, hIntegralNumberOfElements]
099
100            double x = hIntegralInverse(u);
101            int k = (int) (x + 0.5);
102
103            // Limit k to the range [1, numberOfElements] if it would be outside
104            // due to numerical inaccuracies.
105            if (k < 1) {
106                k = 1;
107            } else if (k > numberOfElements) {
108                k = numberOfElements;
109            }
110
111            // Here, the distribution of k is given by:
112            //
113            //   P(k = 1) = C * (hIntegral(1.5) - hIntegralX1) = C
114            //   P(k = m) = C * (hIntegral(m + 1/2) - hIntegral(m - 1/2)) for m >= 2
115            //
116            //   where C = 1 / (hIntegralNumberOfElements - hIntegralX1)
117
118            if (k - x <= s || u >= hIntegral(k + 0.5) - h(k)) {
119
120                // Case k = 1:
121                //
122                //   The right inequality is always true, because replacing k by 1 gives
123                //   u >= hIntegral(1.5) - h(1) = hIntegralX1 and u is taken from
124                //   (hIntegralX1, hIntegralNumberOfElements].
125                //
126                //   Therefore, the acceptance rate for k = 1 is P(accepted | k = 1) = 1
127                //   and the probability that 1 is returned as random value is
128                //   P(k = 1 and accepted) = P(accepted | k = 1) * P(k = 1) = C = C / 1^exponent
129                //
130                // Case k >= 2:
131                //
132                //   The left inequality (k - x <= s) is just a short cut
133                //   to avoid the more expensive evaluation of the right inequality
134                //   (u >= hIntegral(k + 0.5) - h(k)) in many cases.
135                //
136                //   If the left inequality is true, the right inequality is also true:
137                //     Theorem 2 in the paper is valid for all positive exponents, because
138                //     the requirements h'(x) = -exponent/x^(exponent + 1) < 0 and
139                //     (-1/hInverse'(x))'' = (1+1/exponent) * x^(1/exponent-1) >= 0
140                //     are both fulfilled.
141                //     Therefore, f(x) = x - hIntegralInverse(hIntegral(x + 0.5) - h(x))
142                //     is a non-decreasing function. If k - x <= s holds,
143                //     k - x <= s + f(k) - f(2) is obviously also true which is equivalent to
144                //     -x <= -hIntegralInverse(hIntegral(k + 0.5) - h(k)),
145                //     -hIntegralInverse(u) <= -hIntegralInverse(hIntegral(k + 0.5) - h(k)),
146                //     and finally u >= hIntegral(k + 0.5) - h(k).
147                //
148                //   Hence, the right inequality determines the acceptance rate:
149                //   P(accepted | k = m) = h(m) / (hIntegrated(m+1/2) - hIntegrated(m-1/2))
150                //   The probability that m is returned is given by
151                //   P(k = m and accepted) = P(accepted | k = m) * P(k = m) = C * h(m) = C / m^exponent.
152                //
153                // In both cases the probabilities are proportional to the probability mass function
154                // of the Zipf distribution.
155
156                return k;
157            }
158        }
159    }
160
161    /** {@inheritDoc} */
162    @Override
163    public String toString() {
164        return "Rejection inversion Zipf deviate [" + super.toString() + "]";
165    }
166
167    /**
168     * {@code H(x)} is defined as
169     * <ul>
170     *  <li>{@code (x^(1 - exponent) - 1) / (1 - exponent)}, if {@code exponent != 1}</li>
171     *  <li>{@code log(x)}, if {@code exponent == 1}</li>
172     * </ul>
173     * H(x) is an integral function of h(x), the derivative of H(x) is h(x).
174     *
175     * @param x Free parameter.
176     * @return {@code H(x)}.
177     */
178    private double hIntegral(final double x) {
179        final double logX = Math.log(x);
180        return helper2((1 - exponent) * logX) * logX;
181    }
182
183    /**
184     * {@code h(x) = 1 / x^exponent}
185     *
186     * @param x Free parameter.
187     * @return {@code h(x)}.
188     */
189    private double h(final double x) {
190        return Math.exp(-exponent * Math.log(x));
191    }
192
193    /**
194     * The inverse function of {@code H(x)}.
195     *
196     * @param x Free parameter
197     * @return y for which {@code H(y) = x}.
198     */
199    private double hIntegralInverse(final double x) {
200        double t = x * (1 - exponent);
201        if (t < -1) {
202            // Limit value to the range [-1, +inf).
203            // t could be smaller than -1 in some rare cases due to numerical errors.
204            t = -1;
205        }
206        return Math.exp(helper1(t) * x);
207    }
208
209    /**
210     * Helper function that calculates {@code log(1 + x) / x}.
211     * <p>
212     * A Taylor series expansion is used, if x is close to 0.
213     * </p>
214     *
215     * @param x A value larger than or equal to -1.
216     * @return {@code log(1 + x) / x}.
217     */
218    private static double helper1(final double x) {
219        if (Math.abs(x) > 1e-8) {
220            return Math.log1p(x) / x;
221        } else {
222            return 1 - x * (0.5 - x * (0.33333333333333333 - 0.25 * x));
223        }
224    }
225
226    /**
227     * Helper function to calculate {@code (exp(x) - 1) / x}.
228     * <p>
229     * A Taylor series expansion is used, if x is close to 0.
230     * </p>
231     *
232     * @param x Free parameter.
233     * @return {@code (exp(x) - 1) / x} if x is non-zero, or 1 if x = 0.
234     */
235    private static double helper2(final double x) {
236        if (Math.abs(x) > 1e-8) {
237            return Math.expm1(x) / x;
238        } else {
239            return 1 + x * 0.5 * (1 + x * 0.33333333333333333 * (1 + 0.25 * x));
240        }
241    }
242}