Skip to content

Commit 4616070

Browse files
cicirellocicirello
committed
Added RandomVariates.nextCauchy()
1 parent a840567 commit 4616070

File tree

3 files changed

+243
-0
lines changed

3 files changed

+243
-0
lines changed

‎lib/jpt2.jar‎

212 Bytes
Binary file not shown.

‎src/org/cicirello/math/rand/RandomVariates.java‎

Lines changed: 113 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -77,4 +77,117 @@ public static int nextBinomial(int n, double p, Random r) {
7777
public static int nextBinomial(int n, double p, SplittableRandom r) {
7878
return BTPE.nextBinomial(n, p, r);
7979
}
80+
81+
/**
82+
* Generates a pseudorandom number from a Cauchy distribution.
83+
* @param median The median of the Cauchy.
84+
* @param scale The scale parameter of the Cauchy.
85+
* @return a pseudorandom number from a Cauchy distribution
86+
*/
87+
public static double nextCauchy(double median, double scale) {
88+
// Inverse Method:
89+
// I'm tempted to map u==0 to Double.NEGATIVE_INFINITY, and u==1
90+
// to Double.POSTIVE_INFINITY (where u=r.nextDouble()).
91+
// However, u will never equal 1.0 since nextDouble's
92+
// return is exclusive of 1.0, and it doesn't seem right to consider the special case
93+
// on the one side but not the other. The maximum absolute value that the call to tan
94+
// below gives appears to be on the order of 1.6ish * 10 to the power 16, which is
95+
// well off from the min/max double values.
96+
return median + scale * StrictMath.tan(StrictMath.PI * (ThreadLocalRandom.current().nextDouble() - 0.5));
97+
}
98+
99+
/**
100+
* Generates a pseudorandom number from a Cauchy distribution with median 0
101+
* and chosen scale parameter.
102+
* @param scale The scale parameter of the Cauchy.
103+
* @return a pseudorandom number from a Cauchy distribution
104+
*/
105+
public static double nextCauchy(double scale) {
106+
// Inverse Method:
107+
// I'm tempted to map u==0 to Double.NEGATIVE_INFINITY, and u==1
108+
// to Double.POSTIVE_INFINITY (where u=r.nextDouble()).
109+
// However, u will never equal 1.0 since nextDouble's
110+
// return is exclusive of 1.0, and it doesn't seem right to consider the special case
111+
// on the one side but not the other. The maximum absolute value that the call to tan
112+
// below gives appears to be on the order of 1.6ish * 10 to the power 16, which is
113+
// well off from the min/max double values.
114+
return scale * StrictMath.tan(StrictMath.PI * (ThreadLocalRandom.current().nextDouble() - 0.5));
115+
}
116+
117+
/**
118+
* Generates a pseudorandom number from a Cauchy distribution.
119+
* @param median The median of the Cauchy.
120+
* @param scale The scale parameter of the Cauchy.
121+
* @param r The source of randomness.
122+
* @return a pseudorandom number from a Cauchy distribution
123+
*/
124+
public static double nextCauchy(double median, double scale, Random r) {
125+
// Inverse Method:
126+
// I'm tempted to map u==0 to Double.NEGATIVE_INFINITY, and u==1
127+
// to Double.POSTIVE_INFINITY (where u=r.nextDouble()).
128+
// However, u will never equal 1.0 since nextDouble's
129+
// return is exclusive of 1.0, and it doesn't seem right to consider the special case
130+
// on the one side but not the other. The maximum absolute value that the call to tan
131+
// below gives appears to be on the order of 1.6ish * 10 to the power 16, which is
132+
// well off from the min/max double values.
133+
return median + scale * StrictMath.tan(StrictMath.PI * (r.nextDouble() - 0.5));
134+
}
135+
136+
/**
137+
* Generates a pseudorandom number from a Cauchy distribution with median 0
138+
* and chosen scale parameter.
139+
* @param scale The scale parameter of the Cauchy.
140+
* @param r The source of randomness.
141+
* @return a pseudorandom number from a Cauchy distribution
142+
*/
143+
public static double nextCauchy(double scale, Random r) {
144+
// Inverse Method:
145+
// I'm tempted to map u==0 to Double.NEGATIVE_INFINITY, and u==1
146+
// to Double.POSTIVE_INFINITY (where u=r.nextDouble()).
147+
// However, u will never equal 1.0 since nextDouble's
148+
// return is exclusive of 1.0, and it doesn't seem right to consider the special case
149+
// on the one side but not the other. The maximum absolute value that the call to tan
150+
// below gives appears to be on the order of 1.6ish * 10 to the power 16, which is
151+
// well off from the min/max double values.
152+
return scale * StrictMath.tan(StrictMath.PI * (r.nextDouble() - 0.5));
153+
}
154+
155+
/**
156+
* Generates a pseudorandom number from a Cauchy distribution.
157+
* @param median The median of the Cauchy.
158+
* @param scale The scale parameter of the Cauchy.
159+
* @param r The source of randomness.
160+
* @return a pseudorandom number from a Cauchy distribution
161+
*/
162+
public static double nextCauchy(double median, double scale, SplittableRandom r) {
163+
// Inverse Method:
164+
// I'm tempted to map u==0 to Double.NEGATIVE_INFINITY, and u==1
165+
// to Double.POSTIVE_INFINITY (where u=r.nextDouble()).
166+
// However, u will never equal 1.0 since nextDouble's
167+
// return is exclusive of 1.0, and it doesn't seem right to consider the special case
168+
// on the one side but not the other. The maximum absolute value that the call to tan
169+
// below gives appears to be on the order of 1.6ish * 10 to the power 16, which is
170+
// well off from the min/max double values.
171+
return median + scale * StrictMath.tan(StrictMath.PI * (r.nextDouble() - 0.5));
172+
}
173+
174+
/**
175+
* Generates a pseudorandom number from a Cauchy distribution with median 0
176+
* and chosen scale parameter.
177+
* @param scale The scale parameter of the Cauchy.
178+
* @param r The source of randomness.
179+
* @return a pseudorandom number from a Cauchy distribution
180+
*/
181+
public static double nextCauchy(double scale, SplittableRandom r) {
182+
// Inverse Method:
183+
// I'm tempted to map u==0 to Double.NEGATIVE_INFINITY, and u==1
184+
// to Double.POSTIVE_INFINITY (where u=r.nextDouble()).
185+
// However, u will never equal 1.0 since nextDouble's
186+
// return is exclusive of 1.0, and it doesn't seem right to consider the special case
187+
// on the one side but not the other. The maximum absolute value that the call to tan
188+
// below gives appears to be on the order of 1.6ish * 10 to the power 16, which is
189+
// well off from the min/max double values.
190+
return scale * StrictMath.tan(StrictMath.PI * (r.nextDouble() - 0.5));
191+
}
192+
80193
}

