Probability distribution of a step by step process I've come across a "sort of" simple probability problem but I don't know where to start to model it mathematically. The process is the following:

*

*Flip a coin.

*If I get tails I roll a die, write down its value and go back to step 1; otherwise the process goes to step 3.

*Add up all the values obtained from rolling the dice and call it $X$.  (If I don't get tails the first time I flip the coin then $X = 0$.)

Is it possible to derive a PDF for the value $X$?
What would be the name of this process and how can I model it? I've simulated the process in Python and can't make sense of the shape of the resulting histogram of $X$.
 A: The process is known as a compound process.
First, let's start by modeling the number of coin flips that come up heads.   Since we are flipping coins until we get a tail, the distribution of the number of heads before we see the first tail is what is of interest, and this is a geometric distribution.
The number distribution on the die is simple; it's uniformly distributed over the integers $1, \dots, 6$.
What we are interested in finding is the distribution of:
$$y = \sum_{i=1}^Nx_i\,\,\text{where  }N \sim \mathrm{Geometric}(p=1/2),\,x_i \sim \mathrm{Uniform}(1,2,3,4,5,6)$$
Unfortunately, there is no nice formula for this; even the distribution when $N=2\,\,(x_1 + x_2)$ is enumerated.  However, it can be approximated with a little effort.
Our approach is to calculate the distribution of $\sum x_i$ for each $N$, then add them together, weighting by the probability of $N$ calculated from the geometric distribution.  Since I'm not so familiar with Python, I've written the code in R:
foo <- function(x, p) {
  x_out <- rep(0, length(x) + 6)
  p_out <- rep(0, length(x) + 6)
  
  for (i in 1:6) {                         # Iterating over the next dice roll
    for (j in 1:length(x)) {               # Iterating over the current sum
      x_out[i+j] <- i + x[j]               # Adding the dice roll to the sum
      p_out[i+j] <- p_out[i+j] + p[j] / 6  # Calculating the probability
    }
  }
  
  x_out <- x_out[p_out > 0]
  p_out <- p_out[p_out > 0]
  list(x = x_out, p = p_out)
}


p_y <- rep(0, 100)
tmp <- list(x=1:6, p=rep(1/6,6))
p_y[tmp$x] <- 1/24 # 1/6 * 1/4

for (N in 2:16) {
  tmp <- foo(tmp$x, tmp$p)
  p_y[tmp$x] <- p_y[tmp$x] + tmp$p / (2^(N+1))
}

# Add in the 50% probability of a zero
p_y <- c(0.5, p_y)

We have to cut off our iteration at some point, so I've arbitrarily done so at $N = 16$.  This gives pretty good results:
> sum(p_y)  # ideally equal to 1
[1] 0.9999924

We are missing very little probability mass.
And a plot:
plot(p_y ~ c(0:100),
     xlab = "Sum of dice rolls",
     ylab = "Probability",
     type = "b",
     pch = 16)


