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 0.0?

Or equivalently, how I can determine the maximum exponent n such that rand()^(-n) does not result in Inf?

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(nextfloat(0.0))/log(floatmax()) or log(floatmax())/log(floatmin()) respectively. (0.9534450651769087 or -1.0019569471624266)

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?

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/τ)
end
return matkk
end
invsurvival(x) = -log(x)
function simulkk!(matkk, τ; randmat=randmat)
claytonsample!(matkk, τ; randmat=randmat)
matkk .= invsurvival.(matkk)
end
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 randmat.

Wait, is rand() in [0, 1) or [floatmin(), 1) (so basically, (0,1)) or (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.

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 Inf.

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)

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.

using Copulas, Distributions, Random
C = ClaytonCopula(2,300)
X₁ = Exponential()
X₂ = Exponential()
X = SklarDist(C,(X₁,X₂))
Random.rand!(X,matkk')

The same overflows appends of course, but here returns zeros and not NaN. If you dont want the overflows, however, moving to BigFloats sampling is really easy:

Note that \theta=300 is insanely close to a pure comonotony, especially in the lower tail, which is why your inverse Rosenblat transformation is failling when u is too small. You could also just pass bigfloats random numbers to your rosenblat transformation, it’ll work too.

See on the github page for more details: many other copulas are available.