Skip to content

Commit 11873c2

Browse files
cicirellocicirello
committed
nextCauchy: optimized and adjusted for bias
1 parent a9ccd75 commit 11873c2

File tree

2 files changed

+174
-48
lines changed

2 files changed

+174
-48
lines changed

‎lib/jpt2.jar‎

123 Bytes
Binary file not shown.

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

Lines changed: 174 additions & 48 deletions
Original file line numberDiff line numberDiff line change
@@ -86,14 +86,35 @@ public static int nextBinomial(int n, double p, SplittableRandom r) {
8686
*/
8787
public static double nextCauchy(double median, double scale) {
8888
// 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));
89+
// Mathematically, it should be: median + scale * tan(PI * (u - 0.5)),
90+
// where u is uniformly random from the interval [0, 1].
91+
// However, since tan goes through one complete cycle every PI,
92+
// we can replace it with: median + scale * tan(PI * u), going from
93+
// 0 to PI, rather than from -PI/2 to PI/2. This is equivalent
94+
// as far as generating a random Cauchy variate is concerned, but saves
95+
// one arithmetic operation.
96+
// At first glance, it may appear as if we might be
97+
// doubly sampling u == 0 since tan(0)==tan(PI), however, our uniform
98+
// random numbers are generated from [0, 1), so that right endpoint will never
99+
// be sampled.
100+
// We have one special case to consider. When u==0.5, we have tan(PI/2),
101+
// which is undefined. In the limit, however, tan(PI/2) is infinity.
102+
// We could map this to the constant for infinity. However, this would introduce
103+
// a very slight bias in favor of positive results since our interval considers
104+
// from tan(0) through tan(PI-epsilon), which doesn't include tan(-PI/2), though it
105+
// comes close since tan(PI/2+epsilon)==tan(-PI/2+epsilon). In the limit, tan(-PI/2)
106+
// is -infinity. So mapping tan(PI/2) to infinity would result in one extra value that
107+
// leads to a positive result relative to the number of values that lead to negative results.
108+
// We handle this in the following way. First, when u==0.5, we generate a
109+
// random boolean to control whether u==0.5 means PI/2 or -PI/2. Second, rather than
110+
// map to the constants for positive and negative infinity from the Double class, we
111+
// pass these along to the tan method and let it do its thing numerically, which is
112+
// a value around 1.6ish * 10 to the power 16 (and negative of that in the case of -PI/2).
113+
double u = ThreadLocalRandom.current().nextDouble();
114+
if (u == 0.5 && ThreadLocalRandom.current().nextBoolean()) {
115+
u = -0.5;
116+
}
117+
return median + scale * StrictMath.tan(StrictMath.PI * u);
97118
}
98119

99120
/**
@@ -104,14 +125,35 @@ public static double nextCauchy(double median, double scale) {
104125
*/
105126
public static double nextCauchy(double scale) {
106127
// 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));
128+
// Mathematically, it should be: median + scale * tan(PI * (u - 0.5)),
129+
// where u is uniformly random from the interval [0, 1].
130+
// However, since tan goes through one complete cycle every PI,
131+
// we can replace it with: median + scale * tan(PI * u), going from
132+
// 0 to PI, rather than from -PI/2 to PI/2. This is equivalent
133+
// as far as generating a random Cauchy variate is concerned, but saves
134+
// one arithmetic operation.
135+
// At first glance, it may appear as if we might be
136+
// doubly sampling u == 0 since tan(0)==tan(PI), however, our uniform
137+
// random numbers are generated from [0, 1), so that right endpoint will never
138+
// be sampled.
139+
// We have one special case to consider. When u==0.5, we have tan(PI/2),
140+
// which is undefined. In the limit, however, tan(PI/2) is infinity.
141+
// We could map this to the constant for infinity. However, this would introduce
142+
// a very slight bias in favor of positive results since our interval considers
143+
// from tan(0) through tan(PI-epsilon), which doesn't include tan(-PI/2), though it
144+
// comes close since tan(PI/2+epsilon)==tan(-PI/2+epsilon). In the limit, tan(-PI/2)
145+
// is -infinity. So mapping tan(PI/2) to infinity would result in one extra value that
146+
// leads to a positive result relative to the number of values that lead to negative results.
147+
// We handle this in the following way. First, when u==0.5, we generate a
148+
// random boolean to control whether u==0.5 means PI/2 or -PI/2. Second, rather than
149+
// map to the constants for positive and negative infinity from the Double class, we
150+
// pass these along to the tan method and let it do its thing numerically, which is
151+
// a value around 1.6ish * 10 to the power 16 (and negative of that in the case of -PI/2).
152+
double u = ThreadLocalRandom.current().nextDouble();
153+
if (u == 0.5 && ThreadLocalRandom.current().nextBoolean()) {
154+
u = -0.5;
155+
}
156+
return scale * StrictMath.tan(StrictMath.PI * u);
115157
}
116158

117159
/**
@@ -123,14 +165,35 @@ public static double nextCauchy(double scale) {
123165
*/
124166
public static double nextCauchy(double median, double scale, Random r) {
125167
// 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));
168+
// Mathematically, it should be: median + scale * tan(PI * (u - 0.5)),
169+
// where u is uniformly random from the interval [0, 1].
170+
// However, since tan goes through one complete cycle every PI,
171+
// we can replace it with: median + scale * tan(PI * u), going from
172+
// 0 to PI, rather than from -PI/2 to PI/2. This is equivalent
173+
// as far as generating a random Cauchy variate is concerned, but saves
174+
// one arithmetic operation.
175+
// At first glance, it may appear as if we might be
176+
// doubly sampling u == 0 since tan(0)==tan(PI), however, our uniform
177+
// random numbers are generated from [0, 1), so that right endpoint will never
178+
// be sampled.
179+
// We have one special case to consider. When u==0.5, we have tan(PI/2),
180+
// which is undefined. In the limit, however, tan(PI/2) is infinity.
181+
// We could map this to the constant for infinity. However, this would introduce
182+
// a very slight bias in favor of positive results since our interval considers
183+
// from tan(0) through tan(PI-epsilon), which doesn't include tan(-PI/2), though it
184+
// comes close since tan(PI/2+epsilon)==tan(-PI/2+epsilon). In the limit, tan(-PI/2)
185+
// is -infinity. So mapping tan(PI/2) to infinity would result in one extra value that
186+
// leads to a positive result relative to the number of values that lead to negative results.
187+
// We handle this in the following way. First, when u==0.5, we generate a
188+
// random boolean to control whether u==0.5 means PI/2 or -PI/2. Second, rather than
189+
// map to the constants for positive and negative infinity from the Double class, we
190+
// pass these along to the tan method and let it do its thing numerically, which is
191+
// a value around 1.6ish * 10 to the power 16 (and negative of that in the case of -PI/2).
192+
double u = r.nextDouble();
193+
if (u == 0.5 && r.nextBoolean()) {
194+
u = -0.5;
195+
}
196+
return median + scale * StrictMath.tan(StrictMath.PI * u);
134197
}
135198

136199
/**
@@ -142,14 +205,35 @@ public static double nextCauchy(double median, double scale, Random r) {
142205
*/
143206
public static double nextCauchy(double scale, Random r) {
144207
// 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));
208+
// Mathematically, it should be: median + scale * tan(PI * (u - 0.5)),
209+
// where u is uniformly random from the interval [0, 1].
210+
// However, since tan goes through one complete cycle every PI,
211+
// we can replace it with: median + scale * tan(PI * u), going from
212+
// 0 to PI, rather than from -PI/2 to PI/2. This is equivalent
213+
// as far as generating a random Cauchy variate is concerned, but saves
214+
// one arithmetic operation.
215+
// At first glance, it may appear as if we might be
216+
// doubly sampling u == 0 since tan(0)==tan(PI), however, our uniform
217+
// random numbers are generated from [0, 1), so that right endpoint will never
218+
// be sampled.
219+
// We have one special case to consider. When u==0.5, we have tan(PI/2),
220+
// which is undefined. In the limit, however, tan(PI/2) is infinity.
221+
// We could map this to the constant for infinity. However, this would introduce
222+
// a very slight bias in favor of positive results since our interval considers
223+
// from tan(0) through tan(PI-epsilon), which doesn't include tan(-PI/2), though it
224+
// comes close since tan(PI/2+epsilon)==tan(-PI/2+epsilon). In the limit, tan(-PI/2)
225+
// is -infinity. So mapping tan(PI/2) to infinity would result in one extra value that
226+
// leads to a positive result relative to the number of values that lead to negative results.
227+
// We handle this in the following way. First, when u==0.5, we generate a
228+
// random boolean to control whether u==0.5 means PI/2 or -PI/2. Second, rather than
229+
// map to the constants for positive and negative infinity from the Double class, we
230+
// pass these along to the tan method and let it do its thing numerically, which is
231+
// a value around 1.6ish * 10 to the power 16 (and negative of that in the case of -PI/2).
232+
double u = r.nextDouble();
233+
if (u == 0.5 && r.nextBoolean()) {
234+
u = -0.5;
235+
}
236+
return scale * StrictMath.tan(StrictMath.PI * u);
153237
}
154238

155239
/**
@@ -161,14 +245,35 @@ public static double nextCauchy(double scale, Random r) {
161245
*/
162246
public static double nextCauchy(double median, double scale, SplittableRandom r) {
163247
// 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));
248+
// Mathematically, it should be: median + scale * tan(PI * (u - 0.5)),
249+
// where u is uniformly random from the interval [0, 1].
250+
// However, since tan goes through one complete cycle every PI,
251+
// we can replace it with: median + scale * tan(PI * u), going from
252+
// 0 to PI, rather than from -PI/2 to PI/2. This is equivalent
253+
// as far as generating a random Cauchy variate is concerned, but saves
254+
// one arithmetic operation.
255+
// At first glance, it may appear as if we might be
256+
// doubly sampling u == 0 since tan(0)==tan(PI), however, our uniform
257+
// random numbers are generated from [0, 1), so that right endpoint will never
258+
// be sampled.
259+
// We have one special case to consider. When u==0.5, we have tan(PI/2),
260+
// which is undefined. In the limit, however, tan(PI/2) is infinity.
261+
// We could map this to the constant for infinity. However, this would introduce
262+
// a very slight bias in favor of positive results since our interval considers
263+
// from tan(0) through tan(PI-epsilon), which doesn't include tan(-PI/2), though it
264+
// comes close since tan(PI/2+epsilon)==tan(-PI/2+epsilon). In the limit, tan(-PI/2)
265+
// is -infinity. So mapping tan(PI/2) to infinity would result in one extra value that
266+
// leads to a positive result relative to the number of values that lead to negative results.
267+
// We handle this in the following way. First, when u==0.5, we generate a
268+
// random boolean to control whether u==0.5 means PI/2 or -PI/2. Second, rather than
269+
// map to the constants for positive and negative infinity from the Double class, we
270+
// pass these along to the tan method and let it do its thing numerically, which is
271+
// a value around 1.6ish * 10 to the power 16 (and negative of that in the case of -PI/2).
272+
double u = r.nextDouble();
273+
if (u == 0.5 && r.nextBoolean()) {
274+
u = -0.5;
275+
}
276+
return median + scale * StrictMath.tan(StrictMath.PI * u);
172277
}
173278

174279
/**
@@ -180,14 +285,35 @@ public static double nextCauchy(double median, double scale, SplittableRandom r)
180285
*/
181286
public static double nextCauchy(double scale, SplittableRandom r) {
182287
// 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));
288+
// Mathematically, it should be: median + scale * tan(PI * (u - 0.5)),
289+
// where u is uniformly random from the interval [0, 1].
290+
// However, since tan goes through one complete cycle every PI,
291+
// we can replace it with: median + scale * tan(PI * u), going from
292+
// 0 to PI, rather than from -PI/2 to PI/2. This is equivalent
293+
// as far as generating a random Cauchy variate is concerned, but saves
294+
// one arithmetic operation.
295+
// At first glance, it may appear as if we might be
296+
// doubly sampling u == 0 since tan(0)==tan(PI), however, our uniform
297+
// random numbers are generated from [0, 1), so that right endpoint will never
298+
// be sampled.
299+
// We have one special case to consider. When u==0.5, we have tan(PI/2),
300+
// which is undefined. In the limit, however, tan(PI/2) is infinity.
301+
// We could map this to the constant for infinity. However, this would introduce
302+
// a very slight bias in favor of positive results since our interval considers
303+
// from tan(0) through tan(PI-epsilon), which doesn't include tan(-PI/2), though it
304+
// comes close since tan(PI/2+epsilon)==tan(-PI/2+epsilon). In the limit, tan(-PI/2)
305+
// is -infinity. So mapping tan(PI/2) to infinity would result in one extra value that
306+
// leads to a positive result relative to the number of values that lead to negative results.
307+
// We handle this in the following way. First, when u==0.5, we generate a
308+
// random boolean to control whether u==0.5 means PI/2 or -PI/2. Second, rather than
309+
// map to the constants for positive and negative infinity from the Double class, we
310+
// pass these along to the tan method and let it do its thing numerically, which is
311+
// a value around 1.6ish * 10 to the power 16 (and negative of that in the case of -PI/2).
312+
double u = r.nextDouble();
313+
if (u == 0.5 && r.nextBoolean()) {
314+
u = -0.5;
315+
}
316+
return scale * StrictMath.tan(StrictMath.PI * u);
191317
}
192318

193319
}

0 commit comments

Comments
 (0)