# Untitled

By: a guest on May 23rd, 2012  |  syntax: None  |  size: 8.11 KB  |  hits: 14  |  expires: Never
Text below is selected. Please press Ctrl+C to copy to your clipboard. (⌘+C on Mac)
1. # -*- coding: utf-8 -*-
2. begin
3.   require 'ansi'
4.   def go(color)
5.     print ANSI.send(color)
6.   end
8.   \$stderr.puts "to get colored output: gem install ansi"
9.   def go(color); nil; end
10. end
11.
12.
13. # rolls  = an array holding the individual results
14. # low    = minimum element value in rolls; for a dN die, this should be 1
15. # high   = maximum element value in rolls; for a dN die, this should be N
16. # counts = an array indexed by the element values from rolls with the observed count
17. # mean   = actual average roll
18. # ideal_mean = (low + high) / 2
19.
20. rolls = []
21. File.foreach(ARGV[0], 'r') do |line|
22.   rolls.concat line.split.map{|_|_.to_i(10)}
23. end
24. low, high = rolls.minmax
25. puts "#{rolls.size} rolls from #{low} to #{high}"
26. faces = high - low + 1
27.
28. counts = Array.new(high + 1)
29. sum = 0
30. rolls.each {|_|counts[_] = (counts[_] || 0) + 1; sum += _}
31. mean = sum.to_f / rolls.size
32.
33. # mu = (((num_values + 1)*num_values)/2).to_f / num_values
34. variance = rolls.inject(0.0) {|sum,x| sum + (x.to_f - mean)**2 } / rolls.size
35.
36. # assuming that all values will have been seen
37. buckets = counts.compact.size
38.
39. # The actual number in the cell will be approximately normally distributed
40. # with mean equal to the expected number in each cell and standard deviation
41. # approximately the square root of the mean.
42. expected = rolls.size.to_f / buckets
43. exp_std_dev = Math.sqrt(expected)
44. puts "Expected per bucket: %3.1f (σ ≈ %4.2f)"%[expected, exp_std_dev]
45. low_bucket, high_bucket = counts.compact.minmax
46. (low_bucket..high_bucket).each do |count|
47.   if count <= expected - exp_std_dev
48.     go :red_on_black
49.   elsif count <= expected
50.     go :cyan_on_black
51.   elsif count <= expected + exp_std_dev
52.     go :yellow_on_black
53.   else
54.     go :red_on_black
55.   end
56.   print "%3d: "%[count]
57.   counts.map.with_index.select{|(c,i)| c==count }.map{|(c,i)| i }.each do |pos|
58.     print "(%2d)"%[pos]
59.   end
60.   go :white_on_black
61.   puts
62. end
63.
64. puts "Mean     (μ) : %5.2f"%[mean]
65. puts "Variance (σ²): %5.2f"%[variance]
66. puts "Std.Dev  (σ) : %7.4f"%[Math.sqrt(variance)]
67. unless buckets == faces
68.   go :black_on_red
69.   puts " bad rolls, not all values from #{low} to #{high} "
70.   missing = (low..high).to_a - counts.map_with_index{|c,i|i if c}.compact
71.   puts "  missing: #{missing.inspect} "
72.   go :white_on_black
73.
74.   go :red_on_black
75.   puts "Observed #{buckets} bucket#{'s' unless buckets == 1}, but expected #{faces}"
76.   go :white_on_black
77. end
78.
79. # ------------------------------------------------------------------------------
80. # Histogram tests
81.
82. Histogram = {
83.   # d#     5%    1%
84.   4  => [ 2.49, 3.02 ],         # Due to being in the
85.   6  => [ 2.63, 3.14 ],         # extreme tails of the
86.   8  => [ 2.73, 3.22 ],         # distribution, combined
87.   10 => [ 2.80, 3.29 ],         # with slight asymmetry,
88.   12 => [ 2.86, 3.34 ],         # the ranges we get are
89.   20 => [ 3.02, 3.48 ],         # sometimes out a bit.
90. }
91.
92. if Histogram.has_key?(faces)
93.   hist_factor = Math.sqrt(expected * (faces - 1) / faces)
94.   hist_lines = Histogram[faces].reverse.map{|_|-1*_} + Histogram[faces]
95.   hist_lines.map! {|x|
96.     (expected + x * hist_factor).ceil
97.   }
98.   hist_lines.unshift 0
99.   puts "lines at: #{hist_lines.inspect}"
100. else
101.   hist_lines = [0,0,0,expected.ceil,rolls.size]
102.   go :magenta_on_black
103.   puts "Don't know how to interpret these results for a d%d"%[faces]
104. end
105. go :white_on_black
106.
107. hist_colors = [ :magenta_on_black, # outside 99%
108.                 :blue_on_black,    # outside 95%
109.                 :green_on_black,   #  OK!
110.                 :yellow_on_black,  # outside 95%
111.                 :red_on_black,     # outside 99%
112.               ]
113.
114. # ------------------------------------------------------------------------------
115. # The χ² test
116.
117. ChiSquaredDistribution = {
118.   # d# =>   5%     1%             df
119.   4  => [  7.81, 11.34, ],      #  3
120.   6  => [ 11.07, 15.09, ],      #  5
121.   8  => [ 14.07, 18.48, ],      #  7
122.   10 => [ 16.92, 21.67, ],      #  9
123.   12 => [ 19.68, 24.72, ],      # 11
124.   20 => [ 30.14, 36.19, ],      # 19
125. }
126.
127. chi_squared = 0.0
128.
129. go :yellow_on_black
130. puts "Histogram:"
131. counts.each.with_index do |count, num|
132.   next if count.nil?
133.   chi_squared += ((count.to_f - expected)**2)/expected
134.
135.   go :white_on_black
136.   print "%3d: "%num
137.   go hist_colors[0]
138.   count.times do |c|
139.     if c < hist_lines[1]
140.       go hist_colors[0]
141.     elsif c < hist_lines[2]
142.       go hist_colors[1]
143.     elsif c < hist_lines[3]
144.       go hist_colors[2]
145.     elsif c < hist_lines[4]
146.       go hist_colors[3]
147.     else
148.       go hist_colors[4]
149.     end
150.     print '*'
151.   end
152.   puts
153. end
154. go :white_on_black
155.
156. if ChiSquaredDistribution.has_key?(faces)
157.   puts "For your d%d, the χ² value is %7.4f"%[faces, chi_squared]
158.   if chi_squared > ChiSquaredDistribution[faces].last
159.     go :red_on_black
160.     puts "This exceeds %5.2f and is probably biased"%[ChiSquaredDistribution[faces].last]
161.   elsif chi_squared > ChiSquaredDistribution[faces].first
162.     go :yellow_on_black
163.     puts "This is between %5.2f and %5.2f and may be biased"%ChiSquaredDistribution[faces]
164.   else
165.     go :green_on_black
166.     puts "This is less than %5.2f and is probably fair"%[ChiSquaredDistribution[faces].first]
167.   end
168. else
169.   go :magenta_on_black
170.   puts "Don't know how to interpret these results for a d%d"%[faces]
171. end
172. go :white_on_black
173.
174. # ------------------------------------------------------------------------------
175. # The Kolmogorov-Smirnov test:
176.
177. sum_counts = 0
178. differences  = Array.new(counts.size)
179. (low..high).each.with_index do |roll, i|
180.   sum_counts += counts[roll]
181.   differences[roll] = (sum_counts - (i+1)*expected).abs
182. end
183. raise "Kolmogorov-Smirnov bad math? #{differences.inspect}" unless differences[high].zero?
184.
185. diff_max = differences.compact.max
186. d_stat = diff_max.to_f / Math.sqrt(rolls.size)
187.
188. KolmogorovSmirnov = {
189.   # d# =>  5%    1%
190.   4  => [ 1.08, 1.35 ],         #   These values apply pretty well
191.   6  => [ 1.10, 1.37 ],         #   irrespective of the total number of
192.   8  => [ 1.11, 1.38 ],         #   rolls, but I would use at least 10 rolls
193.   10 => [ 1.12, 1.39 ],         #   per face.  Note also that these values
194.   12 => [ 1.12, 1.40 ],         #   come from simulation, and are hence not
195.   20 => [ 1.14, 1.42 ],         #   exact. This doesn't really matter.
196. }
197.
198. if KolmogorovSmirnov.has_key?(faces)
199.   puts "For your d%d, the Kolmogorov-Smirnov value is %4.2f"%[faces, d_stat]
200.   if d_stat > KolmogorovSmirnov[faces].last
201.     go :red_on_black
202.     puts "This exceeds %5.2f and is probably biased"%KolmogorovSmirnov[faces].last
203.   elsif d_stat > KolmogorovSmirnov[faces].first
204.     go :yellow_on_black
205.     puts "This is between %5.2f and %5.2f and may be biased"%KolmogorovSmirnov[faces]
206.   else
207.     go :green_on_black
208.     puts "This is less than %5.2f and is probably fair"%[KolmogorovSmirnov[faces].first]
209.   end
210. else
211.   go :magenta_on_black
212.   puts "Don't know how to interpret these results for a d%d"%[faces]
213. end
214. go :white_on_black
215.
216. # ------------------------------------------------------------------------------
217. # Histogram sorted by counts
218.
219. go :cyan_on_black
220. puts "By count:"
221. counts.each.with_index.select{|count,_|count}.sort.each do |(count, num)|
222.   puts "%3d: %s"%[num, '*'*count]
223. end
224. go :white_on_black
225.
226. if faces == 20
227.   rab_faces = [ [ 1, ],
228.                 [ 7, 19, 13, ],
229.                 [ 15,17, 3,9, 11,5, ],
230.                 [ 12,10, 16,6, 4,18, ],
231.                 [ 8, 14, 2, ],
232.                 [ 20, ],
233.               ]
234.   rab_faces.each do |layer|
235.     str = layer.map {|value| "#{value}:#{counts[value]}" } * ' '
236.     puts str.center(36)
237.   end
238. end
239.
240. __END__
241. =begin
242. References:
243.
244. [1] http://wiretap.area.com/Gopher/Library/Article/Gaming/fairdice.txt
245.     Newsgroups: rec.games.frp.archives
246.     From: barnett@agsm.unsw.oz.au (Glen Barnett)
247.     Subject: PAPER: Testing Dice for bias
248.     Message-ID: <al#23np@rpi.edu>
249.     Organization: The Australian Graduate School of Management
250.     Date: Wed, 2 Dec 1992 04:41:08 GMT
251.     Approved: goldm@rpi.edu
252.     Lines: 749
253.
254.     The following information is intended for distribution over Internet,
255.     and outside of that may be copied for personal use only.