Can a Markov chain be used to compute probabilities of sets of draws?

Question

Problem statement

Consider the following scenario:

We have an extremely large bag of glass marbles and metallic balls. We do not know how many items are in the bag, but we do know the percentage of the total for each type of metal ball:

Color	Percentage
Gold	1%
Silver	2%
Copper	5%

We are drawing a single ball at a time, replacing it, and mixing up the bag. It is not possible to distinguish what we are drawing by feel, so what we draw is random.

Our goal is to draw (see) each metal at least once, and we want to examine the probability of having completed the set over the first thousand draws. Just to give us some specific values to compute, let's say we'll look at the probabilities at 50, 100, 250, 500, and 1000 draws.

Solution, according to my understanding

I believe we can approach this problem using a Markov chain. We can model the problem by having each state be a subset of metals that we have encountered, and the transition matrix will be the probability of drawing any unseen metal. While there is no interesting long term steady state (The probability of each state will asymptotically approach 0 except for the completed set which will approach 1.), the chain will at the very least help organize the computations in a much simpler way than trying to perform them ad hoc.

To get our states, we can take the power set of the metals, and we can use some abbreviations to shorten things:

Seen	Abbreviation
{}	$ \emptyset $
{Gold}	G
{Silver}	S
{Copper}	C
{Gold,Silver}	GS
{Gold,Copper}	GC
{Silver,Copper}	SC
{Gold,Silver,Copper}	GSC

Before we lay out our transition matrix, let's compute the probability of drawing a glass marble. This value will be useful for constructing the transition matrix:

$$ 1 - (0.01 + 0.02 + 0.05) = 0.92 $$

For our transition matrix, we will use the $ X_n = P X_{n+1} $ convention, meaning each column will represent a current state and each row will represent the next state. So our transition matrix $ P $, with some labels to make interpreting it easier, looks like this:

