Untitled

# -*- coding: utf-8 -*-
begin
  require 'ansi'
  def go(color)
    print ANSI.send(color)
  end
rescue LoadError => e
  $stderr.puts "to get colored output: gem install ansi"
  def go(color); nil; end
end


# rolls  = an array holding the individual results
# low    = minimum element value in rolls; for a dN die, this should be 1
# high   = maximum element value in rolls; for a dN die, this should be N
# counts = an array indexed by the element values from rolls with the observed count
# mean   = actual average roll
# ideal_mean = (low + high) / 2

rolls = []
File.foreach(ARGV[0], 'r') do |line|
  rolls.concat line.split.map{|_|_.to_i(10)}
end
low, high = rolls.minmax
puts "#{rolls.size} rolls from #{low} to #{high}"
faces = high - low + 1

counts = Array.new(high + 1)
sum = 0
rolls.each {|_|counts[_] = (counts[_] || 0) + 1; sum += _}
mean = sum.to_f / rolls.size

# mu = (((num_values + 1)*num_values)/2).to_f / num_values
variance = rolls.inject(0.0) {|sum,x| sum + (x.to_f - mean)**2 } / rolls.size

# assuming that all values will have been seen
buckets = counts.compact.size

# The actual number in the cell will be approximately normally distributed
# with mean equal to the expected number in each cell and standard deviation
# approximately the square root of the mean.
expected = rolls.size.to_f / buckets
exp_std_dev = Math.sqrt(expected)
puts "Expected per bucket: %3.1f (σ ≈ %4.2f)"%[expected, exp_std_dev]
low_bucket, high_bucket = counts.compact.minmax
(low_bucket..high_bucket).each do |count|
  if count <= expected - exp_std_dev
    go :red_on_black
  elsif count <= expected
    go :cyan_on_black
  elsif count <= expected + exp_std_dev
    go :yellow_on_black
  else
    go :red_on_black
  end
  print "%3d: "%[count]
  counts.map.with_index.select{|(c,i)| c==count }.map{|(c,i)| i }.each do |pos|
    print "(%2d)"%[pos]
  end
  go :white_on_black
  puts
end

puts "Mean     (μ) : %5.2f"%[mean]
puts "Variance (σ²): %5.2f"%[variance]
puts "Std.Dev  (σ) : %7.4f"%[Math.sqrt(variance)]
unless buckets == faces
  go :black_on_red
  puts " bad rolls, not all values from #{low} to #{high} "
  missing = (low..high).to_a - counts.map_with_index{|c,i|i if c}.compact
  puts "  missing: #{missing.inspect} "
  go :white_on_black

  go :red_on_black
  puts "Observed #{buckets} bucket#{'s' unless buckets == 1}, but expected #{faces}"
  go :white_on_black
end

# ------------------------------------------------------------------------------
# Histogram tests

Histogram = {
  # d#     5%    1%
  4  => [ 2.49, 3.02 ],         # Due to being in the
  6  => [ 2.63, 3.14 ],         # extreme tails of the
  8  => [ 2.73, 3.22 ],         # distribution, combined
  10 => [ 2.80, 3.29 ],         # with slight asymmetry,
  12 => [ 2.86, 3.34 ],         # the ranges we get are
  20 => [ 3.02, 3.48 ],         # sometimes out a bit.
}

if Histogram.has_key?(faces)
  hist_factor = Math.sqrt(expected * (faces - 1) / faces)
  hist_lines = Histogram[faces].reverse.map{|_|-1*_} + Histogram[faces]
  hist_lines.map! {|x|
    (expected + x * hist_factor).ceil
  }
  hist_lines.unshift 0
  puts "lines at: #{hist_lines.inspect}"
else
  hist_lines = [0,0,0,expected.ceil,rolls.size]
  go :magenta_on_black
  puts "Don't know how to interpret these results for a d%d"%[faces]
end
go :white_on_black

hist_colors = [ :magenta_on_black, # outside 99%
                :blue_on_black,    # outside 95%
                :green_on_black,   #  OK!
                :yellow_on_black,  # outside 95%
                :red_on_black,     # outside 99%
              ]

# ------------------------------------------------------------------------------
# The χ² test

ChiSquaredDistribution = {
  # d# =>   5%     1%             df
  4  => [  7.81, 11.34, ],      #  3
  6  => [ 11.07, 15.09, ],      #  5
  8  => [ 14.07, 18.48, ],      #  7
  10 => [ 16.92, 21.67, ],      #  9
  12 => [ 19.68, 24.72, ],      # 11
  20 => [ 30.14, 36.19, ],      # 19
}

chi_squared = 0.0

