# Discrete Gaussian Samplers over the Integers¶

This class realizes oracles which returns integers proportionally to $$\exp(-(x-c)^2/(2σ^2))$$. All oracles are implemented using rejection sampling. See DiscreteGaussianDistributionIntegerSampler.__init__() for which algorithms are available.

AUTHORS:

• Martin Albrecht (2014-06-28): initial version

EXAMPLES:

We construct a sampler for the distribution $$D_{3,c}$$ with width $$σ=3$$ and center $$c=0$$:

sage: from sage.stats.distributions.discrete_gaussian_integer import DiscreteGaussianDistributionIntegerSampler
sage: sigma = 3.0
sage: D = DiscreteGaussianDistributionIntegerSampler(sigma=sigma)


We ask for 100000 samples:

sage: from collections import defaultdict
sage: counter = defaultdict(Integer)
sage: n = 0
sage: def add_samples(i):
....:     global counter, n
....:     for _ in range(i):
....:         counter[D()] += 1
....:         n += 1

sage: add_samples(100000)


These are sampled with a probability proportional to $$\exp(-x^2/18)$$. More precisely we have to normalise by dividing by the overall probability over all integers. We use the fact that hitting anything more than 6 standard deviations away is very unlikely and compute:

sage: bound = (6*sigma).floor()
sage: norm_factor = sum([exp(-x^2/(2*sigma^2)) for x in range(-bound,bound+1)])
sage: norm_factor
7.519...


With this normalisation factor, we can now test if our samples follow the expected distribution:

sage: expected = lambda x : ZZ(round(n*exp(-x^2/(2*sigma^2))/norm_factor))
sage: observed = lambda x : counter[x]

sage: add_samples(10000)
sage: while abs(observed(0)*1.0/expected(0)   - 1.0) > 5e-2: add_samples(10000)
sage: while abs(observed(4)*1.0/expected(4)   - 1.0) > 5e-2: add_samples(10000)
sage: while abs(observed(-10)*1.0/expected(-10) - 1.0) > 5e-2: add_samples(10000)  # long time


We construct an instance with a larger width:

sage: from sage.stats.distributions.discrete_gaussian_integer import DiscreteGaussianDistributionIntegerSampler
sage: sigma = 127
sage: D = DiscreteGaussianDistributionIntegerSampler(sigma=sigma, algorithm='uniform+online')


ask for 100000 samples:

sage: from collections import defaultdict
sage: counter = defaultdict(Integer)
sage: n = 0
sage: def add_samples(i):
....:     global counter, n
....:     for _ in range(i):
....:         counter[D()] += 1
....:         n += 1

sage: add_samples(100000)


and check if the proportions fit:

sage: expected = lambda x, y: (
....:     exp(-x^2/(2*sigma^2))/exp(-y^2/(2*sigma^2)).n())
sage: observed = lambda x, y: float(counter[x])/counter[y]

sage: while not all(v in counter for v in (0, 1, -100)): add_samples(10000)

sage: while abs(expected(0, 1) - observed(0, 1)) > 2e-1: add_samples(10000)
sage: while abs(expected(0, -100) - observed(0, -100)) > 2e-1: add_samples(10000)


We construct a sampler with $$c\%1 != 0$$:

sage: from sage.stats.distributions.discrete_gaussian_integer import DiscreteGaussianDistributionIntegerSampler
sage: sigma = 3
sage: D = DiscreteGaussianDistributionIntegerSampler(sigma=sigma, c=1/2)
sage: s = 0
sage: n = 0
sage: def add_samples(i):
....:     global s, n
....:     for _ in range(i):
....:         s += D()
....:         n += 1
....:
sage: add_samples(100000)
sage: while abs(float(s)/n - 0.5) > 5e-2: add_samples(10000)


REFERENCES:

class sage.stats.distributions.discrete_gaussian_integer.DiscreteGaussianDistributionIntegerSampler

A Discrete Gaussian Sampler using rejection sampling.

