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
EXAMPLE:
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 six.moves import range
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: norm_factor = sum([exp(-x^2/(2*sigma^2)) for x in range(-6*sigma,sigma*6+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: from six.moves import range
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: from six.moves import range
sage: n=100000; l = [D() for _ in range(n)] # long time
sage: mean(l).n() # long time
0.486650000000000
REFERENCES:
[DDLL13] | (1, 2, 3) Léo Ducas, Alain Durmus, Tancrède Lepoint and Vadim Lyubashevsky. Lattice Signatures and Bimodal Gaussians; in Advances in Cryptology – CRYPTO 2013; Lecture Notes in Computer Science Volume 8042, 2013, pp 40-56 http://www.di.ens.fr/~lyubash/papers/bimodal.pdf |
-
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 ofc
does not have to be an integer. However, some algorithms only support integer-valuedc
(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 [DDLL13] 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 [DDLL13] 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 [DDLL13] 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
TESTS:
We are testing invalid inputs:
sage: from sage.stats.distributions.discrete_gaussian_integer import DiscreteGaussianDistributionIntegerSampler sage: DiscreteGaussianDistributionIntegerSampler(-3.0) Traceback (most recent call last): ... ValueError: sigma must be > 0.0 but got -3.000000 sage: DiscreteGaussianDistributionIntegerSampler(3.0, tau=-1) Traceback (most recent call last): ... ValueError: tau must be >= 1 but got -1 sage: DiscreteGaussianDistributionIntegerSampler(3.0, tau=2, algorithm="superfastalgorithmyouneverheardof") Traceback (most recent call last): ... ValueError: Algorithm 'superfastalgorithmyouneverheardof' not supported by class 'DiscreteGaussianDistributionIntegerSampler' sage: DiscreteGaussianDistributionIntegerSampler(3.0, c=1.5, algorithm="sigma2+logtable") Traceback (most recent call last): ... ValueError: algorithm 'uniform+logtable' requires c%1 == 0
We are testing correctness for multi-precision:
sage: from sage.stats.distributions.discrete_gaussian_integer import DiscreteGaussianDistributionIntegerSampler sage: D = DiscreteGaussianDistributionIntegerSampler(1.0, c=0, tau=2) sage: from six.moves import range sage: l = [D() for _ in range(2^16)] sage: min(l) == 0-2*1.0, max(l) == 0+2*1.0, abs(mean(l)) < 0.01 (True, True, True) sage: from sage.stats.distributions.discrete_gaussian_integer import DiscreteGaussianDistributionIntegerSampler sage: D = DiscreteGaussianDistributionIntegerSampler(1.0, c=2.5, tau=2) sage: from six.moves import range sage: l = [D() for _ in range(2^18)] sage: min(l)==2-2*1.0, max(l)==2+2*1.0, mean(l).n() (True, True, 2.45...) sage: from sage.stats.distributions.discrete_gaussian_integer import DiscreteGaussianDistributionIntegerSampler sage: D = DiscreteGaussianDistributionIntegerSampler(1.0, c=2.5, tau=6) sage: from six.moves import range sage: l = [D() for _ in range(2^18)] sage: min(l), max(l), abs(mean(l)-2.5) < 0.01 (-2, 7, True)
We are testing correctness for double precision:
sage: from sage.stats.distributions.discrete_gaussian_integer import DiscreteGaussianDistributionIntegerSampler sage: D = DiscreteGaussianDistributionIntegerSampler(1.0, c=0, tau=2, precision="dp") sage: from six.moves import range sage: l = [D() for _ in range(2^16)] sage: min(l) == 0-2*1.0, max(l) == 0+2*1.0, abs(mean(l)) < 0.05 (True, True, True) sage: from sage.stats.distributions.discrete_gaussian_integer import DiscreteGaussianDistributionIntegerSampler sage: D = DiscreteGaussianDistributionIntegerSampler(1.0, c=2.5, tau=2, precision="dp") sage: from six.moves import range sage: l = [D() for _ in range(2^18)] sage: min(l)==2-2*1.0, max(l)==2+2*1.0, mean(l).n() (True, True, 2.4...) sage: from sage.stats.distributions.discrete_gaussian_integer import DiscreteGaussianDistributionIntegerSampler sage: D = DiscreteGaussianDistributionIntegerSampler(1.0, c=2.5, tau=6, precision="dp") sage: from six.moves import range sage: l = [D() for _ in range(2^18)] sage: min(l)<=-1, max(l)>=6, abs(mean(l)-2.5) < 0.1 (True, True, True) sage: tuple(l.count(i) for i in range(-2,8)) # output random (7, 242, 4519, 34120, 92714, 91700, 33925, 4666, 246, 5)
We plot a histogram:
sage: from sage.stats.distributions.discrete_gaussian_integer import DiscreteGaussianDistributionIntegerSampler sage: D = DiscreteGaussianDistributionIntegerSampler(17.0) sage: S = [D() for _ in range(2^16)] sage: list_plot([(v,S.count(v)) for v in set(S)]) # long time Graphics object consisting of 1 graphics primitive
These generators cache random bits for performance reasons. Hence, resetting the seed of the PRNG might not have the expected outcome. You can flush this cache with
_flush_cache()
:sage: from sage.stats.distributions.discrete_gaussian_integer import DiscreteGaussianDistributionIntegerSampler sage: D = DiscreteGaussianDistributionIntegerSampler(3.0) sage: sage.misc.randstate.set_random_seed(0); D() 3 sage: sage.misc.randstate.set_random_seed(0); D() 3 sage: sage.misc.randstate.set_random_seed(0); D._flush_cache(); D() 3 sage: D = DiscreteGaussianDistributionIntegerSampler(3.0) sage: sage.misc.randstate.set_random_seed(0); D() 3 sage: sage.misc.randstate.set_random_seed(0); D() 3 sage: sage.misc.randstate.set_random_seed(0); D() -3
-
__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
TESTS:
sage: from sage.stats.distributions.discrete_gaussian_integer import DiscreteGaussianDistributionIntegerSampler sage: DiscreteGaussianDistributionIntegerSampler(3.0, algorithm="uniform+logtable", precision="dp")() # random output 13
-
algorithm
¶
-
c
¶
-
sigma
¶
-
tau
¶
-