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?
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/τ)
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.
Herbie gives some useful suggestions here. https://herbie.uwplse.org/demo/e7768dd3187b2711e9ab6825ebe15f0057554539.ea9754c765d6f2e3443236047de835c7bf56f0de/graph.html (although the last branch is incorrect, it should just be return t_0.
This is unrelated, but I get a very different answer from yours:
https://herbie.uwplse.org/demo/1b508faa0f15ca8e4e244eadb39ba19c403cdcf4.ea9754c765d6f2e3443236047de835c7bf56f0de/graph.html
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.
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.
https://numpy.org/devdocs/reference/random/generated/numpy.random.rand.html
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 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)
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 [0,1).
Which numbers are left out?
rand() generates numbers that are i*2^-53 where i is a uniform integer in [0,2^53)
Eh, it’s [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 0.0
Why not use the Copulas.jl package ?
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:
bigC = ClaytonCopula(2,big(300))
bigX = SklarDist(bigC,(X₁,X₂))
Random.rand!(bigX,matkk')
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.