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 */ 017package org.apache.commons.rng.sampling.distribution; 018 019import org.apache.commons.rng.UniformRandomProvider; 020 021/** 022 * Sampling from a <a href="https://en.wikipedia.org/wiki/Geometric_distribution">geometric 023 * distribution</a>. 024 * 025 * <p>This distribution samples the number of failures before the first success taking values in the 026 * set {@code [0, 1, 2, ...]}.</p> 027 * 028 * <p>The sample is computed using a related exponential distribution. If \( X \) is an 029 * exponentially distributed random variable with parameter \( \lambda \), then 030 * \( Y = \left \lfloor X \right \rfloor \) is a geometrically distributed random variable with 031 * parameter \( p = 1 − e^\lambda \), with \( p \) the probability of success.</p> 032 * 033 * <p>This sampler outperforms using the {@link InverseTransformDiscreteSampler} with an appropriate 034 * geometric inverse cumulative probability function.</p> 035 * 036 * <p>Usage note: As the probability of success (\( p \)) tends towards zero the mean of the 037 * distribution (\( \frac{1-p}{p} \)) tends towards infinity and due to the use of {@code int} 038 * for the sample this can result in truncation of the distribution.</p> 039 * 040 * <p>Sampling uses {@link UniformRandomProvider#nextDouble()}.</p> 041 * 042 * @see <a 043 * href="https://en.wikipedia.org/wiki/Geometric_distribution#Related_distributions">Geometric 044 * distribution - related distributions</a> 045 * @since 1.3 046 */ 047public final class GeometricSampler { 048 /** 049 * Sample from the geometric distribution when the probability of success is 1. 050 */ 051 private static class GeometricP1Sampler 052 implements SharedStateDiscreteSampler { 053 /** The single instance. */ 054 static final GeometricP1Sampler INSTANCE = new GeometricP1Sampler(); 055 056 @Override 057 public int sample() { 058 // When probability of success is 1 the sample is always zero 059 return 0; 060 } 061 062 @Override 063 public String toString() { 064 return "Geometric(p=1) deviate"; 065 } 066 067 @Override 068 public SharedStateDiscreteSampler withUniformRandomProvider(UniformRandomProvider rng) { 069 // No requirement for a new instance 070 return this; 071 } 072 } 073 074 /** 075 * Sample from the geometric distribution by using a related exponential distribution. 076 */ 077 private static class GeometricExponentialSampler 078 implements SharedStateDiscreteSampler { 079 /** Underlying source of randomness. Used only for the {@link #toString()} method. */ 080 private final UniformRandomProvider rng; 081 /** The related exponential sampler for the geometric distribution. */ 082 private final SharedStateContinuousSampler exponentialSampler; 083 084 /** 085 * @param rng Generator of uniformly distributed random numbers 086 * @param probabilityOfSuccess The probability of success (must be in the range 087 * {@code [0 < probabilityOfSuccess < 1]}) 088 */ 089 GeometricExponentialSampler(UniformRandomProvider rng, double probabilityOfSuccess) { 090 this.rng = rng; 091 // Use a related exponential distribution: 092 // λ = −ln(1 − probabilityOfSuccess) 093 // exponential mean = 1 / λ 094 // -- 095 // Note on validation: 096 // If probabilityOfSuccess == Math.nextDown(1.0) the exponential mean is >0 (valid). 097 // If probabilityOfSuccess == Double.MIN_VALUE the exponential mean is +Infinity 098 // and the sample will always be Integer.MAX_VALUE (the distribution is truncated). It 099 // is noted in the class Javadoc that the use of a small p leads to truncation so 100 // no checks are made for this case. 101 final double exponentialMean = 1.0 / (-Math.log1p(-probabilityOfSuccess)); 102 exponentialSampler = ZigguratSampler.Exponential.of(rng, exponentialMean); 103 } 104 105 /** 106 * @param rng Generator of uniformly distributed random numbers 107 * @param source Source to copy. 108 */ 109 GeometricExponentialSampler(UniformRandomProvider rng, GeometricExponentialSampler source) { 110 this.rng = rng; 111 exponentialSampler = source.exponentialSampler.withUniformRandomProvider(rng); 112 } 113 114 @Override 115 public int sample() { 116 // Return the floor of the exponential sample 117 return (int) Math.floor(exponentialSampler.sample()); 118 } 119 120 @Override 121 public String toString() { 122 return "Geometric deviate [" + rng.toString() + "]"; 123 } 124 125 @Override 126 public SharedStateDiscreteSampler withUniformRandomProvider(UniformRandomProvider rng) { 127 return new GeometricExponentialSampler(rng, this); 128 } 129 } 130 131 /** Class contains only static methods. */ 132 private GeometricSampler() {} 133 134 /** 135 * Creates a new geometric distribution sampler. The samples will be provided in the set 136 * {@code k=[0, 1, 2, ...]} where {@code k} indicates the number of failures before the first 137 * success. 138 * 139 * @param rng Generator of uniformly distributed random numbers. 140 * @param probabilityOfSuccess The probability of success. 141 * @return the sampler 142 * @throws IllegalArgumentException if {@code probabilityOfSuccess} is not in the range 143 * {@code [0 < probabilityOfSuccess <= 1]}) 144 */ 145 public static SharedStateDiscreteSampler of(UniformRandomProvider rng, 146 double probabilityOfSuccess) { 147 if (probabilityOfSuccess <= 0 || probabilityOfSuccess > 1) { 148 throw new IllegalArgumentException( 149 "Probability of success (p) must be in the range [0 < p <= 1]: " + 150 probabilityOfSuccess); 151 } 152 return probabilityOfSuccess == 1 ? 153 GeometricP1Sampler.INSTANCE : 154 new GeometricExponentialSampler(rng, probabilityOfSuccess); 155 } 156}