Simulation: Dice

# A number (k) of s-sided dice are cast; the sum of the dice values is calculated. # The experiment is repeated n times and the distribution of the sums is calculated. # The distribution is presented graphically with the expected value, variance and # standard deviation also displayed. # The effect of the central limit theorem is very visible: The distribution converges to the normal distribution. # Lauri Ruotsalainen, 2010 import random from sage.plot.bar_chart import BarChart @interact def dice( k = slider(1, 10, 1, default=3, label="Number of dice (K)"), s = slider(2, 20, 1, default=6, label="Greatest die value (S)"), n = input_box(default=10000, label="Number of tosses (N)"), normal_dist = ("Show normal distribution density function", false) ): # Cast dice and add them into the dictionary of seen outcomes. outcomes = {} for i in range(n): Sum = 0 for j in range(k): Sum += random.randint(1, s) outcomes[Sum] = outcomes.get(Sum, 0) + 1 # Create an index list of all possible sums and calculate the relative frequencies. indices = range(k, s*k+1) # For the bar chart, displace the index values by -0.25 to position the bars correctly. displaced_indices = [c - 0.25 for c in indices] relative_heights = [outcomes.get(i, 0) / n for i in indices] # Draw the bar chart. G = Graphics() G.add_primitive(BarChart(displaced_indices, relative_heights, {"width":0.5, "rgbcolor":(0,0,1)})) # Calculate the variance and expected value of the random variable. expected_value = sum([1, 2, .. , s])/s*k i = var("i") variance = (sum(i^2, i, 1, s)/s - (sum([1, 2, .. , s])/s)^2)*k # Draw the normal distribution if required. if normal_dist: normal_dist_density(E,V,x) = 1/sqrt(2*pi*V)*e^(-(x-E)^2/(2*V)) G += plot(normal_dist_density(expected_value, variance, x), xmin=k, xmax=k*s) # Show the distribution with the expected value, variance and standard deviation. G.show(xmin=k, xmax=s*k, ymin=0) html("$\\text{Expected Value = %s}$" %(N(expected_value, digits=6))) html("$\\text{Variance = %s}$" %(N(variance, digits=6))) html("$\\text{Standard Deviation = %s}$" %(N(sqrt(variance), digits=6)))