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}