#!/usr/bin/env python3

import numpy as np
import random
from scipy.stats import norm
from matplotlib import pyplot as plt , rc

########################################################################

# Here are the parameters for our little numerical experiment:

number_of_dice   = 3
sides_on_die     = 6
number_of_trials = 100000

########################################################################

# These two numbers are the mean and variance of the distribution
# function for rolling number_of_dice dice with sides_on_die sides
# and looking at the sum of all eyes.

mean = 
var  = 

########################################################################

# Here we do the actual work of creating the histogram

histogram  = np.zeros(number_of_dice*(sides_on_die-1)+1,int)
eye_values = np.arange(number_of_dice,number_of_dice*sides_on_die+1)

for trial in range(number_of_trials):
   eyes = 0
   for die in range(number_of_dice):
       eyes += int(random.random() * sides_on_die)
   histogram[eyes] += 1

########################################################################

#### Plotting the results and comparing to a Gaussian pdf

# first, just a few parameters to make the figure look nicer
font = {'family' : 'sans-serif',
       'weight' : 'normal',
       'size'   : 18}
rc('font', **font)
plt.figure(figsize=(8,6), dpi=80)

#...and now we can finally plot the estimated probability distribution:
plt.scatter(eye_values,histogram/number_of_trials,c='k',s=50,label="counts")

# here's the Gaussian distribution function for comparison
xx = np.linspace(number_of_dice-2,number_of_dice*sides_on_die+2,100)
plt.plot(xx, norm.pdf(xx,mean,np.sqrt(var)),'r-', lw=3, alpha=0.6,
label='Gaussian')

# and here are axis labels and plot title
plt.xlabel('total number of eyes, $E_N$')
plt.ylabel('probability, $P(E_N)$')
plt.title('Rolling a {}-die {} times. {} trials'
         .format(sides_on_die,number_of_dice,number_of_trials))
plt.legend()
plt.grid(True)
plt.tight_layout()   # this helps with fixing some cropping troubles
plt.show()

# Commenting out the previous line and uncommenting the
# next line will instead save the figure into a file:
#plt.savefig("rolling_dice.png",format="png",dpi=160)