Finding the PDF of a Conditional Binomial-Exponential Compound Distribution

Question

I have the following compound binomial-exponential distribution:

$Z \sim \text{Binom}(N,p)$

$p = (1 - e^{-bY})$

$Y \sim \text{Exp}(a)$

$N=1$

I would like to know the p.d.f of the probability distribution $Y$, given that $Z=1$, i.e. $f_{Y|Z=1}$

Attempt at a solution

In a previous question, here, the answerer was able to guide me to the PDF of the distribution $Z$ being binomially distributed, i.e. unmodified by the sampling of the parameter $p$ from the distribution $Y$. So, what I know is:

$P(Z=1) = \frac{b}{a+b}$

or another way of putting it:

$Z \sim \text{Binom}(N, \frac{b}{a+b})$

which is proved in the previous link by both the answerer and I independently.

My first thought is to use Bayes Theorem.

$f_{Y \mid Z}(y \mid z) f_{Z}(z)=f_{Z \mid Y}(z \mid y) f_{Y}(y)$

However, I am unsure of how to proceed from here.

Is it possible that I could be offered a hint as to how to approach this problem? I would like to try and solve it myself, but this is pushing the limit of my existing knowledge of statistics.

Further Information

Provided below is a sample Python program and plot, with a clear example of the nature of the abstraction that I am dealing with.

In this example, I sample two Z distributions, each defined differently, for the $N=1$ case:

$Y \sim \text{Exp}(a)$

$Z_{0} \sim \text{Binom}(N=1, p=(1 - e^{bY}))$

$Z_{1} \sim \text{Binom}(N=1, \frac{b}{(a+b)})$

The expectation values of $Z_{0}$ and $Z_{1}$ may be seen to be the same. However, the conditional distributions are quite different. I am after the orange distribution.

import scipy.stats as stats
import numpy as np
import matplotlib.pyplot as plt


N = int(1e6)

lbda0 = 1
lbda1 = 5

Y = np.array(stats.expon.rvs(scale=1/lbda0, size=N))


        
Z0 = np.array(stats.binom.rvs(
                n=np.ones(N).astype(int), 
                p=lbda1/(lbda0 + lbda1), 
                size=N
        )
)
            
            
Z1 = np.array(stats.binom.rvs(
                n=np.ones(N).astype(int), 
                p=(1 - np.exp(-lbda1*Y)), 
                size=N
        )
)

p_Z1 = 100.0*(np.sum(Z1)/len(Z1))
p_Z0 = 100.0*(np.sum(Z0)/len(Z0))          
    
    
bins = np.linspace(-1, 5, 100)



plt.figure(figsize=(10,10))
plt.title("$N_{{sim}} = {:3.3e}$".format(N), fontsize=20)
plt.hist(Y[Z0==1], bins=bins, histtype="step", density=True, label="P(Y|Z0=1), P(Z0=1)={:3.3f} %".format(p_Z0))
plt.hist(Y[Z1==1], bins=bins, histtype="step", density=True, label="P(Y|Z1=1), P(Z1=1)={:3.3f} %".format(p_Z1))
plt.xticks(fontsize=20)
plt.yticks(fontsize=20)
plt.ylabel("$\\frac{\partial P}{\partial y}$", fontsize=30)
plt.xlabel("$y$", fontsize=20)
plt.legend(fontsize=20)
plt.show()

Answer 1

You have been given $Z\mid Y\sim\mathcal{Bin}(1, (1-\mathrm e^{-bY}))$ and $Y\sim\mathcal{Exp}(a)$.

$$\begin{align}p_{Z\mid Y}(z\mid y)&=(1-\mathrm e^{-by})^z\mathrm e^{-by(1-z)}~[z\in\{0,1\}]\\[1ex]f_Y(y)&=a\mathrm e^{-ay}~[y\in\Bbb R^+]\end{align}$$

So $$f_Y(y)~p_{Z\mid Y}(1\mid y)= a(\mathrm e^{-ay}-\mathrm e^{-(a+b)y})\quad[y\in\Bbb R^+]$$

And, to confirm, we do have: $p_Z(1)=\dfrac{b}{a+b}$ because:

$$\int_0^\infty a (\mathrm e^{-ay}-\mathrm e^{-(a+b)y})~\,\mathrm d y = \dfrac{b}{a+b } $$

---

My first thought is to use Bayes Theorem.

Or just the definition of conditioning: $$\begin{align}f_{Y\mid Z}(y\mid 1)&=\dfrac{f_Y(y)~p_{Z\mid Y}(1\mid y)}{p_Z(1)}\end{align}$$