$$ \newcommand\pad[1]{\rlap{#1}\phantom{0000}}\begin{matrix} From: & \begin{array}{*8c} \pad{\emptyset} & \pad{\text{G}} & \pad{\text{S}} & \pad{\text{C}} & \pad{\text{GS}} & \pad{\text{GC}} & \pad{\text{SC}} & \pad{\text{GSC}} \end{array} & \\ & & To: \\ & \begin{bmatrix} \pad{0.92} & \pad{0} & \pad{0} & \pad{0} & \pad{0} & \pad{0} & \pad{0} & \pad{0} \\ \pad{0.01} & \pad{0.93} & \pad{0} & \pad{0} & \pad{0} & \pad{0} & \pad{0} & \pad{0} \\ \pad{0.02} & \pad{0} & \pad{0.94} & \pad{0} & \pad{0} & \pad{0} & \pad{0} & \pad{0} \\ \pad{0.05} & \pad{0} & \pad{0} & \pad{0.97} & \pad{0} & \pad{0} & \pad{0} & \pad{0} \\ \pad{0} & \pad{0.02} & \pad{0.01} & \pad{0} & \pad{0.95} & \pad{0} & \pad{0} & \pad{0} \\ \pad{0} & \pad{0.05} & \pad{0} & \pad{0.01} & \pad{0} & \pad{0.98} & \pad{0} & \pad{0} \\ \pad{0} & \pad{0} & \pad{0.05} & \pad{0.02} & \pad{0} & \pad{0} & \pad{0.99} & \pad{0} \\ \pad{0} & \pad{0} & \pad{0} & \pad{0} & \pad{0.05} & \pad{0.02} & \pad{0.01} & \pad{1} \\ \end{bmatrix} & \begin{array}{*8wc100} \pad{\emptyset} \\ \pad{\text{G}} \\ \pad{\text{S}} \\ \pad{\text{C}} \\ \pad{\text{GS}} \\ \pad{\text{GC}} \\ \pad{\text{SC}} \\ \pad{\text{GSC}} \end{array} \end{matrix} \\ $$

The initial state $ X_0 $ is simple: we start with the empty set 100% of the time. So it's just an 8 element matrix where the first element is 1 and the remaining elements are 0.

With these two matrices in hand, computing the probability of having drawn all metals is just a matter of matrix multiplication and extracting the value of interest.

$$ \begin{align} X_{50} &= P^{50} X_{0} \\ X_{100} &= P^{100 - 50} X_{50} = P^{50} X_{50} \\ X_{250} &= P^{250 - 100} X_{100} = P^{150} X_{100} \\ X_{500} &= P^{500 - 250} X_{250} = P^{250} X_{250} \\ X_{1000} &= P^{1000 - 500} X_{500} = P^{500} X_{500} \\ \end{align} $$

And we find that the probabilities for completion are:

Draws	Probability
50	0.22836397 = 22.8364%
100	0.54550092 = 54.5501%
250	0.91302709 = 91.3027%
500	0.99338874 = 99.3389%
1000	0.99995683 = 99.9957%

Is this analysis correct?

Motivation

If anyone is wondering, this is not a homework problem. I'm planning to apply the technique to analyzing gacha game probabilities, and I'd like to be certain that the approach I'm using actually works.

If this approach is correct, I am interested in generalizing to more complex situations. For example, I want to examine the effect of having two bags to choose from (where you decide which bag to draw from based on the seen set you already have). I may also be interested in the probabilities of the various subsets, as well. But before I get into those, I wanted to be certain that the approach is valid.

Additional properties?

Given the context, are there any well known properties that a Markov chain exhibits that might be of interest in this sort of problem space? I can go read about them on my own, but it would help if I know what to look for.

Extra details

I chose to distinguish the items of interest as metals rather than colors because not having to distinguish between colors of interest and colors not of interest made writing the problem simpler.

I've omitted the actual computations of the final probabilities because I used Python and numpy to perform them rather than do them by hand. If you'd like to perform computations yourself, here's the transition matrix in code form (with the assumption you've done the conventional import numpy as np):

np.array([[0.92, 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  ],
          [0.01, 0.93, 0.  , 0.  , 0.  , 0.  , 0.  , 0.  ],
          [0.02, 0.  , 0.94, 0.  , 0.  , 0.  , 0.  , 0.  ],
          [0.05, 0.  , 0.  , 0.97, 0.  , 0.  , 0.  , 0.  ],
          [0.  , 0.02, 0.01, 0.  , 0.95, 0.  , 0.  , 0.  ],
          [0.  , 0.05, 0.  , 0.01, 0.  , 0.98, 0.  , 0.  ],
          [0.  , 0.  , 0.05, 0.02, 0.  , 0.  , 0.99, 0.  ],
          [0.  , 0.  , 0.  , 0.  , 0.05, 0.02, 0.01, 1.  ]])

Answer 1

Since probabilities remain constant due to replacement, why not simply apply inclusion-exclusion, eg for $100$ trials, Wolfram yields an answer of

$\small0.5455009161013792606431705561885349517959643720020683437164727651...$

Answer 2

The approach is fine. Also, just computing the n-th power of the transition matrix is also correct (as far as I have seen). The only thing is that your transition matrix is transposed to what is standard in the literature. I am using the standard notation below, as it is easier to write down.

However, you can extract a lot more information from the transition matrix. Using your notation, we can split the states into transient states ($\emptyset$, {G}, {C}, {S}, {G,S}, {G,C} and {S,C}) and absorbing ones ({G,S,C}). We collect transition probabilities between the transient states in a matrix $S$ and the transition probabilities to the absorbing state with a column vector $R$ (last row of your matrix without 1 and transposed). The transition matrix then reads $$P=\begin{pmatrix}S&R\\0&1\end{pmatrix}.$$

If we call the random variable $T$ the time until all metal marbles have been seen for the first time, we can explicitly give all interesting quantities. For example $$\mathbb E(T)=\nu\,(Id-S)^{-1}R\approx 116.$$ Note that $\nu$ is a column vector that corresponds to the initial condition, which might be for example $\nu=(1,0,0,0,0,0,0,0)$, the case where you start with the empty set. "Id" is the identity matrix.

If you want to extract more information, the maximum principle might be helpful, or cost functionals. It really depends on what you want to know. All of this should be computable exactly, since the matrix is so small.

I think it would be helpful if you could add a list of quantities that you need to compute to your original post. It makes it easier to decide which approach is more useful for your specific problem.