__init__(sigma, c=0, tau=6, algorithm=None, precision='mp')

Construct a new sampler for a discrete Gaussian distribution.

INPUT:

• sigma - samples $$x$$ are accepted with probability proportional to $$\exp(-(x-c)²/(2σ²))$$

• c - the mean of the distribution. The value of c does not have to be an integer. However, some algorithms only support integer-valued c (default: 0)

• tau - samples outside the range $$(⌊c⌉-⌈στ⌉,...,⌊c⌉+⌈στ⌉)$$ are considered to have probability zero. This bound applies to algorithms which sample from the uniform distribution (default: 6)

• algorithm - see list below (default: "uniform+table" for

$$σt$$ bounded by DiscreteGaussianDistributionIntegerSampler.table_cutoff and "uniform+online" for bigger $$στ$$)

• precision - either "mp" for multi-precision where the actual precision used is taken from sigma or "dp" for double precision. In the latter case results are not reproducible. (default: "mp")

ALGORITHMS:

• "uniform+table" - classical rejection sampling, sampling from the uniform distribution and accepted with probability proportional to $$\exp(-(x-c)²/(2σ²))$$ where $$\exp(-(x-c)²/(2σ²))$$ is precomputed and stored in a table. Any real-valued $$c$$ is supported.

• "uniform+logtable" - samples are drawn from a uniform distribution and accepted with probability proportional to $$\exp(-(x-c)²/(2σ²))$$ where $$\exp(-(x-c)²/(2σ²))$$ is computed using logarithmically many calls to Bernoulli distributions. See [DDLL2013] for details. Only integer-valued $$c$$ are supported.

• "uniform+online" - samples are drawn from a uniform distribution and accepted with probability proportional to $$\exp(-(x-c)²/(2σ²))$$ where $$\exp(-(x-c)²/(2σ²))$$ is computed in each invocation. Typically this is very slow. See [DDLL2013] for details. Any real-valued $$c$$ is accepted.

• "sigma2+logtable" - samples are drawn from an easily samplable distribution with $$σ = k·σ_2$$ with $$σ_2 = \sqrt{1/(2\log 2)}$$ and accepted with probability proportional to $$\exp(-(x-c)²/(2σ²))$$ where $$\exp(-(x-c)²/(2σ²))$$ is computed using logarithmically many calls to Bernoulli distributions (but no calls to $$\exp$$). See [DDLL2013] for details. Note that this sampler adjusts $$σ$$ to match $$k·σ_2$$ for some integer $$k$$. Only integer-valued $$c$$ are supported.

EXAMPLES:

sage: from sage.stats.distributions.discrete_gaussian_integer import DiscreteGaussianDistributionIntegerSampler
sage: DiscreteGaussianDistributionIntegerSampler(3.0, algorithm="uniform+online")
Discrete Gaussian sampler over the Integers with sigma = 3.000000 and c = 0
sage: DiscreteGaussianDistributionIntegerSampler(3.0, algorithm="uniform+table")
Discrete Gaussian sampler over the Integers with sigma = 3.000000 and c = 0
sage: DiscreteGaussianDistributionIntegerSampler(3.0, algorithm="uniform+logtable")
Discrete Gaussian sampler over the Integers with sigma = 3.000000 and c = 0


Note that "sigma2+logtable" adjusts $$σ$$:

sage: DiscreteGaussianDistributionIntegerSampler(3.0, algorithm="sigma2+logtable")
Discrete Gaussian sampler over the Integers with sigma = 3.397287 and c = 0

__call__()

Return a new sample.

EXAMPLES:

sage: from sage.stats.distributions.discrete_gaussian_integer import DiscreteGaussianDistributionIntegerSampler
sage: DiscreteGaussianDistributionIntegerSampler(3.0, algorithm="uniform+online")()  # random
-3
sage: DiscreteGaussianDistributionIntegerSampler(3.0, algorithm="uniform+table")()  # random
3

algorithm
c
sigma
tau