Julia random number generation

Julia is a new programming language for scientific computing. From the Julia site:

Julia is a high-level, high-performance dynamic programming language for technical computing, with syntax that is familiar to users of other technical computing environments. It provides a sophisticated compiler, distributed parallel execution, numerical accuracy, and an extensive mathematical function library. …

I just started playing around with it. I didn’t see functions for non-uniform random number generation so I wrote some as a way to get started.

[Update: there are non-uniform random number generators in Julia, but they have not been added to the documentation yet. See details in this comment.]

Here’s a random number generator for normal (Gaussian) random values:

## return a random sample from a normal (Gaussian) distribution
function rand_normal(mean, stdev)
    if stdev <= 0.0
        error("standard deviation must be positive")
    end
    u1 = rand()
    u2 = rand()
    r = sqrt( -2.0*log(u1) )
    theta = 2.0*pi*u2
    mean + stdev*r*sin(theta)
end

From this you can see Julia is a low-ceremony language: Python-like syntax, you can call common mathematical functions without having to do anything special, etc. You can have explicit return statements, but the preferred style seems to be to let the last line of the function be the implicit return statement.

My most common mistake so far has been forgetting to close code blocks with end; Julia’s syntax is similar enough to Python that I suppose I think indentation should be sufficient.

I’ve written random number generators for the following probability distributions:

  • Beta
  • Cauchy
  • Chi square
  • Exponential
  • Inverse gamma
  • Laplace (double exponential)
  • Normal
  • Student t
  • Uniform
  • Weibull

You can find the code here: Non-uniform random number generation in Julia.

Tagged with: ,
Posted in Software development
12 comments on “Julia random number generation
  1. Mark Henderson says:

    Julia also caught my eye. At the end of your post I was hoping to see more of your opinion of the language based upon first impression.

  2. John says:

    First impressions? The syntax is simple and predictable, easy to pick up. I think the designers of the language have good taste. The performance benchmarks are impressive.

    I haven’t tried it for anything more than simple code, so I can’t comment on how advanced features work.

    It’s new, so it has some maturing to do. The documentation is well-written, but it needs to fill in gaps. There’s no Windows version. The Linux build failed for me, though I was able to download and use a binary distribution.

  3. Wow. I’m a huge fan of this blog — we link to several of your articles in our documentation:

    http://julialang.org/manual/integers-and-floating-point-numbers/#Background+and+References
    http://julialang.org/manual/mathematical-operations/#Mathematical+Functions

    There are actually a bunch of non-uniform random-number generation functions already, but the documentation isn’t complete. We currently have the following:

    rand — uniform
    randn — Gaussian
    randchi2 — Chi-squared
    randexp — exponential
    randg — Gamma

    Most of these are implemented such that filling vectors is efficient and where applicable, using fast algorithms like Ziggurat. There’s also a bunch of random integer and permutation functions. You can type “rand” in the repl and hit tab to see all the possible completions. We really should flesh out the standard library reference, but that’s a fair amount of work that’s less fun than coding.

    Thanks for trying it out and for the extra rand functions. I’ll try to integrate these.

  4. Hi John,

    That’s pretty cool. If ok, we’ll import some of those RNGs. We do have normally distributed random numbers – randn, and also randg. Perhaps they are not documented.

  5. I think one thing that’s really nice about Julia is that you were able to write all these functions without worrying about type (like you can in Python) but if you need it it also has a really nice type system with clever function dispatching.

    For example if you found you had an optimized rand_normal that only applied to integral means you could just also define

    function rand_normal(mean::Integer, stdev)

    then if you run rand_normal with an integral mean the optimized version would be called. Otherwise the generic version would be.

  6. John says:

    Stefan and Viral: Thanks for creating Julia. I’m honored to be referenced in the documentation.

    That’s great that you’re including more sophisticated RNG methods than those mentioned here. I chose algorithms for simplicity, not efficiency. Methods like the zigguarat will be more efficient, as will generating multiple values at once for vectors.

  7. Jorge Stolfi says:

    I see that Julia does not require variables to be declared, so typos on the left side of “=” are not detected by the compiler. I have been burned by this “feature” in (g)awk enough times. It saves keystrokes but is not good for large programs.

    Also, the random generators should be objects with rand() etc. as methods. Certain applications (e.g. in graphics) require multiple independent & repeatable pseudorandom streams.

  8. MattyG says:

    Neat. I saw some really wicked benchmarks float around about this laungage, and I am keen to try it in the next couple days.

  9. Viral B. Shah says:

    John, one of the fun things I was able to also do was to implement Ziggurat in julia:

    https://github.com/JuliaLang/julia/blob/master/examples/ziggurat.j

    The last I had tried a while back, it was about twice as slow as native C. That is why we use the native C version in the standard library until Julia gets faster.

  10. SteveBrooklineMA says:

    I was going to ask “why?”, then I read this

    http://julialang.org/blog/2012/02/why-we-created-julia/

    It sounds like they made it for me.

  11. @Jorge:

    Not detecting misspellings on the right-hand-side of = is certainly a tradeoff of not requiring variables to be declared. However, I’ve done a lot of programming in a wide variety of languages that have this sort of behavior and I’ve really never found it to be a big problem. Maybe I’m just a careful typist? When such bugs do occur, it can take a while to track down, but on the whole I think I’ve gained more time by not declaring variables than I’ve lost to such bugs. It may also be possible to heuristically detect and warn of cases where a name is assigned and never used. Matlab’s IDE does this.

    You’re absolutely right about the random number generators needing to be put into objects so that you can have multiple independent streams of random numbers. That’s an excellent “up-for-grabs” project for someone who wants to get their hands dirty with the language. We’ll retain versions that use a single implicit global state, however, for ease of use. Forcing people to always pass a parameter they don’t care about is far too much “ceremony” for us.

  12. Kevin says:

    On declaring variables before use, Stefan writes:

    “on the whole I think I’ve gained more time by not declaring variables than I’ve lost to such bugs”

    Does it really cost that much time to just type “var” or “let” before the first use of a variable?

    “It may also be possible to heuristically detect and warn of cases where a name is assigned and never used. ”

    The problem is the other case: a name is used but never assigned.

3 Pings/Trackbacks for "Julia random number generation"
  1. [...] Julia random number generation uh! il nuovo linguaggio intriga anche John D. ::: The Endeavour [...]

  2. [...] John D Cook, introduced this Julia in this blog post and generously shared the random number generator codes in Julia for probability distribution, I [...]

  3. [...] John D Cook writes some simple non-uniform random number generators in Julia. [...]