Probability problem, roll a n-sided dice for m times, the probablity to have at least d difference between the largest number and the smallest number. 
You're playing a game, you may choose a fair n-sided die.
  n = 4,6,8,12,20
  And pay 1 buck for rolling the die twices.
  You are allowed to roll the die more than twice, up to n times, each extra roll charge 1 buck.
  The total times of rolling the die is denoted as m.
  m = 2,3,4,...,n.The rule is simple, shout out a number d before rolling the die.
  Within the m rolls, record down the largest number max and the smallest number min.
  If max - min >= d then you win d bucks, otherwise you lose all your wager. 



*

*What is the winning probability as a function of n,m, and d?

*And what is the best strategy to choose n,m, and d?


My thought:
Start from the default m = 2 case.
Denote the probablity to get number i from the first trail as $P(die=i)$, where $i \in [1,n]$, $P(die=i)$ satisfy the continuous uniform distribution. $P(die=i)=unifpdf(i,1,n)$.
And then, roll the die for the second time, the result is denoted as j, again $P(die=j)=unifpdf(j,1,n)$.
The probability to win $P_{a}$ can be written as:
\begin{equation}
P_{a} = \sum_{i=1}^{n}\Big(unifpdf(i,1,n) \times \sum_{j=1}^{n}\big(unifpdf(j,1,n) \times I_{A}(i,j)\big)\Big)
\end{equation}
Where $I_A(i,j)$ is an indicator function.
\begin{equation}
  I_A(i,j)=\left\{
  \begin{array}{@{}ll@{}}
    1, & \text{if}\ |i-j| \geq d \\
    0, & \text{otherwise}
  \end{array}\right.
\end{equation} 
Consider m > 2 trails. within $\binom{m}{2}$ combinations if there is one pair satisfies $|i-j|\geq d$ then you win.
\begin{equation}
P_{win} = 1-(1-Pa)^{\binom{m}{2}}
\end{equation}
Where $(1-Pa)$ is the probability that one pair doesn't satisfies $|i-j| \geq d$, to the power of $\binom{m}{2}$ means all combination pairs in m trails are failed, take one mimus means at least one pair successed.
Pwin probablity of winning
Pwin*d expected reword
Pwin*d-(m-1) expected reword deduct wager cost
And the matlab code to obtain the figure above.
close;
clc;
%% constants
N = 20; % a n-sided die
M = N;
D = N;
%% variables
Pwin = NaN(D,M);
%% calculation
for idxd = 1:N % a judgement number n to shout out
    for idxm = 2:N % trail number m up to n 
        Pa = 0; % sum i=0->n
        for i = 1:N % the possible value of first roll
            P1st = unifpdf(i,1,N); %the die follows continous uniform distribution
            Pb = 0; % sum j=0->n
            for j = 1:N % the possible value of the second roll
                if(abs(i-j)>=idxd)
                    P2nd = unifpdf(j,1,N);%the die follows continous uniform distribution
                else
                    P2nd = 0;
                end
                Pb = Pb + P2nd;
            end
            Pa = Pa + P1st*Pb;      
        end
        Pwin(idxd,idxm) = 1-(1-Pa)^(nchoosek(idxm,2));
   end
end

%% random test
poll = 10000; % repeat 100 round of game
dice = randi(N,[1,poll]);

Prand = NaN(D,M);
for idxd = 1:N
    for idxm = 2:M
        k = 0; %k is the round wining
        b = 0; %b is the number of round played
        for i = 1:idxm:length(dice)-(idxm-1)
            [maxHW maxI] = max(dice(i:i+idxm-1));
            [minHW minI] = min(dice(i:i+idxm-1));
            b = b + 1;
            if(maxHW - minHW >= idxd) 
                k = k + 1;
            end
        end
        Prand(idxd,idxm) = (k/b);
    end
end

%% plot the Pwin
figure(1);
vector_m = 2:M;
for idxd = 1:N
    subplot(ceil(D/4),4,idxd);
    plot(vector_m,Pwin(idxd,2:M),vector_m,Prand(idxd,2:M));
    sttr = '';
    sttr = sprintf('d = %d',idxd);
    title(sttr);
    legend('Eqn.','Rand test');
    ylim([0 1]);
    ylabel('Pwin');
    xlabel('m trails');
end

%%plot the Pwin*d
% if you play many times, average money win 
figure(2);
vector_m = 2:M;
for idxd = 1:N
    subplot(ceil(D/4),4,idxd);
    plot(vector_m,Pwin(idxd,2:M)*idxd,vector_m,Prand(idxd,2:M)*idxd);
    sttr = '';
    sttr = sprintf('d = %d',idxd);
    title(sttr);
    legend('Eqn.','Rand test');
    ylabel('Pwin*d');
    xlabel('m trails');
end

%%plot the (Pwin*d) - (m-1) 
% net wining, money win deduct you spent
figure(3);
vector_m = 2:M;
for idxd = 1:N
    subplot(ceil(D/4),4,idxd);
    plot(vector_m,(Pwin(idxd,2:M)*idxd - (vector_m-1)),vector_m,(Prand(idxd,2:M)*idxd) - (vector_m-1));
    sttr = '';
    sttr = sprintf('d = %d',idxd);
    title(sttr);
    legend('Eqn.','Rand test');
    ylabel('(Pwin*d)-(m-1)');
    xlabel('m trails');
end

There is a considerable disagreement between the equation and the random test.
Could anyone help me to find out why?
Many thanks!
 A: The main problem I think is your formula 
$$P_{win}=1-(1-P_a)^{m \choose 2}.$$
It assumes that the probability to not get the desired distance $\ge d$ between a pair of rolls is independent from any other pair. That's incorrect, and I assume in a big way. If you for example start with three different rolls with values $a < b < c$, then the roll that gave you $b$ is not in any way 'helping' you to reach the distance $d$, so all the other $m-3$ pairs of $b$ with the remaining rolls are not contributing to the probability.
What you need to do is calculate the probability that among m rolls ($r_1,r_2,\ldots,r_m$), the minimum roll is any given value $a$ from $1$ to $n$, and the maximum roll is another given value $z$ from $1$ to $n$:
$$P(a,z):={\rm Prob}(\min(r_1,r_2,\ldots,r_m)=a \land\max(r_1,r_2,\ldots,r_m)=z)$$
Once you have that, your requested probability $P_a$ is simply the sum $P(a,z)$ over all cases where $z-a\ge d$.
$$P_a=\sum_{a=1}^{n-d}\sum_{z=a+d}^n P(a,z) \tag{1} \label{eqstar}$$
$P(a,z)$ is a little tricky to calculate, but there is a similar probability that is easy to calculate:
$$Q(a,z):={\rm Prob}(\min(r_1,r_2,\ldots,r_m)\ge a \land\max(r_1,r_2,\ldots,r_m)\le z) = \left(\frac{z-a+1}{n}\right)^m. \tag{2} \label{eq1}$$
Because another way to look at $Q(a,z)$ is to ask: What is the probability that among $m$ rolls all values are at least $a$ and at most $z$? That question can be answered, because it poses an independent question on each die roll $r_i$: Is $a \le r_i \le z$ or not? The answer must be "yes" for all $m$ rolls, and it is "yes" for a given roll with probability $(\frac{z-a+1}{n})$ (the usual formula for the number of good outcomes divided by the number of all outcomes) and because those events are independet, the probabiltiies can be multiplied.
How do we get from $Q(a,z)$ to $P(a,z)$? 
Note that 
$$\{(r_1,r_2,\ldots,r_m): (\min(r_1,r_2,\ldots,r_m)=a \land\max(r_1,r_2,\ldots,r_m)=z)\} = \{(r_1,r_2,\ldots,r_m): (\min(r_1,r_2,\ldots,r_m)\ge a \land\max(r_1,r_2,\ldots,r_m)=z)\} - \{(r_1,r_2,\ldots,r_m): (\min(r_1,r_2,\ldots,r_m)\ge a+1 \land\max(r_1,r_2,\ldots,r_m)=z)\}$$
because for the minimum to be equal to $a$ it is equivalent to ask for the minium to be at least $a$ but not to be at least $a+1$. Also, the set we are subracting is actually a subset of the set we are subtracting from, so we get:
$$P(a,z)= {\rm Prob}\{(r_1,r_2,\ldots,r_m): (\min(r_1,r_2,\ldots,r_m)=a \land\max(r_1,r_2,\ldots,r_m)=z)\} = {\rm Prob}\{(r_1,r_2,\ldots,r_m): (\min(r_1,r_2,\ldots,r_m)\ge a \land\max(r_1,r_2,\ldots,r_m)=z)\} - {\rm Prob}\{(r_1,r_2,\ldots,r_m): (\min(r_1,r_2,\ldots,r_m)\ge a+1 \land\max(r_1,r_2,\ldots,r_m)=z)\}$$
Our two operands are of the same type, we can do the same trick to reduce each probability to something expressable by some $Q(r,s)$:
$$\{(r_1,r_2,\ldots,r_m): (\min(r_1,r_2,\ldots,r_m)\ge r \land\max(r_1,r_2,\ldots,r_m)=z)\} = \{(r_1,r_2,\ldots,r_m): (\min(r_1,r_2,\ldots,r_m)\ge r \land\max(r_1,r_2,\ldots,r_m) \le z)\} -  \{(r_1,r_2,\ldots,r_m): (\min(r_1,r_2,\ldots,r_m)\ge r \land\max(r_1,r_2,\ldots,r_m) \le z-1)\}$$
and we get
$${\rm Prob}\{(r_1,r_2,\ldots,r_m): (\min(r_1,r_2,\ldots,r_m)\ge r \land\max(r_1,r_2,\ldots,r_m)=z)\} = {\rm Prob}\{(r_1,r_2,\ldots,r_m): (\min(r_1,r_2,\ldots,r_m)\ge r \land\max(r_1,r_2,\ldots,r_m) \le z)\} -  {\rm Prob}\{(r_1,r_2,\ldots,r_m): (\min(r_1,r_2,\ldots,r_m)\ge r \land\max(r_1,r_2,\ldots,r_m) \le z-1)\} = Q(r,z) - Q(r,z-1)$$
and finally for $P(a,z)$, by subsituting this in the above forumla:
$$P(a,z)=(Q(a,z)-Q(a,z-1)) - (Q(a+1,z)- Q(a+1,z-1)). \tag{3} \label{eq2}$$
So \eqref{eqstar}, \eqref{eq1} and \eqref{eq2} should allow you to compute via code the theoretical $P_a$. Maybe the difference structure of \eqref{eq2} allows some more simplification when substituting into \eqref{eqstar}, I didn't check.
Note that when deriving equation \eqref{eq1}, I assumed $1\le a \le z \le n$. But adding $1$ to $a$ and subtracting $1$ from $z$ makes that assumption invalid. \eqref{eq1} still gives the correct result (zero) for $z=a-1$, but provides an obviously incorrect negative value when $z < a-1$. That can happen when you start with $a=z=k$ and have to evaluate $Q(a+1,z-1)$ in \eqref{eq2}. 
So when calculating $Q(r,s)$, always check if $1 \le r$ and $s \le n$ first. If $r$ is smaller than $1$, set it to $1$ and if $s$ is bigger than $n$, set it to $n$. Then check if $r \le s$ and substitute $0$ otherwise instead of using \eqref{eq1}.