go :yellow_on_black
puts "Histogram:"
counts.each.with_index do |count, num|
  next if count.nil?
  chi_squared += ((count.to_f - expected)**2)/expected

  go :white_on_black
  print "%3d: "%num
  go hist_colors[0]
  count.times do |c|
    if c < hist_lines[1]
      go hist_colors[0]
    elsif c < hist_lines[2]
      go hist_colors[1]
    elsif c < hist_lines[3]
      go hist_colors[2]
    elsif c < hist_lines[4]
      go hist_colors[3]
    else
      go hist_colors[4]
    end
    print '*'
  end
  puts
end
go :white_on_black

if ChiSquaredDistribution.has_key?(faces)
  puts "For your d%d, the χ² value is %7.4f"%[faces, chi_squared]
  if chi_squared > ChiSquaredDistribution[faces].last
    go :red_on_black
    puts "This exceeds %5.2f and is probably biased"%[ChiSquaredDistribution[faces].last]
  elsif chi_squared > ChiSquaredDistribution[faces].first
    go :yellow_on_black
    puts "This is between %5.2f and %5.2f and may be biased"%ChiSquaredDistribution[faces]
  else
    go :green_on_black
    puts "This is less than %5.2f and is probably fair"%[ChiSquaredDistribution[faces].first]
  end
else
  go :magenta_on_black
  puts "Don't know how to interpret these results for a d%d"%[faces]
end
go :white_on_black

# ------------------------------------------------------------------------------
# The Kolmogorov-Smirnov test:

sum_counts = 0
differences  = Array.new(counts.size)
(low..high).each.with_index do |roll, i|
  sum_counts += counts[roll]
  differences[roll] = (sum_counts - (i+1)*expected).abs
end
raise "Kolmogorov-Smirnov bad math? #{differences.inspect}" unless differences[high].zero?

diff_max = differences.compact.max
d_stat = diff_max.to_f / Math.sqrt(rolls.size)

KolmogorovSmirnov = {
  # d# =>  5%    1%
  4  => [ 1.08, 1.35 ],         #   These values apply pretty well
  6  => [ 1.10, 1.37 ],         #   irrespective of the total number of
  8  => [ 1.11, 1.38 ],         #   rolls, but I would use at least 10 rolls
  10 => [ 1.12, 1.39 ],         #   per face.  Note also that these values
  12 => [ 1.12, 1.40 ],         #   come from simulation, and are hence not
  20 => [ 1.14, 1.42 ],         #   exact. This doesn't really matter.
}

if KolmogorovSmirnov.has_key?(faces)
  puts "For your d%d, the Kolmogorov-Smirnov value is %4.2f"%[faces, d_stat]
  if d_stat > KolmogorovSmirnov[faces].last
    go :red_on_black
    puts "This exceeds %5.2f and is probably biased"%KolmogorovSmirnov[faces].last
  elsif d_stat > KolmogorovSmirnov[faces].first
    go :yellow_on_black
    puts "This is between %5.2f and %5.2f and may be biased"%KolmogorovSmirnov[faces]
  else
    go :green_on_black
    puts "This is less than %5.2f and is probably fair"%[KolmogorovSmirnov[faces].first]
  end
else
  go :magenta_on_black
  puts "Don't know how to interpret these results for a d%d"%[faces]
end
go :white_on_black

# ------------------------------------------------------------------------------
# Histogram sorted by counts

go :cyan_on_black
puts "By count:"
counts.each.with_index.select{|count,_|count}.sort.each do |(count, num)|
  puts "%3d: %s"%[num, '*'*count]
end
go :white_on_black

if faces == 20
  rab_faces = [ [ 1, ],
                [ 7, 19, 13, ],
                [ 15,17, 3,9, 11,5, ],
                [ 12,10, 16,6, 4,18, ],
                [ 8, 14, 2, ],
                [ 20, ],
              ]
  rab_faces.each do |layer|
    str = layer.map {|value| "#{value}:#{counts[value]}" } * ' '
    puts str.center(36)
  end
end

__END__
=begin
References:

[1] http://wiretap.area.com/Gopher/Library/Article/Gaming/fairdice.txt
    Newsgroups: rec.games.frp.archives
    From: barnett@agsm.unsw.oz.au (Glen Barnett)
    Subject: PAPER: Testing Dice for bias
    Message-ID: <al#23np@rpi.edu>
    Organization: The Australian Graduate School of Management
    Date: Wed, 2 Dec 1992 04:41:08 GMT
    Approved: goldm@rpi.edu
    Lines: 749

    The following information is intended for distribution over Internet,
    and outside of that may be copied for personal use only.
    (c) Glen L. Barnett, 1992. All rights reserved.

-------------------------------------------------------------------

Conover, W.J. (1980): Practical nonparametric statistics,
2nd Ed., Wiley, New York.

Neave, H.R. and Worthington, P.L.B. (1988): Distribution-free tests,
Unwin Hyman, London.
=end