Probability of the sum of n dice rolls with three different methods

sum_of_rolling_dice.py

import math, numpy as np
from math import factorial as fact
import statistics as stats

def prob_of_sum_gaussian(s,n):
    mean = 3.5*n
    std = np.sqrt(n*105/36)
    normal_dist = stats.NormalDist(mean, std)
    return normal_dist.pdf(s)

def prob_of_sum_combinatorics(s,n):
    prob = 0
    k = np.floor((s - n)/6).astype(int)
    for i in range(0,k + 1):
        prob += (-1)**(i) * math.comb(n,i) * math.comb(s - 6*i - 1,n - 1)
    return np.divide(prob,6**n)

def prob_of_sum_experimental(s,n, throws = 1_000_000):
    sum_throws = simulate_throwing_n_dice(n, throws)
    return np.divide((sum_throws == s).sum(),throws)

def simulate_throwing_n_dice(n_dice, t_throws):
    throws = np.random.randint(1,7, size = (t_throws, n_dice))
    sum_throws = throws.sum(axis = 1)
    return sum_throws

s = 19
n = 4
p_exp   = prob_of_sum_experimental(s,n, throws = 1_000_000)
p_comb  = prob_of_sum_combinatorics(s,n)
p_gauss = prob_of_sum_gaussian(s, n)
print(p_exp)
print(p_comb)
print(p_gauss)