This is a continuation of a previous post where I wanted to create a random uniform distribution between 0 (which is excluded) and 10 (included). I am making a second post but keeping the range general in the hope it will help someone one day.
So for example:
I would like an array of length 8
In this array I would like to randomly generate 8 values
These values fall uniformly within the ranges of 2 and 7.
I guess this could also be phrased as,
How to generate 8 random uniform values between 2 and 7, and place them within an array?
Another way of thinking of uniformly distributed random values between 2 and 7 is uniformly distributed random values between 0 and 5 and then adding 2 to those. The rand function generates values in [0, 1], so if you multiply that by 5, you get values in [0, 5].
The range of value is from 2 to 7 therefore the span is 7 - 2 = 5.
The offset is 2 because the smallest number is 2.
Create a random number generator with span 5 and offset 2
julia> function myrng()
return 2.0 + 5.0 * rand()
end
myrng (generic function with 1 method)
finally create your array
julia> a = [ myrng() for n in 1:8 ]
8-element Array{Float64,1}:
4.636049027494998
5.061366873363691
6.420185007452647
6.955145662491416
2.708159366136055
6.223899415297151
6.682556947797874
4.927220645178398
If this is a common problem, then perhaps the best way to help others is to make a small PR to Distributions with a variant of Uniform that does this. Or even generalizing the idea to other (half)bounded distributions.
I recall from your previous question that there was a big discussion on your requirement that you want a random number in the range (F1,F2], in other words that F1 is excluded and F2 included. The “problem” is that rand() defaults to producing numbers in the range [0,1), i.e., lower bound included and upper bound excluded. In that thread, someone suggested that 1-rand() produces numbers in the range (0,1]. Thus, I think you can do:
function my_rand(F1,F2,dims...)
return F1 .+ (F2-F1)*(1 .- rand(dims...))
end
If you want a function that returns a scalar when dims is skipped, you can add the method:
function my_rand(F1,F2)
return F1 + (F2-F1)*(1-rand())
end
A couple of examples (assuming the functions have been defined):
Actually, as the last bit of a value produced by rand is 0, you need to call rand() / prevfloat(1.0, 2) in order to have 1.0 included in the possible outcomes.
Stepping back a bit from the immediate problem, the best approach may not be designing a function that can ensure random draws \in (a, b], but an algorithm that can deal with =a robustly.
It is common to have modeling problems where eg some x \in D holds for an open subset of \mathbb{R}^n. Then for continuous functions f_i, it of course follows that f_1(x) \in f_1(D), f_2(f_1(x)) \in f_2(f_1(D)), etc.
But for floating point, numerical error will inevitably creep in, especially for complex transformations. At that point the algorithm should be able to cope with x \in D but f_n(\dots) being “invalid”, eg
by returning a -Inf log likelihood/posterior for likelihood-based inference,
by not recording a simulated sample and trying again for MC simulation,
“regularizing” the value so that it is valid (eg for x < a, replace with a+\epsilon, this is sometimes not innocuous though).