Enable Javascript in your browser and then refresh this page, for a much enhanced experience.
6-liner: scipy.integrate and probabilities power solution in 3rd party category for Probably Dice by Phil15
""" D1, ..., Dn n dices with equivalent probability for 1, 2, ..., sides: 1/sides.
What is the probability distribution of D1 + ... + Dn ?
It is equivalent to know the characteristic function 'phi'.
https://en.wikipedia.org/wiki/Characteristic_function_(probability_theory)
We have this, by definition:
phi(x) = sum(exp(i*k*x)*P(D1+...+Dn = k) for k in range(n, s*n+1))
We have this too, with a little calculus:
phi(x) = E(exp(i(D1+...+Dn)x)) = ...
= exp(i*n*x) * ((1 - exp(i*s*x)) / (1 - exp(i*x)) / s)**n
It's because we have:
E(exp(i*D1*x)) = exp(i*x) * (1 - exp(i*s*x)) / (1 - exp(i*x)) / s
and
E(exp(i(D1+...+Dn)x)) = ( E(exp(i*D1*x)) ) ** n
Tricky part:
So P(D1+...+Dn=t) = Real_part( integral(phi(x)*exp(-i*t*x), 0, pi) / pi )
Because Real_part( integral(exp(i*k*x), 0, pi) ) == 0 if k!=0.
We take real part then compute the integral.
f(x) = (phi(x) * exp(-i*t*x)).real """
from scipy import integrate
from cmath import exp, pi
ei = lambda x: exp(1j * x)
def probability(n, s, t): # dice_Number, Sides, Target
f = lambda x: (ei((n-t)*x) * ((1 - ei(s*x)) / (1 - ei(x)) / s)**n).real
return integrate.quad(f, 0, pi)[0] / pi #[0] give result, [1] an abs error.
Nov. 6, 2018
