Skip to content

Instantly share code, notes, and snippets.

@deanmalmgren

deanmalmgren/README.md

Last active August 29, 2015 14:11
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save deanmalmgren/d33d4e56a5f1dc9ed2cc to your computer and use it in GitHub Desktop.
Save deanmalmgren/d33d4e56a5f1dc9ed2cc to your computer and use it in GitHub Desktop.
Download ZIP
elongated dice analysis
Raw
README.md

This gist has some quick analysis related to the outcomes from an elongated dice experiment. analyze_roll_data.py is a script to conduct the analysis. p16_aspect_ratio.agr contains the results of this analysis in an xmgrace file.

We tried a few quick things to work out a theoretical relationship between the aspect ratio and the probability of a 1 or a 6 being rolled, which is contained in theoretical_curve.py.

One idea we had was to assume that, just prior to the elongated die coming to a final resting state, the die would be at a random angle with the surface of the table. The probability of the die landing 1 or a 6 is therefore related to the probability that the center of mass is at an angle greater than the angle with the table. Its possible we didn't do this correctly, but our calculations did not match very well with the data (see theoretical_curve.py if you need convincing), which leads us to believe that the mechanics of how these dice roll actually plays a large part in biasing the angle between the die and the surface of a table just prior to landing.

Raw
analyze_roll_data.py
"""Run this script to create a plot of the likelihood that each die is rolled with a 1 or a 6.
"""
import csv
import collections
import math
class Die(object):
def __init__(self, label):
self.label = label
self.aspect_ratio = None
self.rolls = []
def p16(self):
"""calculate the probability of rolling a 1 or a 6"""
# count up the number of rolls of each die, and be sure to add a
# pseudocount for each side of the die to properly account for prior
counter = collections.Counter(range(1,7))
counter.update(self.rolls)
# calculate the probability
return float(counter[1]+counter[6]) / sum(counter.values())
def p2345(self):
return 1.0 - self.p16()
def stderr_p16(self):
"""calculate the error bars on p16"""
if self.rolls:
return math.sqrt(self.p16()*(1-self.p16())/len(self.rolls))
return None
# function to convert aspect ratios to floats
def float_or_null(value):
if value:
return float(value)
return None
# instantiate a bunch of dice
dice_list = []
for o in range(ord('A'), ord('Z')+1):
label = chr(o)
dice_list.append(Die(label))
# read in the roll data from the master spreadsheet
#
# TODO: should someday just access directly from google docs
# https://docs.google.com/spreadsheets/d/1cG_XIIq8cH_5sK4FhPIG7iHbP0WNb_e5Rn8tTJCed68/edit#gid=0
with open('roll_data.csv') as stream:
reader = csv.reader(stream)
# set the aspect ratio on each die
aspect_ratio_list = map(float_or_null,reader.next())
for aspect_ratio, die in zip(aspect_ratio_list, dice_list):
die.aspect_ratio = aspect_ratio
# add the rolls
for row in reader:
for roll, die in zip(row, dice_list):
if roll:
die.rolls.append(int(roll))
# for each die, calculate the probability of rolling a 1 or 6 versus a 2-5
for die in dice_list:
print die.label, die.aspect_ratio, die.p16(), die.stderr_p16()
Raw
p16_aspect_ratio.agr
# Grace project file
#
@version 50123
@page size 792, 612
@page scroll 5%
@page inout 5%
@link page off
@map font 0 to "Times-Roman", "Times-Roman"
@map font 1 to "Times-Italic", "Times-Italic"
@map font 2 to "Times-Bold", "Times-Bold"
@map font 3 to "Times-BoldItalic", "Times-BoldItalic"
@map font 4 to "Helvetica", "Helvetica"
@map font 5 to "Helvetica-Oblique", "Helvetica-Oblique"
@map font 6 to "Helvetica-Bold", "Helvetica-Bold"
@map font 7 to "Helvetica-BoldOblique", "Helvetica-BoldOblique"
@map font 8 to "Courier", "Courier"
@map font 9 to "Courier-Oblique", "Courier-Oblique"
@map font 10 to "Courier-Bold", "Courier-Bold"
@map font 11 to "Courier-BoldOblique", "Courier-BoldOblique"
@map font 12 to "Symbol", "Symbol"
@map font 13 to "ZapfDingbats", "ZapfDingbats"
@map color 0 to (255, 255, 255), "white"
@map color 1 to (0, 0, 0), "black"
@map color 2 to (228, 26, 28), "Set1-1"
@map color 3 to (55, 126, 184), "Set1-2"
@map color 4 to (77, 175, 74), "Set1-3"
@map color 5 to (152, 78, 163), "Set1-4"
@map color 6 to (255, 127, 0), "Set1-5"
@map color 7 to (255, 255, 51), "Set1-6"
@map color 8 to (166, 86, 40), "Set1-7"
@map color 9 to (247, 129, 191), "Set1-8"
@map color 10 to (153, 153, 153), "Set1-9"
@map color 11 to (102, 194, 165), "Set2-1"
@map color 12 to (252, 141, 98), "Set2-2"
@map color 13 to (141, 160, 203), "Set2-3"
@map color 14 to (231, 138, 195), "Set2-4"
@map color 15 to (166, 216, 84), "Set2-5"
@map color 16 to (255, 217, 47), "Set2-6"
@map color 17 to (229, 196, 148), "Set2-7"
@map color 18 to (179, 179, 179), "Set2-8"
@reference date 0
@date wrap off
@date wrap year 1950
@default linewidth 2.0
@default linestyle 1
@default color 1
@default pattern 1
@default font 0
@default char size 1.000000
@default symbol size 1.000000
@default sformat "%.8g"
@background color 0
@page background fill on
@timestamp off
@timestamp 0.03, 0.03
@timestamp color 1
@timestamp rot 0
@timestamp font 0
@timestamp char size 1.000000
@timestamp def "Sat Dec 13 04:55:38 2014"
@with string
@ string on
@ string loctype view
@ string 0.97, 0.45
@ string color 10
@ string rot 0
@ string font 6
@ string just 0
@ string char size 1.000000
@ string def "fair die"
@r0 off
@link r0 to g0
@r0 type above
@r0 linestyle 1
@r0 linewidth 1.0
@r0 color 1
@r0 line 0, 0, 0, 0
@r1 off
@link r1 to g0
@r1 type above
@r1 linestyle 1
@r1 linewidth 1.0
@r1 color 1
@r1 line 0, 0, 0, 0
@r2 off
@link r2 to g0
@r2 type above
@r2 linestyle 1
@r2 linewidth 1.0
@r2 color 1
@r2 line 0, 0, 0, 0
@r3 off
@link r3 to g0
@r3 type above
@r3 linestyle 1
@r3 linewidth 1.0
@r3 color 1
@r3 line 0, 0, 0, 0
@r4 off
@link r4 to g0
@r4 type above
@r4 linestyle 1
@r4 linewidth 1.0
@r4 color 1
@r4 line 0, 0, 0, 0
@g0 on
@g0 hidden false
@g0 type XY
@g0 stacked false
@g0 bar hgap 0.000000
@g0 fixedpoint off
@g0 fixedpoint type 0
@g0 fixedpoint xy 0.000000, 0.000000
@g0 fixedpoint format general general
@g0 fixedpoint prec 6, 6
@with g0
@ world 0, 0, 4, 1
@ stack world 0, 0, 0, 0
@ znorm 1
@ view 0.200000, 0.200000, 1.100000, 0.900000
@ title ""
@ title font 4
@ title size 1.500000
@ title color 1
@ subtitle ""
@ subtitle font 0
@ subtitle size 1.000000
@ subtitle color 1
@ xaxes scale Normal
@ yaxes scale Normal
@ xaxes invert off
@ yaxes invert off
@ xaxis on
@ xaxis type zero false
@ xaxis offset 0.000000 , 0.000000
@ xaxis bar on
@ xaxis bar color 1
@ xaxis bar linestyle 1
@ xaxis bar linewidth 2.0
@ xaxis label "Aspect ratio"
@ xaxis label layout para
@ xaxis label place auto
@ xaxis label char size 2.000000
@ xaxis label font 4
@ xaxis label color 1
@ xaxis label place normal
@ xaxis tick on
@ xaxis tick major 1
@ xaxis tick minor ticks 1
@ xaxis tick default 6
@ xaxis tick place rounded true
@ xaxis tick in
@ xaxis tick major size 1.000000
@ xaxis tick major color 1
@ xaxis tick major linewidth 2.0
@ xaxis tick major linestyle 1
@ xaxis tick major grid off
@ xaxis tick minor color 1
@ xaxis tick minor linewidth 2.0
@ xaxis tick minor linestyle 1
@ xaxis tick minor grid off
@ xaxis tick minor size 0.500000
@ xaxis ticklabel on
@ xaxis ticklabel format general
@ xaxis ticklabel prec 5
@ xaxis ticklabel formula ""
@ xaxis ticklabel append ""
@ xaxis ticklabel prepend ""
@ xaxis ticklabel angle 0
@ xaxis ticklabel skip 0
@ xaxis ticklabel stagger 0
@ xaxis ticklabel place normal
@ xaxis ticklabel offset auto
@ xaxis ticklabel offset 0.000000 , 0.010000
@ xaxis ticklabel start type auto
@ xaxis ticklabel start 0.000000
@ xaxis ticklabel stop type auto
@ xaxis ticklabel stop 0.000000
@ xaxis ticklabel char size 1.500000
@ xaxis ticklabel font 4
@ xaxis ticklabel color 1
@ xaxis tick place both
@ xaxis tick spec type none
@ yaxis on
@ yaxis type zero false
@ yaxis offset 0.000000 , 0.000000
@ yaxis bar on
@ yaxis bar color 1
@ yaxis bar linestyle 1
@ yaxis bar linewidth 2.0
@ yaxis label "Probability of rolling 1 or 6"
@ yaxis label layout para
@ yaxis label place auto
@ yaxis label char size 2.000000
@ yaxis label font 4
@ yaxis label color 1
@ yaxis label place normal
@ yaxis tick on
@ yaxis tick major 0.2
@ yaxis tick minor ticks 1
@ yaxis tick default 6
@ yaxis tick place rounded true
@ yaxis tick in
@ yaxis tick major size 1.000000
@ yaxis tick major color 1
@ yaxis tick major linewidth 2.0
@ yaxis tick major linestyle 1
@ yaxis tick major grid off
@ yaxis tick minor color 1
@ yaxis tick minor linewidth 2.0
@ yaxis tick minor linestyle 1
@ yaxis tick minor grid off
@ yaxis tick minor size 0.500000
@ yaxis ticklabel on
@ yaxis ticklabel format decimal
@ yaxis ticklabel prec 1
@ yaxis ticklabel formula ""
@ yaxis ticklabel append ""
@ yaxis ticklabel prepend ""
@ yaxis ticklabel angle 0
@ yaxis ticklabel skip 0
@ yaxis ticklabel stagger 0
@ yaxis ticklabel place normal
@ yaxis ticklabel offset auto
@ yaxis ticklabel offset 0.000000 , 0.010000
@ yaxis ticklabel start type auto
@ yaxis ticklabel start 0.000000
@ yaxis ticklabel stop type auto
@ yaxis ticklabel stop 0.000000
@ yaxis ticklabel char size 1.500000
@ yaxis ticklabel font 4
@ yaxis ticklabel color 1
@ yaxis tick place both
@ yaxis tick spec type none
@ altxaxis off
@ altyaxis off
@ legend on
@ legend loctype view
@ legend 0.85, 0.8
@ legend box color 1
@ legend box pattern 1
@ legend box linewidth 2.0
@ legend box linestyle 0
@ legend box fill color 0
@ legend box fill pattern 1
@ legend font 4
@ legend char size 1.250000
@ legend color 1
@ legend length 4
@ legend vgap 1
@ legend hgap 1
@ legend invert false
@ frame type 0
@ frame linestyle 1
@ frame linewidth 2.0
@ frame color 1
@ frame pattern 1
@ frame background color 0
@ frame background pattern 0
@ s0 hidden false
@ s0 type xy
@ s0 symbol 0
@ s0 symbol size 1.000000
@ s0 symbol color 10
@ s0 symbol pattern 1
@ s0 symbol fill color 10
@ s0 symbol fill pattern 0
@ s0 symbol linewidth 2.0
@ s0 symbol linestyle 1
@ s0 symbol char 65
@ s0 symbol char font 0
@ s0 symbol skip 0
@ s0 line type 1
@ s0 line linestyle 3
@ s0 line linewidth 4.0
@ s0 line color 10
@ s0 line pattern 1
@ s0 baseline type 0
@ s0 baseline off
@ s0 dropline off
@ s0 fill type 0
@ s0 fill rule 0
@ s0 fill color 1
@ s0 fill pattern 1
@ s0 avalue off
@ s0 avalue type 2
@ s0 avalue char size 1.000000
@ s0 avalue font 0
@ s0 avalue color 1
@ s0 avalue rot 0
@ s0 avalue format general
@ s0 avalue prec 3
@ s0 avalue prepend ""
@ s0 avalue append ""
@ s0 avalue offset 0.000000 , 0.000000
@ s0 errorbar on
@ s0 errorbar place both
@ s0 errorbar color 10
@ s0 errorbar pattern 1
@ s0 errorbar size 1.000000
@ s0 errorbar linewidth 2.0
@ s0 errorbar linestyle 1
@ s0 errorbar riser linewidth 2.0
@ s0 errorbar riser linestyle 1
@ s0 errorbar riser clip off
@ s0 errorbar riser clip length 0.100000
@ s0 comment "gg"
@ s0 legend ""
@ s1 hidden false
@ s1 type xydy
@ s1 symbol 2
@ s1 symbol size 0.750000
@ s1 symbol color 1
@ s1 symbol pattern 1
@ s1 symbol fill color 1
@ s1 symbol fill pattern 1
@ s1 symbol linewidth 2.0
@ s1 symbol linestyle 1
@ s1 symbol char 65
@ s1 symbol char font 0
@ s1 symbol skip 0
@ s1 line type 0
@ s1 line linestyle 1
@ s1 line linewidth 2.0
@ s1 line color 1
@ s1 line pattern 1
@ s1 baseline type 0
@ s1 baseline off
@ s1 dropline off
@ s1 fill type 0
@ s1 fill rule 0
@ s1 fill color 1
@ s1 fill pattern 1
@ s1 avalue off
@ s1 avalue type 2
@ s1 avalue char size 1.000000
@ s1 avalue font 0
@ s1 avalue color 1
@ s1 avalue rot 0
@ s1 avalue format general
@ s1 avalue prec 3
@ s1 avalue prepend ""
@ s1 avalue append ""
@ s1 avalue offset 0.000000 , 0.000000
@ s1 errorbar on
@ s1 errorbar place both
@ s1 errorbar color 1
@ s1 errorbar pattern 1
@ s1 errorbar size 1.000000
@ s1 errorbar linewidth 2.0
@ s1 errorbar linestyle 1
@ s1 errorbar riser linewidth 2.0
@ s1 errorbar riser linestyle 1
@ s1 errorbar riser clip off
@ s1 errorbar riser clip length 0.100000
@ s1 comment "/Users/rdm/Dropbox (Datascope Analytics)/Moveo/dice_results/kk"
@ s1 legend ""
@ s2 hidden false
@ s2 type xydy
@ s2 symbol 2
@ s2 symbol size 0.750000
@ s2 symbol color 3
@ s2 symbol pattern 1
@ s2 symbol fill color 3
@ s2 symbol fill pattern 1
@ s2 symbol linewidth 2.0
@ s2 symbol linestyle 1
@ s2 symbol char 65
@ s2 symbol char font 0
@ s2 symbol skip 0
@ s2 line type 0
@ s2 line linestyle 1
@ s2 line linewidth 2.0
@ s2 line color 3
@ s2 line pattern 1
@ s2 baseline type 0
@ s2 baseline off
@ s2 dropline off
@ s2 fill type 0
@ s2 fill rule 0
@ s2 fill color 1
@ s2 fill pattern 1
@ s2 avalue on
@ s2 avalue type 4
@ s2 avalue char size 1.000000
@ s2 avalue font 4
@ s2 avalue color 3
@ s2 avalue rot 0
@ s2 avalue format general
@ s2 avalue prec 3
@ s2 avalue prepend "Die "
@ s2 avalue append ""
@ s2 avalue offset 0.050000 , 0.010000
@ s2 errorbar on
@ s2 errorbar place both
@ s2 errorbar color 3
@ s2 errorbar pattern 1
@ s2 errorbar size 1.000000
@ s2 errorbar linewidth 2.0
@ s2 errorbar linestyle 1
@ s2 errorbar riser linewidth 2.0
@ s2 errorbar riser linestyle 1
@ s2 errorbar riser clip off
@ s2 errorbar riser clip length 0.100000
@ s2 comment "copy of set G0.S1"
@ s2 legend ""
@target G0.S0
@type xy
0 0.333333
4 0.333333
&
@target G0.S1
@type xydy
0.6875 0.67857143 0.066047293 "B"
1.3125 0.18867925 0.039125361 "C"
0.3125 0.93396226 0.024834805 "E"
2.5625 0.053571429 0.031843847 "F"
0.75 0.51515152 0.064520078 "G"
1.875 0.10714286 0.043740888 "H"
1.1875 0.23214286 0.059708048 "I"
1.5 0.16071429 0.051939427 "K"
0.75 0.55357143 0.070303642 "L"
0.9375 0.49056604 0.049991099 "M"
0.3125 0.85714286 0.049487166 "N"
2.25 0.41071429 0.069574142 "O"
0.875 0.035714286 0.026244533 "P"
1 0.46969697 0.064431065 "Q"
0.375 0.83928571 0.051939427 "R"
1.0625 0.33018868 0.047028089 "S"
0.5625 0.74528302 0.043570201 "T"
1.25 0.19642857 0.056186188 "U"
1.625 0.14285714 0.049487166 "V"
3.5 0.018867925 0.013605854 "W"
1.5 0.1969697 0.051344041 "X"
3.5625 0.03030303 0.022130204 "Z"
&
@target G0.S2
@type xydy
0.875 0.035714286 0.026244533 "P"
2.25 0.41071429 0.069574142 "O"
&
Raw
theoretical_curve.py
"""The idea for this came from @vlsd.
Imagine a pyramid that connects the center of mass with each of the four corners
on the square side of an elongated die. Just prior to coming to a final rest, if
the die is at a random angle, we know it will land on the square side if the
center of mass is above the square side.
I'm sure there's a way to theoretically calculate this (in fact, I think @vlsd
found one), but its pretty easy to simulate.
"""
import random
import math
def random_spherical_angle():
"""http://mathworld.wolfram.com/SpherePointPicking.html"""
u = random.random()
v = random.random()
theta = 2 * math.pi * u
phi = math.acos(2*v-1)
return theta, phi
def is_through_pyramid(aspect_ratio, theta, phi, w=1.0):
h = aspect_ratio * w
r = math.sqrt(0.25*h*h + 0.5*w*w)
x = r * math.cos(theta) * math.sin(phi)
y = r * math.sin(theta) * math.sin(phi)
z = r * math.cos(phi)
# without loss of generality, we can force z to always be positive
z = abs(z)
# need to figure out if the segment S that goes from the origin to the point
# at x,y,z also intersects the square Q that is located at z=h/2 with corners at +/-
# w/2
#
# start by finding the point of intersection of S with the plane that
# contains Q (e.g. z=h).
z_hat = 0.5*h
r_hat = z_hat / math.cos(phi)
x_hat = r_hat * math.cos(theta) * math.sin(phi)
y_hat = r_hat * math.sin(theta) * math.sin(phi)
# if the point x,y is inside of the square Q, return True!
return -0.5*w <= x_hat <= 0.5*w and -0.5*w <= y_hat <= 0.5*w
def estimate_p16(aspect_ratio):
p_16 = 0.0
n_trials = 10000
for n in range(n_trials):
theta, phi = random_spherical_angle()
if is_through_pyramid(aspect_ratio, theta, phi):
p_16 += 1
p_16 /= n_trials
return p_16
for a in range(300):
aspect_ratio = a * 0.01
print aspect_ratio, estimate_p16(aspect_ratio)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment