Joint probability of two distributions

Question

If I have one random variable that represents hours worked per job X~exponential($\theta$). I have another random variable that represents how many jobs obtained per month Y~Poisson($\lambda$). Using bayes I have used inverse gamma on the exponential and gamma on the Poisson. I now have two gamma distributions. How would I now use those distributions that obtain the probability that the total number of hours worked in a year is Z? p(Z>2000) for example? I imagine I have to use joint probabilities such that $X \times Y > 2000$ and divide that by the total joint probabilities of the two distributions. Given that I have obtained these two distributions? Does anybody know how I would do this in R?

would I be using dgamma for each distribution? How would I then multiply them together to get the total range of possible total hours worked?

I have to multiply all the possible values of one distribution by all the possible values of the other distribution, so I imagine with my two random variables, my function would look like this.

$$f_{XY}(x,y) = xy$$

$$f_{XY}(x,y) = \text{dgamma()*dgamma()}$$ with appropriate parameters?

So I want something like

$$P(X \times Y > 2000) = \int\int f_{XY}(x,y)dxdy$$

$$P(X \times Y > 2000) = \int\int xy \: dxdy$$

Can anyone help?

Answer 1

Total number of hours per year (Z) is a function dependent on the number of hours per job and the number of jobs per months (and we have 12 months).

Z = X(Y(12))

I would use iterate some this expression with your lambdas and rates:

rexp(rpois(12, lambda = 0.5), rate = 0.2)

To estimate the distribution of Z, and then calculate the probability of Z > 2000 from that distribution. Also, as sanity check, the total amount of hours cannot be above 365*24, so I would delete any estimation above 8760 hours of threshold. Example:

iterations <- 1000
Z <- vapply(seq_len(iterations), function(x){sum(rexp(rpois(12, lambda = 0.5), rate = 0.2))}, numeric(1L))
hist(Z)

And then you can calculate the probability with sum(Z >2000)/length(Z)