Suppose I roll a 4-side dice 25 time to get a final sum between 25 and 100. How I calculated the distribution of probability for each sum between 25 and 100?
Solved – Multinomial distribution for 4 side dice roll
dicemultinomial-distributionprobability
Related Solutions
The relationship is simple: $S = Z_{(x-y+1):x} + ... + Z_{x:x}$. This should make logical sense: The sum of the top $k$ dice is the highest die plue the second highest die etc. down to the $k$th highest die. If you still believe $S \ne Z_{(x-y+1):x} + ... + Z_{x:x}$, try to exhibit a roll of the dice so that the sum of the top $k$ dice is not equal to the highest die plus the second highest etc. plus the $k$th highest.
Order statistics are only independent for trivial constant distributions. They are not independent here. Since the summands are not independent, you can't use TransformedDistribution in Mathematica. See the documentation which says "TransformedDistribution represents a transformed distribution where , , ... are independent and follow the distributions , , ...." This is why the distribution you calculate for the right hand side is not correct.
Because of the dependence, you can't determine the distribution of the sum from the distributions of the summands. The same is true in much simpler cases. If $X_0$ and $X_1$ are both $1$ with probability $1/2$ and $0$ with probability $1/2$, then it is possible that $X_0 = 1-X_1$ so that $X_0 + X_1$ is the constant $1$. It is also possible that $X_0 = X_1$ so that $X_0+X_1$ is $0$ with probability $1/2$ and $2$ with probability $1/2$. It is also possible that $X_0$ and $X_1$ are independent so that $X_0+X_1$ takes the values $0,1,2$ with probabilities $1/4,1/2,1/4$, respectively.
Nevertheless, $S = Z_{(x-y+1):x} + ... + Z_{x:x}$. Expectation is linear regardless of whether the terms are independent, so $E(S) = E(Z_{(x-y+1):x}) + ... + E(Z_{x:x})$.
I've thought about expressing the distribution of $S$, and I keep getting the same expression that techmologist did (with the corrected upper bound I edited in). Order statistics for discrete distributions are messy, so I don't expect to find a big simplification by using them.
The distribution of 2d8 is discrete triangular.
x 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
64.P(X=x) 1 2 3 4 5 6 7 8 7 6 5 4 3 2 1
% 1.56 3.13 4.69 6.25 7.81 9.38 10.9 12.5 10.9 9.38 7.81 6.25 4.69 3.13 1.56
If you need an algebraic expression for it, for 2d8 it's:
$p(x) = P(X=x) = \frac{1}{64} \min(x-1,17-x)$
As you add more dice, the cdf becomes closer and closer to a normal distribution, but if you want to use normal distributions to approximate probabilities for it, I'd suggest using a continuity correction. If you don't get too far into the tail, it should work pretty well for more than 3 dice. However, it's not all that hard to do the convolution - or even complete enumeration - by hand for small numbers of dice to get exact answers. e.g. here's me doing 3d8 in R:
o2d8=c(outer(1:8,1:8,"+"))
o3d8=c(outer(o2d8,1:8,"+"))
table(o3d8)
o3d8
3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24
1 3 6 10 15 21 28 36 42 46 48 48 46 42 36 28 21 15 10 6 3 1
The table at the end shows the number of ways (out of $8^3=256$) of getting each result on 3d8. You convert them to probabilities by dividing by $8^3$:
> print(table(o3d8)/8^3,d=2)
o3d8
3 4 5 6 7 8 9 10 11 12 13
0.0020 0.0059 0.0117 0.0195 0.0293 0.0410 0.0547 0.0703 0.0820 0.0898 0.0938
14 15 16 17 18 19 20 21 22 23 24
0.0938 0.0898 0.0820 0.0703 0.0547 0.0410 0.0293 0.0195 0.0117 0.0059 0.0020
The algebraic formula becomes relatively complicated past 3 dice, and I wouldn't bother with it, but the probabilities are easy to work out in either R or Excel (or any other tool with the relevant ability to do the kind of calculations you need)
Let's say I want to compute probability of rolling at least 9 on 3d8 from a normal approximation (I suggested more than 3 dice, but let's try it anyway). The exact answer is easily computed by summing probabilities in the above table (it's 0.890625).
> pnorm(8.5,3*(8+1)/2,sqrt(3*(8+1)*(8-1)/12),lower.tail=FALSE)
[1] 0.896144
(The use of 8.5 rather than 9 is because of the continuity correction).
That's not too bad, a relative error of a little over half a percent. But the exact answer takes only a few seconds longer to generate.
Incidentally, a simulation script for rolling 2d8 in R is as simple as:
r2d8=replicate(10000,sum(sample(8,2,TRUE)))
And to display a table of the results as proportions:
table(r2d8)/length(r2d8)
The results of the simulation can be seen as red circles, compared with the exact values (black dots):
--
There's detailed instructions on how to use Excel to do similar calculations to the ones I did in R to compute the exact probabilities here
Best Answer
There are several approaches:
Expand[(x/4+x^2/4+x^3/4+x^4/4)^25]into Wolfram Alpha and press "=" then "show more terms" several times