Distribution of the sum of $N$ loaded dice rolls I would like to calculate the probability distribution of the sum of all the faces of $N$ dice rolls. The face probabilities ${p_i}$ are know, but are not $1 \over 6$.
I have found answers for the case of a fair dice (i.e. $p_i={1 \over 6}$) here and here
For large $N$ I could apply the central limit theorem and use a normal distribution, but I don't know how to proceed for small $N$. (In particular, $N=2,4, 20$)
 A: It seems to me the spirit of the original context is experimental.
In order to get an analytic answer, even in a simple case with only two rolls, you need to specify precisely how the die is loaded.
Here is a simulation in R of an experiment with $n = 3$ rolls of a die loaded so that
faces 1 to 6 appear with probabilities $(1/12, 1/12, 1/12, 3/12, 3/12, 3/12)$. In practice, such a bias might be achieved by embedding
a heavy weight just below the corner where faces 1, 2, and 3 meet.
A million three-roll experiments are simulated. Each pass through
the loop simulates one three-roll experiment. The random variable
$T$ is the total on the three dice (values between 3 and 18, inclusive). 
For the mean and SD of the best-fitting normal distribution, you should be able to
find $E(T) = 12.75$ and $SD(T) = ?$, which are approximated in the simulation. (Also, the exact computation is shown for 
$E(T) = n\sum_{i=1}^6 ip_i.$)
 m = 10^6;  n = 3;  p = c(1, 1, 1, 3, 3, 3)/12
 t = numeric(m) # vector to receive totals
 for (i in 1:m) {
   faces = sample(1:6, n, repl=T, prob = p)
   t[i] = sum(faces)  }
 mean(t);  sd(t)
 ## 12.74981  # Approx E(T) = 12.75
 ## 2.657790  # Approx SD(T)
 hist(t, br=(2:18)+.5, prob=T, main="Sum of 3 Rolls", col="wheat", ylim=c(0,.15))
 curve(dnorm(x, mean(t), sd(t)), lwd=2, col="blue", add=T)
 3*sum((1:6)*p)
 ## 12.75     # Exact E(T)


The normal curve with matching mean and SD is not yet a good
fit for only three rolls. Perhaps much better for 20.
Using code $mutatis\;mutandis$: Yes 20 looks better.

A: As an alternative to Matthew Conroy's suggestion to use a computer algebra system, one can also code this with Python using numpy class of Polynomials. The convolution power can be simply obtained by calculating the power of a polynomial. The same example he suggests can be coded as follows:
In [14]: from numpy.polynomial.polynomial import Polynomial

In [15]: p=Polynomial((1/7, 1/7, 1/7, 1/7,1/7, 2/7))

In [16]: p**3  # or alternatively: np.power(p,3)
Out[16]: 
Polynomial([ 0.00291545,  0.00874636,  0.01749271,  0.02915452,  0.04373178,
        0.06997085,  0.09037901,  0.10495627,  0.11370262,  0.11661808,
        0.12244898,  0.09620991,  0.0728863 ,  0.05247813,  0.03498542,
        0.02332362], [-1.,  1.], [-1.,  1.])

Of course, one has to work with floats, but otherwise the results agree.
In numpy, the convolution can also be coded with np.convolve , but the application of successive convolutions is more cumbersome.
