29

What is some code to generate normally distributed random numbers in ruby?

(Note: I answered my own question, but I'll wait a few days before accepting to see if anyone has a better answer.)

EDIT:

Searching for this, I looked at all pages on SO resulting from the two searches:

+"normal distribution" ruby

and

+gaussian +random ruby

Eponymous
  • 6,143
  • 4
  • 43
  • 43
  • Did your check related question (see right side panel)? – Gumbo Apr 28 '11 at 22:19
  • Yes, I checked and though there are places that have the algorithm, no one has coded it up in Ruby. It is such a common task that it really should be in the standard library. But failing that, I think copy-paste code should be findable on StackOverflow. – Eponymous Apr 28 '11 at 22:26
  • It might be a good idea to mention what you've checked, so that people thinking of answering won't check them, unless they think you missed something. – Andrew Grimm Apr 28 '11 at 23:34

4 Answers4

53

Python's random.gauss() and Boost's normal_distribution both use the Box-Muller transform, so that should be good enough for Ruby too.

def gaussian(mean, stddev, rand)
  theta = 2 * Math::PI * rand.call
  rho = Math.sqrt(-2 * Math.log(1 - rand.call))
  scale = stddev * rho
  x = mean + scale * Math.cos(theta)
  y = mean + scale * Math.sin(theta)
  return x, y
end

The method can be wrapped up in a class that returns the samples one by one.

class RandomGaussian
  def initialize(mean, stddev, rand_helper = lambda { Kernel.rand })
    @rand_helper = rand_helper
    @mean = mean
    @stddev = stddev
    @valid = false
    @next = 0
  end

  def rand
    if @valid then
      @valid = false
      return @next
    else
      @valid = true
      x, y = self.class.gaussian(@mean, @stddev, @rand_helper)
      @next = y
      return x
    end
  end

  private
  def self.gaussian(mean, stddev, rand)
    theta = 2 * Math::PI * rand.call
    rho = Math.sqrt(-2 * Math.log(1 - rand.call))
    scale = stddev * rho
    x = mean + scale * Math.cos(theta)
    y = mean + scale * Math.sin(theta)
    return x, y
  end
end

CC0 (CC0)

To the extent possible under law, antonakos has waived all copyright and related or neighboring rights to the RandomGaussian Ruby class. This work is published from: Denmark.

idmean
  • 14,540
  • 9
  • 54
  • 83
antonakos
  • 8,261
  • 2
  • 31
  • 34
  • 1
    Could you add a permissive license to your code (BSD/CC-0 or some-such) (since it is intended for cut-paste reuse) – Eponymous May 31 '11 at 20:57
  • @antonakos in my openoffice calc I have a function called NORMDIST with parameters (Number, Mean, STDEV, C), how I can give the 4 parameters using your code ? (your code accepts `mean, stddev, rand` which is different than what I use in my excel/openoffice calc) – medBouzid Mar 01 '18 at 16:57
  • I wonder why you do `Math.log(1 - rand.call)` instead of just `Math.log(rand.call)` since `rand.call` returns a number between 0 and 1. I've checked the python and the boost code, and they too do that... but why? – kaikuchn Jul 26 '18 at 07:39
  • The smallest nonzero float in Ruby is `0.0.next_float => 5.0e-324` the biggest float below `1.0` is `1.0.prev_float => 0.9999999999999999`. I don't know if `Kernel.rand` can return `5.0e-324` but if it does, then the subtraction would truncate the tail of the random variable. – kaikuchn Jul 26 '18 at 08:52
21

The original question asked for code, but the author's followup comment implied an interest in using existing libraries. I was interested in the same, and my searches turned up these two ruby gems:

gsl - "Ruby interface to the GNU Scientific Library" (requires you to install GSL). The calling sequence for normally distributed random numbers with mean = 0 and a given standard deviation is

 rng = GSL::Rng.alloc
 rng.gaussian(sd)      # a single random sample
 rng.gaussian(sd, 100) # 100 random samples

rubystats - "a port of the statistics libraries from PHPMath" (pure ruby). The calling sequence for normally distributed random numbers with a given mean and standard deviation is

 gen = Rubystats::NormalDistribution.new(mean, sd)
 gen.rng               # a single random sample
 gen.rng(100)          # 100 random samples
ronen
  • 2,529
  • 26
  • 23
12

Another option, this one using the distribution gem, written by one of the SciRuby fellows.

It is a little simpler to use, I think.

require 'distribution'
normal = Distribution::Normal.rng(1)
norm_distribution = 1_000.times.map {normal.call}
fny
  • 31,255
  • 16
  • 96
  • 127
Ryanmt
  • 3,215
  • 3
  • 22
  • 23
  • Looking at this code, I have no idea what the expected standard deviation is supposed to be. – Paul Brannan Dec 10 '13 at 20:10
  • 2
    Go here: http://rubydoc.info/gems/distribution/0.7.0/Distribution/Normal/Ruby_ You'll see that the rng(1) is specifying the mean, and that you can specify the standard deviation you want by passing the additional parameter `Distribution::Normal.rng(mean, standard_deviation)` which you call to provide a random value in that distribution. – Ryanmt Dec 19 '13 at 20:24
12

+1 on @antonakos's answer. Here's the implementation of Box-Muller that I've been using; it's essentially identical but slightly tighter code:

class RandomGaussian
  def initialize(mean = 0.0, sd = 1.0, rng = lambda { Kernel.rand })
    @mean, @sd, @rng = mean, sd, rng
    @compute_next_pair = false
  end

  def rand
    if (@compute_next_pair = !@compute_next_pair)
      # Compute a pair of random values with normal distribution.
      # See http://en.wikipedia.org/wiki/Box-Muller_transform
      theta = 2 * Math::PI * @rng.call
      scale = @sd * Math.sqrt(-2 * Math.log(1 - @rng.call))
      @g1 = @mean + scale * Math.sin(theta)
      @g0 = @mean + scale * Math.cos(theta)
    else
      @g1
    end
  end
end

Of course, if you really cared about speed, you should implement the Ziggurat Algorithm :).

fearless_fool
  • 33,645
  • 23
  • 135
  • 217