In some simulations I am encountering a problem due to machine precision. I have to conduct an exponentiation operation on a number that lies in (0, 1). That number comes from a random number generator, while the exponent is determined by an optimization algorithm.
My question is how I can determine the maximum exponent
n such that
rand()^n does not result in
Or equivalently, how I can determine the maximum exponent
n such that
rand()^(-n) does not result in
These aren’t actually equivalent questions due to the way subnormals work. The key things to realize are that
rand always generates a number in
(floatmin(),1), and therefore the answer is the power which will cause
floatmin() to under/overflow. This number is
log(floatmax())/log(floatmin()) respectively. (
Thank you. It makes sense that the numbers are so low.
I was thinking of using
clamp to limit the range of
n, but ~0.95 is very low for my application. I am expecting
n be up to 100 to 150. I understand that the operation depends on a random number so I can never know for sure what number will be drawn, but do you have any suggestions regarding that?
What you actually want to do is generate a random number with the desired distribution.
Specifically, what I am trying to do is generate a correlated exponential distribution using the Clayton copula. The input is simply a matrix of random numbers in (0, 1):
randmat = rand(100, 2)
matkk = similar(randmat)
function claytonsample!(matkk, τ; randmat=randmat)
matkk .= randmat
τ == 0 && return matkk
n = size(matkk, 1)
for i in 1:n
v = matkk[i, 2]
u = matkk[i, 1]
matkk[i, 2] = (1 - u^(-τ) + u^(-τ)*v^(-(τ/(1 + τ))))^(-1/τ)
invsurvival(x) = -log(x)
function simulkk!(matkk, τ; randmat=randmat)
claytonsample!(matkk, τ; randmat=randmat)
matkk .= invsurvival.(matkk)
simulkk!(matkk, 300; randmat=randmat)
The problem is in the terms
u^(-τ) above. I was thinking of clamping
τ. I’d rather not mess with the “seed” random numbers in
Herbie gives some useful suggestions here. Result for (1 - x^(-y) + x^(-y)*z^(-(y/(1 + y))))^(-1/y) (although the last branch is incorrect, it should just be return
I might have typed something in wrong, never mind, herbie just got drunk when evaluating yours.
Oh man, that is some serious drinking problem. Hopefully both are equivalent.
[0, 1) or
[floatmin(), 1) (so basically,
(floatmin(), 1)? I can’t find it in the docs, but numpy’s
np.random.rand() claims to be
[0,1) and I imagine Julia would follow the same convention.
Posting in the same topic because I have a fundamentally similar question–I want to know whether I can safely do
1 / rand(Float16).
Numerically, no matter what you’re doing, division by anything even close to
floatmin will be unstable. It will however be “safe”, in that it will not error out. Even division by
0.0 is not an error, it just returns
In my case, it arises in a unit test that uses randomly generated input data, so while we would expect higher numbers in “real” input, all the test requires is that
1/rand() be finite.
Don’t you want your unit tests to use realistic numbers? Or just add
1.0 to the random numbers to be safe (or maybe put in some handling of
Inf in your program)
Sure, depends on the context.
A secondary motivation for my question is so that I can write a sentence like "in this experiment,
x was drawn randomly from [0,1)" in a scientific paper describing a computational experiment (one that doesn’t involve computing
1/x). I would like my sentence to be accurate.
rand() is uniform over
[0,1) but the posible return values are not all floating point numbers in
Which numbers are left out?
rand() generates numbers that are
i is a uniform integer in
[0, 1) with
p<0.05. You’re more likely to have a cosmic bit flip place you outside the interval than for it to literally pick