Code to generate Gaussian (normally distributed) random numbers in Ruby

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, yend

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  endend

(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.

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

+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  endend

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