Partially-Sampled Exponential Random Numbers
Peter Occil - 13/Jul/2020
Peter Occil - 13/Jul/2020
[SHOWTOGROUPS=4,20]
A Python implementation of partially-sampled exponential random numbers, or e-rands.
This page introduces partially-sampled exponential random numbers. Called e-rands, they represent incomplete numbers whose contents are determined only when necessary, making them have potentially arbitrary precision. This page includes an implementation of e-rands in Python, correctness test results, and an application to weighted random sampling of a data stream of unknown size.
Introduction
This page introduces an implementation of partially-sampled exponential random numbers. Called e-rands in this document, they represent incomplete numbers whose contents are determined only when necessary, making them have potentially arbitrary precision.
Moreover, this document includes methods that operate on e-rands in a way that uses only uniform random bits, and without relying on floating-point arithmetic (except for conversion purposes in the example method exprand). Also, the methods support e-rands with an arbitrary rate parameter (λ) greater than 0.
There are papers that discuss generating exponential random numbers using random bits (Flajolet and Saheb 1982)(1), (Karney 2014)(2), (Devroye and Gravel 2015)(3), (Thomas and Luk 2008)(4), but none I am aware of deal with generating partially-sampled exponential random numbers using an arbitrary rate, not just 1.
About the Exponential Distribution
The exponential distribution takes a parameter λ. Informally speaking, a random number that follows an exponential distribution is the number of units of time between one event and the next, and λ is the expected average number of events per unit of time. Usually, λ is equal to 1.
An exponential random number is commonly generated as follows: -ln(1 - RNDU01()) / lamda, where RNDU01() is a uniform random number in the interval [0, 1). (This particular formula is not robust, though, for reasons that are outside the scope of this document, but see (Pedersen 2018)(5).) This page presents an alternative way to sample exponential random numbers.
Code
The following Python code implements partially-sampled exponential random numbers, called e-rands in this document for convenience (similarly to Karney's "u-rands" for partially-sampled uniform random numbers (Karney 2014)(2)). It implements a way to determine whether one e-rand is less than another, as well as a way to fill an e-rand as necessary to give its fractional part a given number of bits.
It makes use of two observations (based on the parameter λ of the exponential distribution):
To implement these probabilities using just random bits, the code uses two algorithms:
[/SHOWTOGROUPS]
A Python implementation of partially-sampled exponential random numbers, or e-rands.
This page introduces partially-sampled exponential random numbers. Called e-rands, they represent incomplete numbers whose contents are determined only when necessary, making them have potentially arbitrary precision. This page includes an implementation of e-rands in Python, correctness test results, and an application to weighted random sampling of a data stream of unknown size.
Introduction
This page introduces an implementation of partially-sampled exponential random numbers. Called e-rands in this document, they represent incomplete numbers whose contents are determined only when necessary, making them have potentially arbitrary precision.
Moreover, this document includes methods that operate on e-rands in a way that uses only uniform random bits, and without relying on floating-point arithmetic (except for conversion purposes in the example method exprand). Also, the methods support e-rands with an arbitrary rate parameter (λ) greater than 0.
There are papers that discuss generating exponential random numbers using random bits (Flajolet and Saheb 1982)(1), (Karney 2014)(2), (Devroye and Gravel 2015)(3), (Thomas and Luk 2008)(4), but none I am aware of deal with generating partially-sampled exponential random numbers using an arbitrary rate, not just 1.
About the Exponential Distribution
The exponential distribution takes a parameter λ. Informally speaking, a random number that follows an exponential distribution is the number of units of time between one event and the next, and λ is the expected average number of events per unit of time. Usually, λ is equal to 1.
An exponential random number is commonly generated as follows: -ln(1 - RNDU01()) / lamda, where RNDU01() is a uniform random number in the interval [0, 1). (This particular formula is not robust, though, for reasons that are outside the scope of this document, but see (Pedersen 2018)(5).) This page presents an alternative way to sample exponential random numbers.
Code
The following Python code implements partially-sampled exponential random numbers, called e-rands in this document for convenience (similarly to Karney's "u-rands" for partially-sampled uniform random numbers (Karney 2014)(2)). It implements a way to determine whether one e-rand is less than another, as well as a way to fill an e-rand as necessary to give its fractional part a given number of bits.
It makes use of two observations (based on the parameter λ of the exponential distribution):
- While a coin flip with probability of heads of exp(-λ) is heads, the exponential random number is increased by 1.
- If a coin flip with probability of heads of 1/(1+exp(λ/2k)) is heads, the exponential random number is increased by 2-k, where k > 0 is an integer.
To implement these probabilities using just random bits, the code uses two algorithms:
- One to simulate a probability of the form exp(-x/y) (zero_or_one_exp_minus; (Canonne et al. 2020)(6)).
- One to simulate a probability of the form 1/(1+exp(x/(y*pow(2, prec)))) (logisticexp (Morina et al. 2019)(7)).
Code:
import random
def logisticexp(ln, ld, prec):
"'" Returns 1 with probability 1/(1+exp(ln/(ld*2^prec))). """
denom=ld*2**prec
while True:
if random.randint(0,1)==0: return 0
if zero_or_one_exp_minus(ln, denom) == 1: return 1
def exprandnew(lamdanum=1, lamdaden=1):
""" Returns an object to serve as a partially-sampled
exponential random number with the given
rate 'lamdanum'/'lamdaden'. The object is a list of five numbers:
the first is a multiple of 2^X, the second is X, the third is the integer
part (initially -1 to indicate the integer part wasn't sampled yet),
and the fourth and fifth are the lamda parameter's
numerator and denominator, respectively. Default for 'lamdanum'
and 'lamdaden' is 1.
The number created by this method will be "empty"
(no bits sampled yet).
"""
return [0, 0, -1, lamdanum, lamdaden]
def exprandfill(a, bits):
""" Fills the unsampled bits of the given exponential random number
'a' as necessary to make a number whose fractional part
has 'bits' many bits. If the number's fractional part already has
that many bits or more, the number is rounded using the round-to-nearest,
ties to even rounding rule. Returns the resulting number as a
multiple of 2^'bits'. """
# Fill the integer if necessary.
if a[2]==-1:
a[2]=0
while zero_or_one_exp_minus(a[3], a[4]) == 1:
a[2]+=1
if a[1] > bits:
# Shifting bits beyond the first excess bit.
aa = a[0] >> (a[1] - bits - 1)
# Check the excess bit; if odd, round up.
ret=aa >> 1 if (aa & 1) == 0 else (aa >> 1) + 1
return ret|(a[2]<<bits)
# Fill the fractional part if necessary.
while a[1] < bits:
index = a[1]
a[1]+=1
a[0]=(a[0]<<1)|logisticexp(a[3], a[4], index+1)
return a[0]|(a[2]<<bits)
def exprandless(a, b):
""" Determines whether one partially-sampled exponential number
is less than another; returns
true if so and false otherwise. During
the comparison, additional bits will be sampled in both numbers
if necessary for the comparison. """
# Check integer part of exponentials
if a[2] == -1:
a[2] = 0
while zero_or_one_exp_minus(a[3], a[4]) == 1:
a[2] += 1
if b[2] == -1:
b[2] = 0
while zero_or_one_exp_minus(b[3], b[4]) == 1:
b[2] += 1
if a[2] < b[2]:
return True
if a[2] > b[2]:
return False
index = 0
while True:
# Fill with next bit in a's exponential number
if a[1] < index:
raise ValueError
if b[1] < index:
raise ValueError
if a[1] <= index:
a[1] += 1
a[0] = logisticexp(a[3], a[4], index + 1) | (a[0] << 1)
# Fill with next bit in b's exponential number
if b[1] <= index:
b[1] += 1
b[0] = logisticexp(b[3], b[4], index + 1) | (b[0] << 1)
aa = (a[0] >> (a[1] - 1 - index)) & 1
bb = (b[0] >> (b[1] - 1 - index)) & 1
if aa < bb:
return True
if aa > bb:
return False
index += 1
def zero_or_one_exp_minus(x, y):
""" Generates 1 with probability exp(-px/py); 0 otherwise.
Reference:
Canonne, C., Kamath, G., Steinke, T., "[The Discrete Gaussian
for Differential Privacy](https://arxiv.org/abs/2004.00010v2)", arXiv:2004.00010v2 [cs.DS], 2020. """
if y <= 0 or x < 0:
raise ValueError
if x > y:
xf = int(x / y) # Get integer part
x = x % y # Reduce to fraction
if x > 0 and zero_or_one_exp_minus(x, y) == 0:
return 0
for i in range(1, xf + 1):
if zero_or_one_exp_minus(1, 1) == 0:
return 0
return 1
r = 1
oy = y
while True:
# NOTE: randint is used in this example, but
# it could be replaced by a random-bit-only
# algorithm such as the Fast Dice Roller (Lumbroso 2013),
# or the Bernoulli method given in the same paper.
if random.randint(0, y-1) >= x:
return r
if r == 1:
r = 0
else:
r = 1
y = y + oy
# Example of use
def exprand(lam):
return exprandfill(exprandnew(lam),53)*1.0/(1<<53)
[/SHOWTOGROUPS]