Class rename (to conform with existing similar functionality).
Project: http://git-wip-us.apache.org/repos/asf/commons-rng/repo Commit: http://git-wip-us.apache.org/repos/asf/commons-rng/commit/8d216eb2 Tree: http://git-wip-us.apache.org/repos/asf/commons-rng/tree/8d216eb2 Diff: http://git-wip-us.apache.org/repos/asf/commons-rng/diff/8d216eb2 Branch: refs/heads/master Commit: 8d216eb2dcfe697d4781c4d04d3c378766ad0429 Parents: 2baae6c Author: Gilles <er...@apache.org> Authored: Wed Oct 11 10:10:08 2017 +0200 Committer: Gilles <er...@apache.org> Committed: Wed Oct 11 10:10:08 2017 +0200 ---------------------------------------------------------------------- .../distribution/ZigguratGaussianSampler.java | 157 ------------------- .../ZigguratNormalizedGaussianSampler.java | 156 ++++++++++++++++++ .../distribution/ContinuousSamplersList.java | 2 +- 3 files changed, 157 insertions(+), 158 deletions(-) ---------------------------------------------------------------------- http://git-wip-us.apache.org/repos/asf/commons-rng/blob/8d216eb2/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/distribution/ZigguratGaussianSampler.java ---------------------------------------------------------------------- diff --git a/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/distribution/ZigguratGaussianSampler.java b/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/distribution/ZigguratGaussianSampler.java deleted file mode 100644 index a4830c6..0000000 --- a/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/distribution/ZigguratGaussianSampler.java +++ /dev/null @@ -1,157 +0,0 @@ -/* - * Licensed to the Apache Software Foundation (ASF) under one or more - * contributor license agreements. See the NOTICE file distributed with - * this work for additional information regarding copyright ownership. - * The ASF licenses this file to You under the Apache License, Version 2.0 - * (the "License"); you may not use this file except in compliance with - * the License. You may obtain a copy of the License at - * - * http://www.apache.org/licenses/LICENSE-2.0 - * - * Unless required by applicable law or agreed to in writing, software - * distributed under the License is distributed on an "AS IS" BASIS, - * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. - * See the License for the specific language governing permissions and - * limitations under the License. - */ - -package org.apache.commons.rng.sampling.distribution; - -import org.apache.commons.rng.UniformRandomProvider; - -/** - * Gaussian Sampling by - * <a href="https://en.wikipedia.org/wiki/Ziggurat_algorithm";>Ziggurat algorithm</a> - * - * <p>Based on - * "The Ziggurat Method for Generating Random Variables"<br> - * by George Marsaglia and Wai Wan Tsang</p> - * - * @see <a href="http://www.jstatsoft.org/article/view/v005i08/ziggurat.pdf";>Ziggurat Method for Generating Random Variables</a> - * - * @since 1.1 - */ - -public class ZigguratGaussianSampler - extends SamplerBase - implements NormalizedGaussianSampler { - - /** - * Generates values from Gaussian (normal) probability distribution - * It uses two tables, integers KN and reals WN. Some 99% of the time, - * the required x is produced by: - * generate a random 32-bit integer j and let i be the index formed from - * the rightmost 8 bits of j. If j < k_i return x = j * w_i. - */ - - /** An auxiliary table of integers, k_i = 2^32*(x_{i-1}/x_i) */ - private static final int[] KN = new int[128]; - /** The table of doubles, w_i = x_i/2^32 */ - private static final double[] WN = new double[128]; - /** The function values table (normalized gaussian in this implementation) f_i(x_i) = exp(-x_i^2/2) */ - private static final double[] FN = new double[128]; - - /** - * Initialize tables. - */ - static { - /** - * Filling the tables. - * k_0 = 2^32 * r * f(dn) / vn - * k_i = 2^32 * ( x_{i-1} / x_i ) - * w_0 = .5^32 * vn / f(dn) - * w_i = .5^32 * x_i - * where dn - the rightmost x_i - * vn - the area of the rectangle - * f(dn) = exp(-.5 * dn * dn) - */ - final double m = 2147483648.0; // 2^31 - - /* provides z(r) = 0, where z(r): x_255 = r, vn = r*f(r)+integral_r^inf f(x)dx */ - final double vn = 9.91256303526217e-3; - - double dn = 3.442619855899; - double tn = dn; - double e = Math.exp(-.5 * dn * dn); - final double q = vn / e; - - KN[0] = (int) ((dn / q) * m); - KN[1] = 0; - - WN[0] = q / m; - WN[127] = dn / m; - - FN[0] = 1.0; - FN[127] = e; - - for (int i = 126; i >= 1; i--){ - e = Math.exp(-.5 *dn * dn); - dn = Math.sqrt(-2. * Math.log(vn / dn + e)); - KN[i+1] = (int) ((dn / tn) * m); - tn = dn; - FN[i] = e; - WN[i] = dn / m; - } - } - - /** - * @param rng Generator of uniformly distributed random numbers. - */ - public ZigguratGaussianSampler(UniformRandomProvider rng) { - super(rng); - } - - /** {@inheritDoc} */ - @Override - public double sample() { - int j = nextInt(); - int i = j & 127; - return (j < KN[i]) ? j * WN[i] : nfix(j,i); - } - - /** - * Get the value from the tail of the distribution - * @param hz - start random integer - * @param iz - corresponding to hz cell's number - * @return the requested random value - */ - private double nfix(int hz, int iz) { - /* The start of the right tail */ - final double r = 3.442619855899; - - double uni; - double x; - double y; - - while (true) { - uni = .5 + hz * .2328306e-9; - x = hz * WN[iz]; - /* iz==0, handles the base strip */ - if (iz == 0) { - /* return x from the tail */ - do { - y = -Math.log(uni); - x = y * 0.2904764; - uni = .5 + nextInt() * .2328306e-9; - } while (y + y < x * x); - return (hz > 0) ? r + x : -r - x; - } - /* iz>0, handle the wedges of other strips */ - if (FN[iz] + uni * (FN[iz - 1] - FN[iz]) < Math.exp(-.5 * x * x)) { - return x; - } - hz = nextInt(); - iz = hz & 127; - if (Math.abs(hz) < KN[iz]) { - return hz * WN[iz]; - } - } - } - - /** {@inheritDoc} */ - @Override - public String toString() { - return "Ziggurat Gaussian deviate [" + super.toString() + "]"; - } - -} http://git-wip-us.apache.org/repos/asf/commons-rng/blob/8d216eb2/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/distribution/ZigguratNormalizedGaussianSampler.java ---------------------------------------------------------------------- diff --git a/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/distribution/ZigguratNormalizedGaussianSampler.java b/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/distribution/ZigguratNormalizedGaussianSampler.java new file mode 100644 index 0000000..2ed8057 --- /dev/null +++ b/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/distribution/ZigguratNormalizedGaussianSampler.java @@ -0,0 +1,156 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one or more + * contributor license agreements. See the NOTICE file distributed with + * this work for additional information regarding copyright ownership. + * The ASF licenses this file to You under the Apache License, Version 2.0 + * (the "License"); you may not use this file except in compliance with + * the License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +package org.apache.commons.rng.sampling.distribution; + +import org.apache.commons.rng.UniformRandomProvider; + +/** + * Gaussian Sampling by + * <a href="https://en.wikipedia.org/wiki/Ziggurat_algorithm";>Ziggurat algorithm</a> + * + * <p>Based on + * "The Ziggurat Method for Generating Random Variables"<br> + * by George Marsaglia and Wai Wan Tsang</p> + * + * @see <a href="http://www.jstatsoft.org/article/view/v005i08/ziggurat.pdf";>Ziggurat Method for Generating Random Variables</a> + * + * @since 1.1 + */ + +public class ZigguratNormalizedGaussianSampler + extends SamplerBase + implements NormalizedGaussianSampler { + + /** + * Generates values from Gaussian (normal) probability distribution + * It uses two tables, integers KN and reals WN. Some 99% of the time, + * the required x is produced by: + * generate a random 32-bit integer j and let i be the index formed from + * the rightmost 8 bits of j. If j < k_i return x = j * w_i. + */ + + /** An auxiliary table of integers, k_i = 2^32*(x_{i-1}/x_i) */ + private static final int[] KN = new int[128]; + /** The table of doubles, w_i = x_i/2^32 */ + private static final double[] WN = new double[128]; + /** The function values table (normalized gaussian in this implementation) f_i(x_i) = exp(-x_i^2/2) */ + private static final double[] FN = new double[128]; + + /** + * Initialize tables. + */ + static { + /** + * Filling the tables. + * k_0 = 2^32 * r * f(dn) / vn + * k_i = 2^32 * ( x_{i-1} / x_i ) + * w_0 = .5^32 * vn / f(dn) + * w_i = .5^32 * x_i + * where dn - the rightmost x_i + * vn - the area of the rectangle + * f(dn) = exp(-.5 * dn * dn) + */ + final double m = 2147483648.0; // 2^31 + + /* provides z(r) = 0, where z(r): x_255 = r, vn = r*f(r)+integral_r^inf f(x)dx */ + final double vn = 9.91256303526217e-3; + + double dn = 3.442619855899; + double tn = dn; + double e = Math.exp(-.5 * dn * dn); + final double q = vn / e; + + KN[0] = (int) ((dn / q) * m); + KN[1] = 0; + + WN[0] = q / m; + WN[127] = dn / m; + + FN[0] = 1.0; + FN[127] = e; + + for (int i = 126; i >= 1; i--){ + e = Math.exp(-.5 *dn * dn); + dn = Math.sqrt(-2. * Math.log(vn / dn + e)); + KN[i+1] = (int) ((dn / tn) * m); + tn = dn; + FN[i] = e; + WN[i] = dn / m; + } + } + + /** + * @param rng Generator of uniformly distributed random numbers. + */ + public ZigguratNormalizedGaussianSampler(UniformRandomProvider rng) { + super(rng); + } + + /** {@inheritDoc} */ + @Override + public double sample() { + int j = nextInt(); + int i = j & 127; + return (j < KN[i]) ? j * WN[i] : nfix(j,i); + } + + /** + * Get the value from the tail of the distribution + * @param hz - start random integer + * @param iz - corresponding to hz cell's number + * @return the requested random value + */ + private double nfix(int hz, int iz) { + /* The start of the right tail */ + final double r = 3.442619855899; + + double uni; + double x; + double y; + + while (true) { + uni = .5 + hz * .2328306e-9; + x = hz * WN[iz]; + /* iz==0, handles the base strip */ + if (iz == 0) { + /* return x from the tail */ + do { + y = -Math.log(uni); + x = y * 0.2904764; + uni = .5 + nextInt() * .2328306e-9; + } while (y + y < x * x); + return (hz > 0) ? r + x : -r - x; + } + /* iz>0, handle the wedges of other strips */ + if (FN[iz] + uni * (FN[iz - 1] - FN[iz]) < Math.exp(-.5 * x * x)) { + return x; + } + hz = nextInt(); + iz = hz & 127; + if (Math.abs(hz) < KN[iz]) { + return hz * WN[iz]; + } + } + } + + /** {@inheritDoc} */ + @Override + public String toString() { + return "Ziggurat normalized Gaussian deviate [" + super.toString() + "]"; + } +} http://git-wip-us.apache.org/repos/asf/commons-rng/blob/8d216eb2/commons-rng-sampling/src/test/java/org/apache/commons/rng/sampling/distribution/ContinuousSamplersList.java ---------------------------------------------------------------------- diff --git a/commons-rng-sampling/src/test/java/org/apache/commons/rng/sampling/distribution/ContinuousSamplersList.java b/commons-rng-sampling/src/test/java/org/apache/commons/rng/sampling/distribution/ContinuousSamplersList.java index dfda24e..bf6e925 100644 --- a/commons-rng-sampling/src/test/java/org/apache/commons/rng/sampling/distribution/ContinuousSamplersList.java +++ b/commons-rng-sampling/src/test/java/org/apache/commons/rng/sampling/distribution/ContinuousSamplersList.java @@ -54,7 +54,7 @@ public class ContinuousSamplersList { meanNormal, sigmaNormal)); // Gaussian ("Ziggurat"). add(LIST, new org.apache.commons.math3.distribution.NormalDistribution(meanNormal, sigmaNormal), - new GaussianSampler(new ZigguratGaussianSampler(RandomSource.create(RandomSource.MT)), + new GaussianSampler(new ZigguratNormalizedGaussianSampler(RandomSource.create(RandomSource.MT)), meanNormal, sigmaNormal)); // Beta ("inverse method").
