Bayesian Inference: Modeling checkout times at a store

I am currently learning how to use Bayesian inference. I have been making up problems (by defining some population parameters) and then trying to infer those values from samples. I recently made up a problem that I cannot figure out how to solve.

I would like to model the time that it takes to checkout at a grocery store. In this model, there are two main parts. The first is, the time it takes the cashier to ring up and bag the items. The second is, the time it takes to pay the cashier.

I have chosen to model the time it takes the cashier to scan and bag $$k$$ items (assume for simplicity that the same number of items are always purchased) as a gamma distribution:

$$S \sim \textrm{Gamma}(k, \theta_1)$$

The $$\theta_1$$ parameter is the average time it takes to scan and bag one item. So this is modeling the sum of $$k$$ exponential random variables (with rate parameter $$1/\theta_1$$).

I have chosen to model the time it takes to pay the cashier as an exponential variable:

$$P \sim \textrm{Exponential}(\theta_2)$$

The total time it takes to checkout is the sum of the time it takes to scan and bag plus the time it takes to pay:

$$W_i \sim S + P$$

If I directly observe the random variables $$s_i$$ and $$p_i$$ (and I know the number of items k), then the inference problem for $$\theta_1$$ and $$\theta_2$$ is easy. The likelihood function for $$s_i$$ is the gamma distribution and the likelihood function for $$p_i$$ is the exponential distribution. Standard Monte Carlo sampling techniques can be used to obtain their posterior distributions.

However, if all that I know are the total waiting times $$w_i$$ (and $$k$$), then I do not know what the likelihood function should be. From the research I have done, it appears that in general sums of random variables do not have closed form PDFs.

How can this kind of problem be solved? Are there certain approximations that could be made to obtain a closed form likelihood?

I am currently using the pymc3 library to do Bayesian inference. I realize that it is possible to perform a numerical integration to evaluate the convolution of two PDFs, but this is computationally expensive and I'm not sure it would be possible to obtain the gradient of such a function, which would rule out the use of Hamiltonian Monte Carlo sampling (e.g. NUTS).

You are not trying to infer $$\theta_1$$ and $$\theta_2$$ individually (which is easy), so I assume you are attempting to infer them from their convolution, and working backward, given your assumptions on the distributions taken by $$S$$ and $$P$$.