1 /*
2 * Licensed to the Apache Software Foundation (ASF) under one or more
3 * contributor license agreements. See the NOTICE file distributed with
4 * this work for additional information regarding copyright ownership.
5 * The ASF licenses this file to You under the Apache License, Version 2.0
6 * (the "License"); you may not use this file except in compliance with
7 * the License. You may obtain a copy of the License at
8 *
9 * http://www.apache.org/licenses/LICENSE-2.0
10 *
11 * Unless required by applicable law or agreed to in writing, software
12 * distributed under the License is distributed on an "AS IS" BASIS,
13 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14 * See the License for the specific language governing permissions and
15 * limitations under the License.
16 */
17 package org.apache.commons.rng.sampling.distribution;
18
19 import org.apache.commons.rng.UniformRandomProvider;
20 import org.apache.commons.rng.core.source64.SplitMix64;
21 import org.apache.commons.rng.sampling.RandomAssert;
22 import org.junit.jupiter.api.Assertions;
23 import org.junit.jupiter.api.Test;
24
25 /**
26 * Test for the {@link LevySampler}.
27 */
28 class LevySamplerTest {
29 /**
30 * Test the constructor with a negative scale.
31 */
32 @Test
33 void testConstructorThrowsWithNegativeScale() {
34 final UniformRandomProvider rng = RandomAssert.seededRNG();
35 final double location = 1;
36 final double scale = -1e-6;
37 Assertions.assertThrows(IllegalArgumentException.class,
38 () -> LevySampler.of(rng, location, scale));
39 }
40
41 /**
42 * Test the constructor with a zero scale.
43 */
44 @Test
45 void testConstructorThrowsWithZeroScale() {
46 final UniformRandomProvider rng = RandomAssert.seededRNG();
47 final double location = 1;
48 final double scale = 0;
49 Assertions.assertThrows(IllegalArgumentException.class,
50 () -> LevySampler.of(rng, location, scale));
51 }
52
53 /**
54 * Test the SharedStateSampler implementation.
55 */
56 @Test
57 void testSharedStateSampler() {
58 final UniformRandomProvider rng1 = RandomAssert.seededRNG();
59 final UniformRandomProvider rng2 = RandomAssert.seededRNG();
60 final double location = 4.56;
61 final double scale = 1.23;
62 final LevySampler sampler1 = LevySampler.of(rng1, location, scale);
63 final LevySampler sampler2 = sampler1.withUniformRandomProvider(rng2);
64 RandomAssert.assertProduceSameSequence(sampler1, sampler2);
65 }
66
67 /**
68 * Test the support of the standard distribution is {@code [0, inf)}.
69 */
70 @Test
71 void testSupport() {
72 final double location = 0.0;
73 final double scale = 1.0;
74 // Force the underlying ZigguratSampler.NormalizedGaussian to create 0
75 final LevySampler s1 = LevySampler.of(
76 new SplitMix64(0L) {
77 @Override
78 public long next() {
79 return 0L;
80 }
81 }, location, scale);
82 Assertions.assertEquals(Double.POSITIVE_INFINITY, s1.sample());
83
84 // Force the underlying ZigguratSampler.NormalizedGaussian to create a large
85 // sample in the tail of the distribution.
86 // The first two -1,-1 values enters the tail of the distribution.
87 // Here an exponential is added to 3.6360066255009455861.
88 // The exponential also requires -1,-1 to recurse. Each recursion adds 7.569274694148063
89 // to the exponential. A value of 0 stops recursion with a sample of 0.
90 // Two exponentials are required: x and y.
91 // The exponential is multiplied by 0.27502700159745347 to create x.
92 // The condition 2y >= x^x must be true to return x.
93 // Create x = 4 * 7.57 and y = 16 * 7.57
94 final long[] sequence = {
95 // Sample the Gaussian tail
96 -1, -1,
97 // Exponential x = 4 * 7.57... * 0.275027001597525
98 -1, -1, -1, -1, -1, -1, -1, -1, 0,
99 // Exponential y = 16 * 7.57...
100 -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
101 -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, 0,
102 };
103 final LevySampler s2 = LevySampler.of(
104 new SplitMix64(0L) {
105 private int i;
106 @Override
107 public long next() {
108 if (i++ < sequence.length) {
109 return sequence[i - 1];
110 }
111 return super.next();
112 }
113 }, location, scale);
114 // The tail of the zigguart should be approximately s=11.963
115 final double s = 4 * 7.569274694148063 * 0.27502700159745347 + 3.6360066255009455861;
116 // expected is 1/s^2 = 0.006987
117 // So the sampler never achieves the lower bound of zero.
118 // It requires an extreme deviate from the Gaussian.
119 final double expected = 1 / (s * s);
120 Assertions.assertEquals(expected, s2.sample());
121 }
122 }