How to create a random Uniform Distribution between (but excluding) 0 and 10?

New to Julia
1

Hello,

How can I create a random uniform distribution between 0 (which is excluded) and 10 (included).
That is, the range is between 0 and 10, but 0 itself cannot be included.

using Distributions
rand(Uniform(0, 10)

Using rand(Uniform(>0, 10) does not work.

2 Likes
2

Technically, the endpoints have a very small probability (theoretically 0, in practice \approx0 because of floating point). Depending on your application, you can reject and redraw.

Minimal example:

using Random

"`Uniform(0,b)`, with `0` excluded for sure, and we really mean it."
struct PositiveUniform{T}
    b::T
end

function Base.rand(rng::Random.AbstractRNG, pu::PositiveUniform)
    while true
        r = rand(rng)
        r > 0 && return r * pu.b
    end
end

rand(PositiveUniform(10))

I am curious why you expected it to work. Adding a closing ),

julia> rand(Uniform(>0, 10))
ERROR: syntax: ">" is not a unary operator

informs you that it fails at parsing >0. You are dealing with a programming language, and can’t just make up arbitrary expressions and hope it will do something sensible.

11 Likes
3

rand() samples a random uniform in [0,1), so 1-rand() is a random uniform in (0,1], and 10*(1-rand()) isa rnadom uniform in (0,10].

This is because the PRNG used (and most PRNGS) generate random unsigned integers, and then mask them to produce a number between 1.0 and prevfloat(2.0).

julia> bitstring(1.0)
"0011111111110000000000000000000000000000000000000000000000000000"

julia> bitstring(prevfloat(2.0))
"0011111111111111111111111111111111111111111111111111111111111111"

The last 52 bits are the fraction component. By generating a random unsigned integer, and setting the first 12 bits to 001111111111, you get a random Float64 between 1.0 and prevfloat(2.0).
Then, most random number generators subtract 1, to get the [0,1) range.
In principle, you could skip a step and use:

using Random
2.0 - rand(Random.GLOBAL_RNG, Random.SamplerTrivial(Random.CloseOpen12(Float64)))

instead. Or:

@inline mask(x::UInt64) = reinterpret(Float64,(x & 0x000fffffffffffff) | 0x3ff0000000000000)
@inline mask(x::UInt32) = reinterpret(Float32,(x & 0x007fffff) | 0x3f800000)

2.0 - mask(rand(UInt64))

You can look at the Julia code for rand functions here.

12 Likes
5

Thank you for your replies. I confess I don’t really understand them - they are at a higher level of coding/Julia than I am currently at, but I understand it is not as straightforward as I had hoped. Nonetheless I will go over and try to understand them.

Could someone please explain what the difference between [0, 1) and (0, 1] is?

Do you recommend I use the Random instead of Distributions package, because the rand() in Random already samples a random uniform number (or did I misunderstand)?

6

Maybe you tell us what you are up to. You didn’t even tell if you look for integer or floating point numbers, maybe that’s the reason why you are now confused (all the floating point intrinsics)?

2 Likes
7

It’s not code, but standard math notation for open and closed intervals.

2 Likes
8

I am writing code from a research paper, but they just state “…animal and plant biomass densities were initialized with random values uniformly distributed between 0 (exclusive) and 10 (inclusive).” I assume the numbers can contain decimals.

9

You might want to give a link to the article.

10

https://www.nature.com/articles/ncomms12718.pdf?origin=ppub

A Google search :wink:

11

Correct. I am trying to translate the code into Julia. The code they used is in C++, but that does not help me because I don’t know that language at all.

12

hahaha! Should have thought about it.

13

It might be easier to understand if you go from the article as the prime inspiration and look at the C code for confirmation? (just an academic suggestion).

You might want to try defining your dynamical system with DynamicalSystems.jl

edit - For your initial question @Tamas_Papp provided an elegant solution.

1 Like
14

This seems to answer your question though? As @Tamas_Papp says (0, 10] is just standard mathematical notation for a set excluding 0 and including 10.

So:

julia> using Plots

julia> a = 10 .* (1 .- rand(100_000))

julia> minimum(a)
1.1877405667881646e-5

julia> histogram(a)

randuni

Looks good to me!?

1 Like
15

If you want a random umber generator with (0,1]
but all you got is a rng with [0,1) then why not create your own

my_rng = 1 - rng

so my_rng gives the distribution (0,1]

to get (0,10] you just multiply it by ten

RNG = 10 * my_rng

16

rand() samples uniform fixed-point numbers in [0, 1) that are then converted to floating point. It achieves this by sampling a uniform floating point number in [1,2) and then subtracting 1.

This means that the output distribution of rand() is “ragged”: Instead of the naively expected (sample real number in [0,1), round to the next smaller or equal non-subnormal float), which would make exact 0 vanishingly rare, zero results are very likely.

E.g. rand(Float32) has a probability of eps(Float32) ~ 1e-7 of returning exactly zero, compared to the correct probability of 1e-38 (using normal floats only) or 1e-45 (permitting subnormals, which would be a bad idea). A correct distribution close to 0 is probably prohibitively expensive, but 2e-10 ~ 1/(1<<32) could be done branch-free by spending 32 bits of entropy instead of 23 bits.

5 Likes
17

Let me make a little contribution. Can’t you just run

rand(0.0+eps(),10.0-eps())

?

Simple and works.

18

Nope, it doesn’t work. rand has no method for this — it would need something like Distributions.Uniform.

Also, even though sometimes one is forced to make regularizations like this as a quick & dirty fix, correct solutions are actually easy in this case.

1 Like
19

I imagine this would also be computationally expensive to do with only 32 bits, or?
Converting a random string of 32 or 64 bits into 23 or 52 random bits of a uniform [1,2) distribution only requires two vectorizable and fast CPU instructions.
Just set the sign and exponential bits correctly via a bitwise & and then a bitwise |.

How would you actually set the exponential bits for a non-ragged unfirom (0,1)?
There are an equal number of floating point numbers in the range [0.25, 0.5) as there is [0.5, 1) or [0.125, 0.25).
So it sounds like you’d need quote some logic in actually processing your remaining bits of entropy to create that correct distribution in the upper bits.

While this uses twice as many bits as sizeof(T), it should still be relatively fast:

using Random
function nonraggedrand(rng = Random.GLOBAL_RNG)
    zero_sgnexp = 0x000fffffffffffff
    zero_frac = 0xfff0000000000000
    random_exponent = reinterpret(UInt64,rand(rng)) & zero_frac
    random_fraction = rand(rng, UInt64) & zero_sgnexp
    reinterpret(Float64, random_exponent | random_fraction)
end

Looks correct from a superficial glance:

julia> using RNGTest

julia> RNGTest.smallcrushJulia(nonraggedrand)
10-element Array{Any,1}:
 0.5336216850647422                                                                                   
 0.9230989723619558                                                                                   
 0.49280958521636054                                                                                  
 0.6834132908358886                                                                                   
 0.8940787394692569                                                                                   
  (0.5333583536888169, 0.5896538779593443)                                                            
 0.37459169125750136                                                                                  
 0.14003119092065008                                                                                  
 0.7706765198832467                                                                                   
  (0.5617563416418929, 0.9608245561823281, 0.37088087074827114, 0.839032775363961, 0.6894691545291662)

The idea is: the ragged randoms have the correct distribution of exponential bits. So, just combine those exponential bits with random fraction bits.

This uses 128 bits for a 64 bit random number, or 64 bits for a 32 bit random number. But I think that’s rather reasonable.
Many random number generators intentionally discard extra random bits anyway, to help with things like the “birthday problem”, or just for the sake of having a larger state for a longer period.

Maybe it’s a better idea to use a PCG RXS-M-XS generator, and discard state via combining floating points like this, than it is to use a PCG XSH-RS generator, for example?
Those generators can be fast and are vectorizable. If you’re only sampling serially, vectorization means the cost will be less than 2x for taking this approach.

20

Just to try to prevent this thread being overwhelming for a newcomer, let me just highlight the straightforward answer:

10*(1-rand())

Analysis:

  • rand() is uniform from [0, 1)
  • 1-rand() is uniform from (0, 1]
  • 10*(1-rand()) is uniform from (0, 10]

This is a totally fine, simple answer. Is it perfectly uniformly distributed? Meh, maybe, maybe not but the chance that it matters for this use case is negligible.

22 Likes
21

If you want to exclude both 0 and 10, you can use prevfloat(10.0)*(1-rand()) instead.
Same argument you gave applies, the answer being uniform from (0, prevfloat(10.0)], i.e. (0,10).

If you don’t, the probability of getting a 10 on any given draw is 1/2^52 = 2.220446049250313e-16.

3 Likes