Rolling a single dice repeatedly will result in a uniform distribution. But if we roll multiple dice the sum would be the Normal distribution due the Central Limit Theorem (CLT).
To verify this, simulated multiple dice rolls number of times to get a distribution and then used the Shapiro-Wilk test to calculate p-value.
What seems strange is that if we increase the number of rolls, the histogram becomes smoother and apparently closer to normal by p- value tends to decline until it stays below 0.05 for any number of repetitions.
Both Python and R gives the same results. Below are two pieces of code in Python. One plots p value vs. the number of rolls. The other repeats the experiment 1000 times and calculates the percentage of results supporting the null hypothesis with p > 0.05.
E.g., with 20 dice, the percentage for 100 rolls is about 94; 700 rolls - 75%, 1500 rolls - 11%, 2000 rolls - 0%.
How can we explain this phenomenon?
import numpy as np
import matplotlib.pyplot as plt
from scipy import stats
def shapiro_p(n_dice, n_throws):
rng = np.random.default_rng()
throwing = rng.integers(1, high=6, size=(n_throws, n_dice), endpoint=True)
X = np.sum(throwing, axis=1) / n_dice
shapiro_test = stats.shapiro(X)
return shapiro_test.pvalue
n_dice = int(input('Input the number of dice '))
n_max = int(input('Input the number of rolls '))
x = [n for n in range(3, n_max + 1)]
y = [shapiro_p(n_dice, n) for n in range(3, n_max + 1)]
figure_title = 'Number of dice: {0}'.format(n_dice)
plt.plot(x, y)
plt.hlines(0.05, 3, n_max + 1, colors='red')
plt.xlabel('Number of rolls')
plt.ylabel('p')
plt.title(figure_title)
plt.grid(True)
plt.show()
The second piece of code:
from scipy import stats
import random
dice = [i for i in range(1, 7)]
def shapiro_p(n_dice, n_rolls):
sampl = [sum(random.choices(dice, k=n_dice)) for j in range(n_rolls)]
shapiro_test = stats.shapiro(sampl)
return shapiro_test.pvalue
n_dice = int(input('Number of dice '))
n_rolls = int(input('Number of rolls in each test '))
# n_tests = int(input('Number of tests '))
n_tests = 1000
null_h = 0
for n in range(n_tests):
if shapiro_p(n_dice, n_rolls) >= 0.05:
null_h += 1
print('Percentage of tests supporting the null hypothesis: ', 100*null_h/n_tests)