%matplotlib notebook
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import binom, multinomial
Diablo 2 Drop Calculations¶
© youbetterdont 2019
This article looks at a few of the more challenging probability calculations in Diablo 2's drop system.
n-picks, m-drops¶
Problem Statement¶
In Diablo 2, a monster can drop up to $m$ items, but monsters have up to $n$ picks, where $n$ is possibly greater than $m$. In practice, $m=6$ and the maximum value of $n$ is 7. We would like to know the probability that a monster drops at least $c$ of a desired item, where the probabilty is $p$ on a single pick. The chance of nodrop is $z$ on each pick. The conditional probability of dropping the item is $\frac{p}{1-z}$ on a single pick if nodrop is not selected. The monster will continue to drop items from its $n$ picks until at most $m$ items are selected.
The analysis will make use of binomial random variables. A binomial random variable is denoted as $B(n,p)$, where $n$ is the number of trials and $p$ is the probability of success.
Anaylsis¶
The number of drops before the max limitaiton is $X \sim B(n, 1-z)$. The number of drops is limited to $m$, so we get the actual number of drops $Y = \mathrm{min}(X, m)$.
The number of desired items given $Y$ drops is $U|Y \sim B(Y, \frac{p}{(1-z)})$. We would like to find $\mathrm{Pr}(U \geq c)$ and $\mathrm{E}[U]$. We generate some example data below.
z = 0.3 # probability of nodrop
p = 0.01 # probability of wanted drop
m = 6 # maximum number of drops (drops that aren't nodrop)
n = 7 # maximum number of picks
num_samples = 50000 # number of samples to generate
X = binom(n,1-z)
Xsamples = X.rvs(size=num_samples)
Ysamples = np.minimum(Xsamples, m)
Uarr = [binom(y,p/(1-z)) for y in Ysamples]
Usamples = np.array([u.rvs(size=1)[0] for u in Uarr], dtype=int)
#Usamples = np.array([np.average(u) for u in Usamples_arr], dtype=int)
Plotted below, the histogram of the generated data reveals its distribution.
plt.figure()
plt.hist((Xsamples,Ysamples,Usamples), bins=range(np.maximum(n,m)+2), log=2)
plt.legend(['Xsamples','Ysamples','Usamples'])
plt.ylabel('Occurences')
plt.title('Histogram of generated data\n$n={},z={},p={}$'.format(n,z,p))
plt.xlabel('Number of drops or desired items')
Using scipy.stats.binom.expect, it is easy to get an accurate estimate on the expected value of $Y$.
print('Expected value of Y using binom.expect: {:.6f}'.format(X.expect(lambda x: np.minimum(x, m))))
#print('Intermediate formula: {}'.format(X.expect(ub=m, conditional=True)*X.cdf(m)+m*(1-X.cdf(m))))
print('Average of Ysamples: {:.6f}'.format(np.average(Ysamples)))
Probability¶
By the law of total probability,
\begin{align} \mathrm{Pr}(U \geq c) &= \sum_{k=0}^{m} \mathrm{Pr}(U \geq c| Y=k) \cdot \mathrm{Pr}(Y=k), \end{align}
where $U | Y=k \sim B(k, \frac{p}{1-z})$, that is, the random variable $U$ given $Y=k$ is itself a binomial random variable. The probability of $Y=k$ is
\begin{align} \mathrm{Pr}(Y=k)&= \begin{cases} \mathrm{Pr}(X=k) & k < m \\ \mathrm{Pr}(X \geq m) & k = m \\ 0 & k > m \end{cases}. \end{align}
The probability of $U \geq c$ is then
\begin{align} \mathrm{Pr}(U \geq c) &= \mathrm{Pr}(U \geq c | Y=m) \cdot \mathrm{Pr}(X \geq m) + \sum_{k=0}^{m-1} \mathrm{Pr}(U \geq c| Y=k) \cdot \mathrm{Pr}(X=k) \\ &= \mathrm{Pr}(U \geq c | Y=m) \cdot \mathrm{Pr}(X > m) + \sum_{k=0}^{m} \mathrm{Pr}(U \geq c| Y=k) \cdot \mathrm{Pr}(X=k). \end{align}
The above form of the solution can be evaluated easily using the CDFs of standard binomial random variables $U|Y=k$ and $X$. We can check this against the empirical results from the generated data.
for c in range(1, m+1):
EPrUsamples_geq_c = np.sum(Usamples>=c)/len(Usamples)
print('Empirical probability of U>={}: {:.9f}'.format(c, EPrUsamples_geq_c))
PrUsamples_geq_c = binom.sf(c-1, m, p/(1-z))*X.sf(m)+np.sum([binom.sf(c-1, k, p/(1-z))*X.pmf(k) for k in range(m+1)])
print('Pr(U>={})={:.9f}'.format(c,PrUsamples_geq_c))
The results from the formulas closely match the empirical results, within the margin of error we would expect given the sample size.
As a sanity check, we can consider the case where $n \leq m$.
If $n \leq m$, then $X = Y$ and $U|X \sim B(X, \frac{p}{1-z})$. It follows that $U \sim B(n,p)$, which is simply the number of desired items in $n$ picks. This is the distribution we would expect without the complication of the $m$ pick limit. We can expand the first form to show that that this result is equivalent to the general result above when $n \leq m$.
\begin{align} \mathrm{if}~n \leq m\mathrm{:}\quad\quad&\\ \mathrm{Pr}(U \geq c) &= \sum_{k=0}^{n} \mathrm{Pr}(U \geq c| X=k) \cdot \mathrm{Pr}(X=k) \\ &= \sum_{k=0}^{m} \mathrm{Pr}(U \geq c| Y=k) \cdot \mathrm{Pr}(Y=k) \end{align}
Expectation¶
By the law of total expectation, the expected number of wanted drops is
\begin{align} \mathrm{E}[U] &= \sum_{k=0}^{m} \mathrm{E}[U | Y=k] \cdot \mathrm{Pr}(Y=k) \\ &= \mathrm{E}[U | Y=m] \cdot \mathrm{Pr}(X \geq m) + \sum_{k=0}^{m-1} \mathrm{E}[U | Y=k] \cdot \mathrm{Pr}(X=k) \\ &= \mathrm{E}[U | Y=m] \cdot \mathrm{Pr}(X > m) + \sum_{k=0}^{m} \mathrm{E}[U | Y=k] \cdot \mathrm{Pr}(X=k) \\ &= \frac{p}{1-z} \left (m \cdot \mathrm{Pr}(X > m) + \sum_{k=0}^{m} k \cdot \mathrm{Pr}(X=k) \right ) \end{align}
If $n = m+1$, then
\begin{align} \mathrm{E}[U] &= \frac{p}{1-z} \left (m \cdot \mathrm{Pr}(X = m + 1) + \sum_{k=0}^{m} k \cdot \mathrm{Pr}(X=k) \right ) \\ &= \frac{p}{1-z} \left (- \mathrm{Pr}(X = m + 1) + \sum_{k=0}^{m+1} k \cdot \mathrm{Pr}(X=k) \right ) \\ &= \frac{p}{1-z} \left (\mathrm{E}[X] - \mathrm{Pr}(X=m+1) \right ) \\ &= p \left (n - \frac{\mathrm{Pr}(X=m+1)}{1-z} \right ) \\ &= p \left (m+1 - (1-z)^m \right ) \end{align}
We can easily check the $\mathrm{E}[U]$ result by comparing to the average of the generated Usamples. The resulting relative error is very small and within the margin of error we would expected given the sample size.
Usamples_mean = np.average(Usamples)
print('Average of Usamples: {:.6f}'.format(Usamples_mean))
EU_formula = p*(m+1-(1-z)**m)
print('Formula for expected value of U with n=m+1: {:.6f} (exact solution)'.format(EU_formula))
print('Relative error: {:.6f} (due to variance of E[U])'.format((EU_formula-Usamples_mean)/Usamples_mean))
#print(p/(1-z)*(m*X.sf(m-1)+np.sum([k*X.pmf(k) for k in range(m)])))
#print(binom.expect(args=(m,p/(1-z)))*X.sf(m-1)+np.sum([binom.expect(args=(k,p/(1-z)))*X.pmf(k) for k in range(m)]))
Negative Picks¶
Problem Statement¶
If the picks column is negative, then the selected TCs are not random. The game will drop the specified TCs in order and quantity listed. For example, the TC Countess has Picks=-2 and TCs ordered as Countess Item with Prob=1 and Countess Rune with Prob=1. This means Countess will first drop from Countess Item once, then from Countess Rune once.
This means that the TCs with Picks<0 should not be included in the probability graph. Rather, the monsters that drop from such TCs should reference the derivative TCs with Picks>0 directly. The process is equivalent to killing two or more monsters, where each monster drops from a TC with Picks>0, with the exception that if more than $m$ drops are generated, we can only take the first $m$. The number of monsters is $q$.
The first monster has probability $p_1$ to drop the item and $z_1$ of NoDrop, the second $p_2$ and $z_2$, and so on. The first monster has $n_1$ picks, the second $n_2$ picks, and so on.
The number of drops from the first monster before the $m$ drop limitation is $X_1 \sim B(n_1, 1-z_1)$. Applying the limitation, the number of drops is then $Y_1 = \mathrm{min}(X_1, m)$. For the second monster, $X_2 \sim B(n_2, 1-z_2)$ and $Y_2 = \mathrm{min}(X_2, m-Y_1)$. For the $q$th monster, $X_q \sim B(n_q, 1-z_q)$ and $Y_q = \mathrm{min}(X_q, m-\sum_{k=1}^{q-1} Y_{k})$. The number of wanted items from each monster given $Y_i$ is $U_i|Y_i \sim B(Y_i, \frac{p_i}{1-z_i})$ where $i \in \{1,\dotsc,q\}$. The total number of wanted items is then $U = \sum_{k=1}^{q} U_k$.
We would like to find both $\mathrm{Pr}(U \geq c)$ and $\mathrm{E}[U]$.
Analysis for $q=2$¶
The probabilty and expectation of random variables $X_1$, $Y_1$, and $U_1$ has already been completed in the n-drops, m-picks analysis.
Probability¶
We will focus on $X_2$, $Y_2$, and $U_2$. The number of desired items given $Y_2$ is $U_2|Y_2 \sim B(Y_2, \frac{p_2}{1-z_2})$. By the law of total probability,
\begin{align} \mathrm{Pr}(U_2 \geq c) = \sum_{k=0}^{m} \mathrm{Pr}(U_2 \geq c | Y_2 = k) \cdot \mathrm{Pr}(Y_2=k), \end{align}
where $U_2|Y_2=k \sim B(k,\frac{p_2}{1-z_2})$. The number of drops is $Y_2 = \mathrm{min}(X_2, m-Y_1)$. The probability that $Y_2=k$ is
\begin{align} \mathrm{Pr}(Y_2=k)&= \begin{cases} \mathrm{Pr}(X_2=k)\cdot\mathrm{Pr}(Y_1 < m-k) + \mathrm{Pr}(Y_1=m-k)\cdot\mathrm{Pr}(X_2 > k) + \mathrm{Pr}(X_2=k) \cdot \mathrm{Pr}(Y_1=m-k) & k < m \\ \mathrm{Pr}(X_2 \geq m)\cdot\mathrm{Pr}(Y_1=0) & k = m \\ 0 & k > m \end{cases}. \end{align}
At this point, it probably makes sense to assume $n_i \leq m$. There is no situation in the game where a TC with negative picks leads to a TC with $n_i > m$. If $n_i \leq m$, then $Y_1 = X_1$ and
\begin{align} \mathrm{Pr}(Y_2=k)&= \mathrm{Pr}(X_2=k)\cdot\mathrm{Pr}(X_1 < m-k) + \mathrm{Pr}(X_1=m-k)\cdot\mathrm{Pr}(X_2 > k) + \mathrm{Pr}(X_2=k) \cdot \mathrm{Pr}(X_1=m-k) . \end{align}
Substituting above, we now have an expression for $\mathrm{Pr}(U_2 \geq c)$ that can be evaluated using the PMFs and CDFs of standard binomial random variables. Because of the $n_1 \leq m$ assumption, we have that $U_1 \sim B(n_1, p_1)$. The probability $U_1 \geq c$ is then also easily found using the CDF. The probability of the $U \geq c$ is found as
\begin{align} \mathrm{Pr}(U \geq c) &= \mathrm{Pr}(U_1 + U_2 \geq c). \end{align}
The PMF of the above may be very difficult to compute analytically. At this point, let's restrict $c=1$, then
\begin{align} \mathrm{Pr}(U \geq 1) &= \mathrm{Pr}(U_1 + U_2 \geq 1) \\ &= 1 - \mathrm{Pr}(U_1 = 0 \cap U_2 = 0) \\ &= 1 - \mathrm{Pr}(U_1 = 0 | U_2 = 0) \cdot \mathrm{Pr}( U_2 = 0), \end{align}
where
\begin{align} \mathrm{Pr}(U_2 = 0) = \sum_{k=0}^{m} \mathrm{Pr}(U_2 = 0 | Y_2 = k) \cdot \mathrm{Pr}(Y_2=k). \end{align}
Now the problem is that $U_1$ and $U_2$ are not independent. However, the error introduced by assuming independence is small given that the product of the probabilities $p_1p_2$ is small. Either $U_1=0$ most of the time, or $U_2=0$ most of the time, and either condition results in little effect on the distribution of the other. In other words, we are still sampling a large percentage of the total population of $U_1$ when $U_2=0$. We can then write
\begin{align} \mathrm{Pr}(U \geq 1) \approx 1 - \mathrm{Pr}(U_1 = 0) \cdot \mathrm{Pr}( U_2 = 0). \end{align}
Next, we generate some data to validate the assumption.
z1 = 0.3 # probability of nodrop
p1 = 0.01 # probability of wanted drop
n1 = 5 # maximum number of picks
X1 = binom(n1,1-z1)
X1samples = X1.rvs(size=num_samples)
Y1samples = np.minimum(X1samples, m)
U1arr = [binom(y,p1/(1-z1)) for y in Y1samples]
U1samples = np.array([u.rvs(size=1)[0] for u in U1arr], dtype=int)
#U1samples = np.array([np.average(u) for u in U1samples_arr], dtype=int)
n2 = 4
p2 = 0.05
z2 = 0.35
X2 = binom(n2, 1-z2)
X2samples = X2.rvs(size=num_samples)
Y2samples = np.minimum(X2samples, m-Y1samples)
U2arr = [binom(y,p2/(1-z2)) for y in Y2samples]
U2samples = np.array([u.rvs(size=1)[0] for u in U2arr], dtype=int)
#U2samples = np.array([np.average(u) for u in U2samples_arr], dtype=int)
We plot the histograms of the distributions and compare them to the histograms with the $U_i=0$ condition. We also plot $U_1$ and $U_2$ on a scatterplot to see if there is any obvious correlation. Additionally, we compute the Pearson correlation coefficient and relative error introduced by the unconditional probability estimate.
plt.figure()
plt.hist((U1samples, U1samples[U2samples==0]), density=True, bins=range(m+2), log=True)
plt.legend(['U1', 'U1|U2=0'])
plt.title('Distribution of $U1$ and $U1|U2=0$\n$n_1={},z_1={},p_1={}$'.format(n1,z1,p1))
plt.ylabel('Normalized occurences')
plt.xlabel('Wanted item drops from 1st monster')
plt.figure()
plt.hist((U2samples, U2samples[U1samples==0]), density=True, bins=range(m+2), log=True)
plt.legend(['U2', 'U2|U1=0'])
plt.title('Distribution of $U2$ and $U2|U1=0$\n$n_2={},z_2={},p_2={}$'.format(n2,z2,p2))
plt.ylabel('Normalized occurences')
plt.xlabel('Wanted item drops from 2nd monster')
#plt.hist(Usamples[U2samples==0], density=True, bins=range(m+2))
plt.figure()
plt.scatter(U1samples, U2samples, alpha=0.05)
plt.xlabel('U1')
plt.ylabel('U2')
plt.title('Scatterplot of U2 vs U1')
print("Pearson correlation, corr(U1, U2): {:.6f}".format(np.corrcoef(U1samples, U2samples, rowvar=False)[0,1]))
EPU1eq0 = np.sum(U1samples==0)/len(U1samples)
print('Empirical probability of U1samples=0: {:.6f}'.format(EPU1eq0))
EPU1eq0_U2eq0 = np.sum(U1samples[U2samples==0]==0)/len(U1samples[U2samples==0])
print('Empirical probability of U1samples=0 when U2samples=0: {:.6f}'.format(EPU1eq0_U2eq0))
print('Relative error of emprical probabilities: {:.6f}'.format((EPU1eq0_U2eq0-EPU1eq0)/EPU1eq0))
EPU2eq0 = np.sum(U2samples==0)/len(U2samples)
print('Empirical probability of U2samples=0: {:.6f}'.format(EPU2eq0))
EPU2eq0_U1eq0 = np.sum(U2samples[U1samples==0]==0)/len(U2samples[U1samples==0])
print('Empirical probability of U2samples=0 when U1samples=0: {:.6f}'.format(EPU2eq0_U1eq0))
print('Relative error of emprical probabilities: {:.6f}'.format((EPU2eq0_U1eq0-EPU2eq0)/EPU2eq0))
sum_U1_U2 = U1samples + U2samples
EPUsum_geq_1 = np.sum(sum_U1_U2>=1)/len(sum_U1_U2)
print('Empirical probabilty of U>=1: {:.6f}'.format(EPUsum_geq_1))
#plt.hist(Usamples)
From the above, we can see that the impact of the independence assumption is very small.
We can now check the overall probability estimate.
def PrY2_eq_k(k):
"""
Return the probability that Y2=k, based on the formula above.
"""
return X2.pmf(k)*X1.cdf(m-k-1)+X1.pmf(m-k)*X2.sf(k)+X2.pmf(k)*X1.pmf(m-k)
PrU1_eq_0 = binom.pmf(0,n1,p1)
PrU2_eq_0 = np.sum([binom.pmf(0, k, p2/(1-z2))*PrY2_eq_k(k) for k in range(m+1)])
PrUsum_geq_1 = 1-PrU1_eq_0*PrU2_eq_0
print('Pr(U_1=0)={:.6f} (exact solution)'.format(PrU1_eq_0))
print('Pr(U_2=0)={:.6f}, relative error: {:.6f}'.format(PrU2_eq_0, (PrU2_eq_0-EPU2eq0)/EPU2eq0))
print('Pr(U>=1)={:.6f}, relative error: {:.6f}'.format(PrUsum_geq_1, (PrUsum_geq_1-EPUsum_geq_1)/EPUsum_geq_1))
The probability estimates from the derived formulas agree very well with the empirical probabilities from the generated data.
Expectation¶
Again we will assume $n_i \leq m$. By the law of total expectation,
\begin{align} \mathrm{E}[U_2] &= \sum_{k=0}^{m} \mathrm{E}[U_2 | Y_2 = k] \cdot \mathrm{Pr}(Y_2=k) \\ &= \frac{p_2}{1-z_2}\sum_{k=0}^{m} k \mathrm{Pr}(Y_2=k), \end{align}
with $\mathrm{Pr}(Y_2=k)$ given above in the Probability section. The expected value of $U_1$ is
\begin{align} \mathrm{E}[U_1] = n_1p_1. \end{align}
Then the expected value of the sum is
\begin{align} \mathrm{E}[U] = \mathrm{E}[U_1 + U_2] = \mathrm{E}[U_1] + \mathrm{E}[U_2] = n_1p_1 + \frac{p_2}{1-z_2}\sum_{k=0}^{m} k \mathrm{Pr}(Y_2=k). \end{align}
Next, we compare the results of this formula with the empirical results from the generated data.
EEU1samples = np.average(U1samples)
EEU2samples = np.average(U2samples)
EEUsum = np.average(U1samples+U2samples)
print('Empirical expected value of U1: {:.6f}'.format(EEU1samples))
print('Empirical expected value of U2: {:.6f}'.format(EEU2samples))
print('Empirical expected value of U: {:.6f}'.format(EEUsum))
EU1 = n1*p1
EU2 = p2/(1-z2)*np.sum([k*PrY2_eq_k(k) for k in range(m+1)])
EUsum = EU1 + EU2
print('E[U1]: {:.6f} (exact solution)'.format(EU1))
print('E[U2]: {:.6f} (exact solution)'.format(EU2))
print('E[U]: {:.6f} (exact solution)'.format(EUsum))
The results from the formulas match the empirical results very well.
Duriel drops¶
Problem statement¶
The Duriel TC is the only TC in the game that randomly drops into a TC that itself has more than 1 pick. Other TCs, such as Countess, drop into secondary TCs with multiple picks, but they will always drop a set number of picks from each. This is also true with Super TCs. They drop potions from the Cpot TCs, which have 2 picks each, but they drop a fixed amount of times from Cpot, and their total drops never exceed 6. The point is that if a TC drops a fixed number of times $\alpha$ into a secondary TC with $\beta$ picks, then it's no different than the primary TC dropping $\alpha \beta$ of those picks directly.
With Duriel, the number of times we drop from Duriel - Base is random, and this TC has 7 picks itself. We have no choice but to first roll the dice to see if Duriel - Base is selected, then evaluate the 7 picks if it was. This procedure must repeat another 4 times, or until the maximum of 6 drops.
We assume that the desired item can only come from Duriel - Base, as the only other item Duriel drops is a scroll of town portal (tsc). The probability of getting a Duriel - Base roll is $p_b$. The probability of no-drop from a single Duriel - Base drop is $z$. The probability of getting the desired item in a single Duriel - Base drop is $p$.
Analysis¶
We can use a Markov chain model to calculate the probability of terminating in any of a number of finite states. The system state must be two dimensional because it's not enough to track only the total number of drops (which would include the tsc drops). It's also not enough to use only the number of Duriel - Base drops because we do not know how many drops are remaining. Therefore, the system state is chosen as the ordered pair $(T_n, Y_n)$, where $T_n$ is the number of tsc drops at step $n$ and $Y_n$ is the number of Duriel - Base drops at step $n$.
There are 28 possible states in the general $m=6$ case. They are enumerated below. We choose to include the $(5,1)$ and $(6,0)$ states for completeness, but in reality, only up to five Duriel Drops will occur, and these states are inaccessible.
| T | Y | count |
|---|---|---|
| 0 | {0, 1, 2, 3, 4, 5, 6} | 7 |
| 1 | {0, 1, 2, 3, 4, 5} | 6 |
| 2 | {0, 1, 2, 3, 4} | 5 |
| 3 | {0, 1, 2, 3} | 4 |
| 4 | {0, 1, 2} | 3 |
| 5 | {0, 1} | 2 |
| 6 | 0 | 1 |
What remains is to calculate the probability of advancing from one state to the next in order to populate the transition matrix. The transition matrix is 28 $\times$ 28. It is upper triangular, leaving 406 entries to populate. The upper triangular property is because the transitions are directed. It's possible to advance the state to a greater or equal number of items, but not to go in the reverse direction.
We must order the state vector in a manner to enforce the upper triangular condition. The sum of $T$ and $Y$ must not decrease as we advance though the vector. This means that $T$ and $Y$ are never permitted to both decrease from one element to the next, otherwise the sum would decrease. Equivalently, this means that traversing the elements in reverse order necessitates either $T$ or $Y$ to decrease, which is not a permitted transition, thus has zero probability. Therefore, the matrix is upper triangular.
We order the states the following way.
states = sorted([(x,y) for x in range(7) for y in range(7-x)], key=lambda x: x[0]+x[1])
print("Number of states: {}".format(len(states)))
print("Ordered states: {}".format(states))
There are two types of possible transitions. They are
\begin{align} (t, y) &\rightarrow (t, y+k) \\ (t, y) &\rightarrow (t+1, y), \end{align}
where $k \geq 0$. The first possibility is that Duriel - Base drops zero or more items. The second is the possibility that Duriel drops a tsc directly. Any transitions that do not fit the above description have zero probability. For the first case, we find
\begin{align} &\mathrm{Pr}(T_{n+1} = t \cap Y_{n+1} = y+k | T_{n} = t \cap Y_{n} = y) \\ & \qquad = \mathrm{Pr}(T_{n+1}-T_n = 0 \cap Y_{n+1}-Y_{n} = k | T_{n} + Y_{n} = t + y) \\ & \qquad = \mathrm{Pr}(T_{n+1}-T_n = 0 | T_{n} + Y_{n} = t + y) \cdot \mathrm{Pr}(Y_{n+1}-Y_{n} = k | T_{n+1}-T_n = 0 \cap T_{n} + Y_{n} = t + y). \end{align}
The first term is the probability that Duriel will not drop another tsc given how many total drops have occured. We find
\begin{align} \mathrm{Pr}(T_{n+1}-T_n = 0 | T_{n} + Y_{n} = t + y) &= \begin{cases} p_b & t+y < 6 \\ 1 & t+y \geq 6 \end{cases}, \end{align}
The above states that Duriel cannot drop any more tscs once the drop limitation is reached.
Much earlier, we found that the number of drops $Y$ from a monster with $n$ picks and $m$ maximum drops has the PMF
\begin{align} \mathrm{Pr}(Y=k)&= \begin{cases} \mathrm{Pr}(X=k) & k < m \\ \mathrm{Pr}(X \geq m) & k = m \\ 0 & k > m \end{cases}, \end{align}
where $X \sim B(n, 1-z)$ is the number of drops before the maximum limitation. We can modify this slightly here to say
\begin{align} \mathrm{Pr}(Y_{n+1}-Y_{n} = k | T_{n+1}-T_n = 0 \cap T_{n} + Y_{n} = t + y) &= \begin{cases} \mathrm{Pr}(X=k) & k < 6-t-y \\ \mathrm{Pr}(X \geq 6-t-y) & k = 6-t-y \\ 0 & k > 6-t-y \end{cases}. \end{align}
The probability of the second type of transition is
\begin{align} &\mathrm{Pr}(T_{n+1} = t+1 \cap Y_{n+1} = y | T_{n} + Y_{n} = t + y) \\ & \qquad = \mathrm{Pr}(T_{n+1}-T_n = 1 \cap Y_{n+1}-Y_{n} = 0 | T_{n} + Y_{n} = t + y) \\ & \qquad = \mathrm{Pr}(T_{n+1}-T_n = 1 | T_{n}+Y_{n}=t+y) \cdot \mathrm{Pr}(Y_{n+1}-Y_{n} = 0 | T_{n+1}-T_n = 1 \cap T_{n} + Y_{n} = t + y) \\ & \qquad = \mathrm{Pr}(T_{n+1}-T_n = 1 | T_{n} + Y_{n} = t + y) \\ & \qquad = \begin{cases} (1-p_b) & t+y < 6 \\ 0 & t+y \geq 6 \end{cases} . \end{align}
We now have all the necessary PMFs to fill out the matrix. We do so programmatically below.
pb = 2./3
Td = np.zeros((len(states),len(states)))
for i, row in enumerate(Td):
for j, _ in enumerate(row):
tn, yn = states[i]
tn1, yn1 = states[j]
k = yn1 - yn
if tn==tn1 and k >= 0:
Td[i, j] = (X.sf(6-tn-yn-1) if k==6-tn-yn else 0 if k > 6-tn-yn else X.pmf(k))*(pb if tn+yn < 6 else 1)
elif tn1-tn==1 and k == 0:
Td[i, j] = (1-pb) if tn+yn < 6 else 0
print('The sum of each row of the transition matrix (should be all ones):\n {}'.format(np.sum(Td, axis=1)))
We can now find the probability of ending in any state at step $n$ by multiplying the initial state row vector by the transition matrix, that is,
\begin{align} \mathbf{x}_n = \mathbf{x}_o \mathbf{P}^n. \end{align}
The initial state $\mathbf{x}_o$ has a 1 in the first column and is zero everywhere else, meaning we start in state $(0,0)$. The probability of ending in each state is found as above. The PMF of the number of Duriel - Base drops can be computed by summing the probabilities of the states with identical Duriel - Base drops. Below, we perform this calculation and also generate a large sample of data for comparison. The results are presented on a combined histogram.
x_o = np.zeros((1,len(states)))
x_o[0,0] = 1
x5 = np.dot(x_o, np.linalg.matrix_power(Td, 5))
y5 = [np.sum(x5[0,np.array(states)[:,1]==i]) for i in range(7)]
num_d = 50000
Td_samples = np.zeros(num_d).astype('int')
Yd_samples = np.zeros(num_d).astype('int')
for i in range(len(Yd_samples)):
for j in range(5):
t = binom.rvs(1, 1-pb)
Td_samples[i] += t
if t == 0:
Yd_samples[i] += np.minimum(X.rvs(), m-Td_samples[i]-Yd_samples[i])
if Td_samples[i] + Yd_samples[i] >= m:
if Td_samples[i] + Yd_samples[i] > m: print('This should not happen.')
break
plt.figure()
plt.bar(np.arange(7), y5, color='C0', width=0.5)
plt.hist(Yd_samples, bins=np.arange(8), density=True, rwidth=0.5, color='C1')
plt.legend(['calculated', 'simulated (N={})'.format(num_d)])
plt.title('PMF of "Duriel - Base" drops\nwith NoDrop={:.3f}'.format(z))
plt.xlabel('Number of "Duriel - Base" drops')
plt.ylabel('Probability or normalized occurences')
From the above, we can see that the results calculated through the Markov chain model are almost exactly matching the simulated results.
The complete transition matrix is given below.
import sys
print('Transition matrix for z={:.3f}'.format(z))
np.savetxt(sys.stdout, Td, fmt='%.3e', delimiter=',')