Normally distributed random values

gistfile1.rb

#!/usr/bin/env ruby

require 'securerandom'
require 'pp'

# NormalRand is an attempt to port the Java implementation from Wikipedia
# <https://en.wikipedia.org/wiki/Marsaglia_polar_method>
#
# private static double spare;
# private static boolean isSpareReady = false;
#
# public static synchronized double getGaussian(double mean, double stdDev) {
#     if (isSpareReady) {
#         isSpareReady = false;
#         return spare * stdDev + mean;
#     } else {
#         double u, v, s;
#         do {
#             u = Math.random() * 2 - 1;
#             v = Math.random() * 2 - 1;
#             s = u * u + v * v;
#         } while (s >= 1 || s == 0);
#         double mul = Math.sqrt(-2.0 * Math.log(s) / s);
#         spare = v * mul;
#         isSpareReady = true;
#         return mean + stdDev * u * mul;
#     }
# }

class NormalRand
  include Enumerable

  RAND_GEN = SecureRandom.method(:random_number)

  def initialize(options = {})
    @mean = options.fetch(:mean)
    @stddev = options.fetch(:stddev)
    @rand_gen = options.fetch(:rand_gen, RAND_GEN)
  end

  def next
    enumerator.next
  end

  def each
    if block_given?
      enumerator.each { |value| yield value }
    else
      enumerator.each
    end
  end

  # Create a hash of the frequency of n random values classified by the
  # corresponding rounded integer.
  # Useful for testing purposes, to check if the random values actually
  # satisfy a normal distribution.
  def rounded_frequencies(n_values)
    pairs = self.take(n_values)
      .group_by { |float| float.round }
      .map { |rounded_value, close_values| [rounded_value, close_values.size] }
      .sort

    Hash[pairs]
  end

  private

  def enumerator
    @enumerator ||= Enumerator.new do |yielder|
      loop do
        if @next_gaussian_rand
          yielder.yield @next_gaussian_rand
          @next_gaussian_rand = nil
        else
          x, y, s = 0, 0, 0
          loop do
            x = rand_between_minus_1_and_1
            y = rand_between_minus_1_and_1
            s = x**2 + y**2
            break if 0 < s && s < 1
          end

          coeff = Math.sqrt(-2.0 * Math.log(s) / s)

          @next_gaussian_rand = @mean + @stddev * y * coeff

          yielder.yield @mean + @stddev * x * coeff
        end
      end
    end
  end

  def rand_between_minus_1_and_1
    @rand_gen.() * 2 - 1
  end
end

# This is to verify that the gaussian random generator is working reasonably
class NormalDistribution
  def initialize(options={})
    @mean = options.fetch(:mean)
    @stddev = options.fetch(:stddev)
    @sample_size = options.fetch(:sample_size, 1)
  end

  def f(x)
    coeff = 1 / (@stddev * Math.sqrt(2 * Math::PI))
    exp = - (x - @mean)**2 / (2 * @stddev**2)
    coeff * (Math::E**exp)
  end

  def for_range(range, step=1)
    {}.tap do |result|
      range.step(step).each { |x| result[x] = f(x) * @sample_size }
    end
  end
end

# Made up values
mean = 23.098
stddev = 4.420

random_gaussian = NormalRand.new(mean: mean, stddev: stddev)
random_gaussian_frequencies = random_gaussian.rounded_frequencies(30_000)

pp random_gaussian_frequencies
pp random_gaussian_frequencies.values.inject(:+)

################################################################################

normal_distribution = NormalDistribution.new(mean: mean, stddev: stddev, sample_size: 30_000)

pp normal_distribution.for_range(0..50)
pp normal_distribution.for_range(0..50).values.inject(:+)