elongated dice analysis

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.

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

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"
&

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)