Myujth

The problem is basically as follows

There are four cards with points [1,2,3,4], every time I draw one card randomly and then put it back to the deck. As I keep drawing I record the sum of the points on the cards, I stop drawing if the sum is greater than or equal to five.

If the final points is equal to five, I win, otherwise I lose. What's the probability of me winning?

My attempt to solve this problem is, let $p_n$ be the probability of the sum EVER equal to $n$ after drawing enough number of times, specifically
$$p_1=frac14 text (draw 1 once)$$
$$p_2=frac14 + frac14p_1 text (draw 2 once + draw 1 twice) $$
$$p_3=frac14 + frac14p_1 + frac14p_2 $$
$$p_4=frac14 + frac14p_1 + frac14p_2 + frac14p_3 $$
then $p_5=frac14p_1 + frac14p_2 + frac14p_3 + frac14p_4 = 369/1024$ should be the probability of winning at five.

However I was told that the answer should be around 0.33 so I'm confused, is my answer correct?

Update
It turns out 0.33 (1/3) is the answer for sampling without replacement.

Whenever there are conflicting answers from combinatorial methods, I
like to simulate the result to see what happens.

My method of simulation (in R) is to sample five cards with replacement from among 1, 2, 3, 4. Then to take cumulative sums of the results to see if the total 5 is
present at any step. If so, that's a win. (Any cards drawn after a total of 5
are ignored.)

After a million iterations the vector w has a million
TRUEs and FALSEs. The mean of this logical vector is the proportion
of its TRUEs. With a million iterations this proportion should approximate
the probability of winning to at least two or three places.

It seems that
simulation results match $369/1024 = 0.3603516,$ the answer proposed in your
Question and computed in the Comment
of @a_statistician.

set.seed(1112); m = 10^6; deck = 1:4
w = replicate( m, 5 %in% cumsum(sample(deck,5, rep=T)) )
mean(w); 369/1024; 15/49
[1] 0.360286 # aprx P(Win) = 0.3604
[1] 0.3603516 # exact P(Win) = 369/1024
[1] 0.3061224 # 15/49

I think another way to approach this problem is to model the game as a Markov chain with the following transition probability matrix, $M$:

$$
beginmatrix
& textFrom state \
textTo state &
beginarrayccccccc
& 0 & 1 & 2 & 3 & 4 & textwin & textlose \
hline
0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \
1 & 1/4 & 0 & 0 & 0 & 0 & 0 & 0 \
2 & 1/4 & 1/4 & 0 & 0 & 0 & 0 & 0 \
3 & 1/4 & 1/4 & 1/4 & 0 & 0 & 0 & 0 \
4 & 1/4 & 1/4 & 1/4 & 1/4 & 0 & 0 & 0 \
textwin (5) & 0 & 1/4 & 1/4 & 1/4 & 1/4 & 1 & 0 \
textlose (5+) & 0 & 0 & 1/4 & 2/4 & 3/4 & 0 & 1
endarray
endmatrix
$$

Note that this is a reducible Markov chain: once we reach the win or lose state, we will be stuck there forever.

Our starting state is:
$$
v = beginbmatrix
1 & 0 & 0 & 0 & 0 & 0 & 0
endbmatrix^top
$$

Since we know that the game will end after five draws at most, we can calculate the probability distribution over the states after five draws by left multiplying $M$ with the initial state vector $v$ five times:
$$
M times M times M times M times M times v
$$

Implementation in R:

library(expm)

M <- matrix(nrow = 7,
 c( 0, 0, 0, 0, 0, 0, 0, 
 .25, 0, 0, 0, 0, 0, 0, 
 .25, .25, 0, 0, 0, 0, 0, 
 .25, .25, .25, 0, 0, 0, 0, 
 .25, .25, .25, .25, 0, 0, 0, 
 0, .25, .25, .25, .25, 1, 0, 
 0, 0, .25, .5, .75, 0, 1),
 byrow = T, 
 dimnames = list(c("0", "1", "2", "3", "4", "win", "lose"), 
 c("0", "1", "2", "3", "4", "win", "lose")))

M %^% 5 %*% c(1, 0, 0, 0, 0, 0, 0)

Output:

 [,1]
0 0.0000000
1 0.0000000
2 0.0000000
3 0.0000000
4 0.0000000
win 0.3603516
lose 0.6396484

The probability for the win state is consistent with $frac3691024$.

P.S. I do wonder though if there is a way to derive the distributions over the final states analytically without carrying out the matrix multiplications.