‎tests/org/cicirello/math/rand/RandomVariatesTests.java‎

Lines changed: 130 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -152,6 +152,136 @@ public void testNextBinomialThreadLocal() {
152152
}
153153
}
154154

155+
@Test
156+
public void testNextCauchyRandom() {
157+
final int TRIALS = 8000;
158+
Random r = new Random(42);
159+
double[] scale = {1, 2, 4};
160+
final double sqrt2 = Math.sqrt(2);
161+
for (double s : scale) {
162+
double[] boundaries = { -s*(sqrt2+1), -s, -s*(sqrt2-1), 0 , s*(sqrt2-1), s, s*(sqrt2+1)};
163+
int[] buckets = new int[8];
164+
for (int j = 0; j < TRIALS; j++) {
165+
double v = RandomVariates.nextCauchy(s, r);
166+
int i = 0;
167+
for (; i < boundaries.length; i++) {
168+
if (v < boundaries[i]) break;
169+
}
170+
buckets[i]++;
171+
}
172+
double chi = chiSquare(buckets);
173+
// critical value for 8-1=7 degrees of freedom is 14.07 (at the .95 level).
174+
assertTrue("Verifying chi square statistic below .95 critical value: chi="+chi+" crit=14.07", chi <= 14.07);
175+
}
176+
double median = 5.0;
177+
double s = 1.0;
178+
double[] boundaries = { median - s*(sqrt2+1), median - s, median - s*(sqrt2-1), median , median + s*(sqrt2-1), median + s, median + s*(sqrt2+1)};
179+
int[] buckets = new int[8];
180+
for (int j = 0; j < TRIALS; j++) {
181+
double v = RandomVariates.nextCauchy(median, s, r);
182+
int i = 0;
183+
for (; i < boundaries.length; i++) {
184+
if (v < boundaries[i]) break;
185+
}
186+
buckets[i]++;
187+
}
188+
double chi = chiSquare(buckets);
189+
// critical value for 8-1=7 degrees of freedom is 14.07 (at the .95 level).
190+
assertTrue("Verifying chi square statistic below .95 critical value: chi="+chi+" crit=14.07", chi <= 14.07);
191+
}
192+
193+
@Test
194+
public void testNextCauchyThreadLocalRandom() {
195+
// Since we cannot seed ThreadLocalRandom, this test case is
196+
// not 100% replicable. Additionally, we know that this version of
197+
// the method simply calls the version that takes a Random as a parameter.
198+
// Since we did test that version for goodness of fit,
199+
// and since replication is not possible without a seed, we simply verify
200+
// that over a large number of trials that samples are generated in each of the buckets
201+
// used for the goodness of fit tests in the class Random case.
202+
final int TRIALS = 800;
203+
double[] scale = {1, 2, 4};
204+
final double sqrt2 = Math.sqrt(2);
205+
for (double s : scale) {
206+
double[] boundaries = { -s*(sqrt2+1), -s, -s*(sqrt2-1), 0 , s*(sqrt2-1), s, s*(sqrt2+1)};
207+
int[] buckets = new int[8];
208+
for (int j = 0; j < TRIALS; j++) {
209+
double v = RandomVariates.nextCauchy(s);
210+
int i = 0;
211+
for (; i < boundaries.length; i++) {
212+
if (v < boundaries[i]) break;
213+
}
214+
buckets[i]++;
215+
}
216+
for (int i = 0; i < buckets.length; i++) {
217+
assertTrue("verifying at least 1 sample in each bucket, scale="+s + " bucket="+i, buckets[i] > 0);
218+
}
219+
}
220+
double median = 5.0;
221+
double s = 1.0;
222+
double[] boundaries = { median - s*(sqrt2+1), median - s, median - s*(sqrt2-1), median , median + s*(sqrt2-1), median + s, median + s*(sqrt2+1)};
223+
int[] buckets = new int[8];
224+
for (int j = 0; j < TRIALS; j++) {
225+
double v = RandomVariates.nextCauchy(median, s);
226+
int i = 0;
227+
for (; i < boundaries.length; i++) {
228+
if (v < boundaries[i]) break;
229+
}
230+
buckets[i]++;
231+
}
232+
for (int i = 0; i < buckets.length; i++) {
233+
assertTrue("verifying at least 1 sample in each bucket, median=5 bucket="+i, buckets[i] > 0);
234+
}
235+
}
236+
237+
@Test
238+
public void testNextCauchySplittableRandom() {
239+
final int TRIALS = 8000;
240+
SplittableRandom r = new SplittableRandom(53);
241+
double[] scale = {1, 2, 4};
242+
final double sqrt2 = Math.sqrt(2);
243+
for (double s : scale) {
244+
double[] boundaries = { -s*(sqrt2+1), -s, -s*(sqrt2-1), 0 , s*(sqrt2-1), s, s*(sqrt2+1)};
245+
int[] buckets = new int[8];
246+
for (int j = 0; j < TRIALS; j++) {
247+
double v = RandomVariates.nextCauchy(s, r);
248+
int i = 0;
249+
for (; i < boundaries.length; i++) {
250+
if (v < boundaries[i]) break;
251+
}
252+
buckets[i]++;
253+
}
254+
double chi = chiSquare(buckets);
255+
// critical value for 8-1=7 degrees of freedom is 14.07 (at the .95 level).
256+
assertTrue("Verifying chi square statistic below .95 critical value: chi="+chi+" crit=14.07", chi <= 14.07);
257+
}
258+
double median = 5.0;
259+
double s = 1.0;
260+
double[] boundaries = { median - s*(sqrt2+1), median - s, median - s*(sqrt2-1), median , median + s*(sqrt2-1), median + s, median + s*(sqrt2+1)};
261+
int[] buckets = new int[8];
262+
for (int j = 0; j < TRIALS; j++) {
263+
double v = RandomVariates.nextCauchy(median, s, r);
264+
int i = 0;
265+
for (; i < boundaries.length; i++) {
266+
if (v < boundaries[i]) break;
267+
}
268+
buckets[i]++;
269+
}
270+
double chi = chiSquare(buckets);
271+
// critical value for 8-1=7 degrees of freedom is 14.07 (at the .95 level).
272+
assertTrue("Verifying chi square statistic below .95 critical value: chi="+chi+" crit=14.07", chi <= 14.07);
273+
}
274+
275+
private double chiSquare(int[] buckets) {
276+
int n = 0;
277+
int v = 0;
278+
for (int i = 0; i < buckets.length; i++) {
279+
n += buckets[i];
280+
v += (buckets[i]*buckets[i]);
281+
}
282+
return v * buckets.length * 1.0 / n - n;
283+
}
284+
155285
private double chiSquare(int[] buckets, double[] dist) {
156286
double v = 0;
157287
int n = 0;

0 commit comments

Comments
 (0)