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: n=100000; l = [D() for _ in range(n)]

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: x=0; l.count(x), ZZ(round(n*exp(-x^2/(2*sigma^2))/norm_factor))
(13355, 13298)
sage: x=4; l.count(x), ZZ(round(n*exp(-x^2/(2*sigma^2))/norm_factor))
(5479, 5467)
sage: x=-10; l.count(x), ZZ(round(n*exp(-x^2/(2*sigma^2))/norm_factor))
(53, 51)

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: n=100000; l = [D() for _ in range(n)] # long time

and check if the proportions fit:

sage: x=0;   y=1; float(l.count(x))/l.count(y), exp(-x^2/(2*sigma^2))/exp(-y^2/(2*sigma^2)).n() # long time
(1.0, 1.00...)
sage: x=0; y=-100; float(l.count(x))/l.count(y), exp(-x^2/(2*sigma^2))/exp(-y^2/(2*sigma^2)).n() # long time
(1.32..., 1.36...)

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: n=100000; l = [D() for _ in range(n)] # long time
sage: mean(l).n() # long time
0.486650000000000

REFERENCES:

class sage.stats.distributions.discrete_gaussian_integer.DiscreteGaussianDistributionIntegerSampler

Bases: sage.structure.sage_object.SageObject

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")()
-3
sage: DiscreteGaussianDistributionIntegerSampler(3.0, algorithm="uniform+table")()
3
algorithm
c
sigma
tau