ML Interview Q Series: Density of Random Variable Product XY via Cumulative Distribution Function Integration.

Browse all the Probability Interview Questions here.

11E-12. The joint density function of the random variables X and Y is given by f(x,y) = x e^(-x(y+1)) for x>0,y>0, and 0 otherwise. What is the density function of the random variable XY?

Short Compact solution

From the integral calculation of P(XY ≤ z) over the region 0 < x ≤ z and 0 < y ≤ z/x (for z>0), and then differentiating with respect to z, the result is:

and f_Z(z) = 0 otherwise.

Comprehensive Explanation

Overview of the approach

We have two random variables X and Y with joint pdf f(x,y) = x e^(-x(y+1)), defined for x>0 and y>0. We define Z = X·Y and want the pdf of Z, denoted f_Z(z).

A common strategy to find the distribution of a function of two random variables is:

1. Express P(Z ≤ z) = ∫∫ I(x y ≤ z) f(x,y) dx dy, where I(...) is the indicator function of the region x·y ≤ z.

2. Differentiate F_Z(z) = P(Z ≤ z) to obtain f_Z(z).

Because X>0 and Y>0, we only integrate over the positive quadrant. The condition x y ≤ z further constrains Y ≤ z/x (with x>0). Hence, for z>0:

F_Z(z) = P(Z ≤ z) = ∫(x from 0 to z) ∫(y from 0 to z/x) x e^(-x(y+1)) dy dx.

We then differentiate F_Z(z) with respect to z to find f_Z(z).

Detailed derivation

1. Set up the integral for P(Z ≤ z):

   F_Z(z) = ∫[0 to z] ∫[0 to z/x] x e^(-x(y+1)) dy dx, (for z>0).

2. Perform the inner integral over y:

   • The inner integral ∫[0 to z/x] x e^(-x(y+1)) dy = x ∫[0 to z/x] e^(-x(y+1)) dy.

   • Let’s integrate w.r.t y: ∫ e^(-x(y+1)) dy = (1 / (-x)) e^(-x(y+1)) (evaluated from 0 to z/x).

   • Carefully evaluating: e^(-x(y+1)) from y=0 to y=z/x gives e^(-x((z/x)+1)) - e^(-x(0+1)). That is e^(-z - x) - e^(-x).

   • Multiplying by x / (-x) = -1, we get: x [ (e^(-z - x) - e^(-x)) / (-x) ] = e^(-x) - e^(-z - x).

   So the inner integral = e^(-x) - e^(-z - x).

3. Now the outer integral:

   F_Z(z) = ∫[0 to z] ( e^(-x) - e^(-(z + x)) ) dx = ∫[0 to z] e^(-x) dx - ∫[0 to z] e^(-(z + x)) dx.

   • First term: ∫[0 to z] e^(-x) dx = 1 - e^(-z).

   • Second term: ∫[0 to z] e^(-(z + x)) dx = e^(-z) ∫[0 to z] e^(-x) dx = e^(-z) [1 - e^(-z)].

   Substituting these, we get: F_Z(z) = (1 - e^(-z)) - e^(-z) [1 - e^(-z)] = (1 - e^(-z)) (1 - e^(-z)) = (1 - e^(-z))^2.

4. Differentiate F_Z(z):

   f_Z(z) = d/dz [ (1 - e^(-z))^2 ] = 2 (1 - e^(-z)) · d/dz [1 - e^(-z)] = 2 (1 - e^(-z)) ( e^(-z) ) = 2 e^(-z) (1 - e^(-z)).

5. Domain: Because X and Y are strictly positive, Z = X·Y is also strictly positive. Thus f_Z(z) = 2 e^(-z) (1 - e^(-z)) for z>0, and f_Z(z) = 0 for z ≤ 0.

Hence the final pdf of Z = X·Y is:

Why this pdf makes sense

• As z→0+, we expect the density to be 0 because X·Y cannot be less than or equal to a value extremely close to 0 unless X or Y is very small. Indeed, 2 e^(-0) (1 - e^( - 0)) = 0, so it matches the intuition of the limit at z=0+.

• As z grows large, the factor e^(-z) drives the pdf to 0, ensuring the distribution has finite total integral over (0,∞).

Potential alternative approaches

• One can also attempt a transformation approach with U = X, V = XY, and then find the joint distribution of (U, V) and integrate out U. However, the direct cumulative distribution approach here is typically more straightforward.

• Sometimes a Laplace transform or moment generating function approach can be used for the product XY, but that might be more involved than the direct integral.

Possible Follow-up Questions

1) Could we have used a change of variables?

Yes. Define U = X and V = X·Y. Then Y = V / U. For x>0 and y>0, that means u>0 and v>0. One would compute the Jacobian of the transformation, find the joint distribution of (U, V), and then integrate out U to get f_V(v). The integral setup can become very similar, but the direct approach used above is often more straightforward.

2) How do we verify this is a valid pdf for Z?

To verify:

1. Check non-negativity: 2 e^(-z) (1 - e^(-z)) ≥ 0 for z>0, which is true because both e^(-z) > 0 and (1 - e^(-z)) ≥ 0 for z≥0.

2. Integrate over the domain (0,∞) and ensure it integrates to 1.

Quick outline of the integral check: ∫[0 to ∞] 2 e^(-z)(1 - e^(-z)) dz. We can split it: 2 ∫[0 to ∞] e^(-z) dz - 2 ∫[0 to ∞] e^(-2z) dz = 2[1] - 2[1/2] = 1. So it indeed integrates to 1, confirming that f_Z(z) is a valid pdf.

3) What if z <= 0?

Because X and Y are strictly positive, Z = XY > 0 always. Hence the pdf for z ≤ 0 must be 0. This is consistent with the domain of the random variable Z.

4) Practical implementation or simulation in Python

To simulate samples from Z = X·Y when X and Y have the given joint density, one straightforward approach is:

1. Recognize that marginal distribution of X can be found, then condition on X to find Y.

2. Draw X from its marginal distribution.

3. Draw Y from the conditional distribution given X.

4. Multiply X and Y to get Z.

For example:

import numpy as np

def sample_XY(n_samples=10_000):
    # Step 1: Find marginal of X and sample X
    # The marginal of X is ∫[0 to ∞] x e^(-x(y+1)) dy
    # = ∫[0 to ∞] x e^(-x) e^(-x y) dy over y>0
    # By standard integrals, the marginal of X is x e^(-x) * ∫[0 to ∞] e^(-x y) dy
    # = x e^(-x) * (1 / (x)) = e^(-x), for x>0.
    # So X ~ Exponential(1). We can sample from np.random.exponential(scale=1, size=n_samples).

    X = np.random.exponential(scale=1.0, size=n_samples)

    # Step 2: Given X=x, the conditional pdf of Y is:
    # f_{Y|X=x}(y) = [x e^(-x(y+1))] / [marginal of X at x]
    # = [x e^(-x(y+1))] / [e^(-x)] = x e^(-x y)
    # This is the pdf of an exponential with parameter x.
    # So Y|X=x ~ Exponential(x).

    Y = np.array([np.random.exponential(scale=1/x_i) for x_i in X])

    # Return samples
    return X, Y

def sample_Z(n_samples=10_000):
    X, Y = sample_XY(n_samples)
    Z = X * Y
    return Z

# Example usage:
Z_samples = sample_Z(100000)
# We can then plot or check histogram to compare with the theoretical pdf 2 e^(-z)(1 - e^(-z)).

This simulation approach verifies that if you histogram Z and compare it to 2 e^(-z)(1 - e^(-z)) for z>0, you should see a close match as sample size grows large.

5) Could we interpret Z = XY in a real-world scenario?

Yes. If X and Y were, for example, waiting times in certain exponential-related processes, their product might arise in models of “product of waiting times.” The derived distribution can appear in specialized reliability or queueing contexts. The ability to derive it ensures we understand how to manipulate distributions for transformations in theoretical or applied settings.

All of these considerations demonstrate not just the final formula but the deeper reasoning and checks that confirm correctness.

Below are additional follow-up questions

1) How can we compute the expected value of Z = XY directly from its PDF, and are there any alternative methods to find E[Z]?

One way is to use the derived pdf of Z directly. Recall that Z takes values in (0,∞), and its pdf is f_Z(z) = 2 e^(-z) (1 - e^(-z)) for z>0. Then:

E[Z] = ∫(from 0 to ∞) z f_Z(z) dz = ∫(from 0 to ∞) z (2 e^(-z) (1 - e^(-z))) dz.

We can split this integral:

E[Z] = 2 ∫(from 0 to ∞) z e^(-z) dz - 2 ∫(from 0 to ∞) z e^(-2z) dz.

Each integral can be computed using known results for the gamma function or by integration by parts:

• ∫(0 to ∞) z e^(-az) dz = 1 / a^2 (for a>0).

Hence:

• ∫(0 to ∞) z e^(-z) dz = 1.

• ∫(0 to ∞) z e^(-2z) dz = 1/4.

Substitute back:

E[Z] = 2·(1) - 2·(1/4) = 2 - 1/2 = 3/2.

An alternative approach is to use knowledge of X and Y directly. We can compute E[Z] = E[X Y] by noticing that X and Y are not independent, but we can try to use the law of total expectation:

E[X Y] = E[ X E(Y|X) ].

We know from the joint pdf that the marginal of X is e^(-x) for x>0 (i.e., X ~ Exponential(1)). Then given X=x, Y ~ Exponential(x). Hence E[Y|X=x] = 1/x. So:

E[X Y] = E[ X (1/x) ] = E[1] = 1,

which might seem contradictory. However, caution: one must carefully verify the correct conditional distribution. Actually, from the joint distribution f(x,y) = x e^(-x(y+1)), we do have Y|X=x ~ Exponential(x). But we must confirm:

E[Y|X=x] = 1/x, thus E[X·Y] = E[X(1/x)] = E[1], and that yields 1 if X is a positive random variable taking values in (0,∞). But we need to confirm the exact distribution of X. If X is indeed Exponential(1), E[1/x] is not 1/X in a naive sense. In fact, for X~Exponential(1), E[1/X] is infinite. That means a direct approach to E[X Y] via conditioning can get tricky due to the tail behavior of 1/X.

Hence a more robust approach is to rely on the direct pdf of Z or the double integral:

E[XY] = ∫∫ x·y f(x,y) dx dy, making sure all integrals converge. In this problem, the easiest correct path is to use the derived f_Z(z) and integrate. That yields 3/2.

The key pitfall is that E[1/X] diverges when X ~ Exp(1), so one cannot just do X * E[Y|X] = X*(1/X) in expectation without carefully handling that X is random. It’s a subtle demonstration that direct conditioning can lead to formal manipulations that don’t hold if integrals are not finite.

2) Are there any potential convergence or numerical issues when estimating this distribution from simulated data?

Yes. A few typical pitfalls:

• Small x values: Because X is Exponential(1), there is a non-negligible probability that X will be very small. Then Y|X=x ~ Exp(x) can have a fairly large mean for tiny x, potentially yielding numerical instability if one tries to invert or transform these samples incorrectly.

• Large product: Sometimes, if you sample a moderately large X and a moderately large Y, their product Z = XY can become quite large, and in floating-point arithmetic you might lose some precision when dealing with exponentials of large negative arguments (e^(-z) might underflow for large z).

• Histogram binning: When building a histogram for the product Z, one must choose the bin ranges carefully. If you pick bins that are too coarse near z=0, you might fail to capture the quick rise in the pdf. If bins are too coarse in the tail, you might get poor estimates for the distribution tail.

A good practical practice is:

• Use double-precision or higher to reduce floating-point issues.

• Use appropriate transformations if you suspect extremely large or extremely small values (e.g., analyzing log(Z) can sometimes help).

• Ensure the number of samples is large enough to get stable estimates in both the lower and upper tails.

3) Could this product distribution Z = X·Y be related to known “product of exponential” distributions or generalized gamma distributions?

Yes, in certain cases. The distribution for a product of two independent exponential random variables with parameters λ1 and λ2 has known forms. However, in our problem, X and Y are not independent (their joint pdf couples them). If X and Y were independent Exp(1), the product XY would have a known pdf involving the modified Bessel function of the second kind (related to the distribution known as the product of two independent exponential(1)). But here, the presence of x in the exponent and the extra factor x in f(x,y) indicates correlation.

In short, while it looks similar to a “product-of-exponentials” scenario, we can’t simply apply standard product-of-independent-Exponential distributions. We found a different pdf, 2 e^(-z)(1 - e^(-z)), which is unique to the specific correlation structure in f(x,y).

4) What is the cumulative distribution function F_Z(z) for Z = XY, and can it help in bounding or confidence interval calculations?

We found that:

F_Z(z) = (1 - e^(-z))^2 for z>0.

Hence if you want to find a probability of the form P(Z ≤ c), you use (1 - e^(-c))^2. For bounding or confidence intervals, you might invert that relationship:

For a desired quantile q in (0,1), solve:

(1 - e^(-z))^2 = q,

which implies:

1 - e^(-z) = sqrt(q), e^(-z) = 1 - sqrt(q), z = - ln(1 - sqrt(q)).

That expression yields the z-value at which the distribution function is q. This can be used for confidence bounds (e.g., a 95% upper bound for Z) or acceptance thresholds in hypothesis testing. A subtlety is that (1 - sqrt(q)) must remain positive, so q must be strictly less than 1 for that expression to make sense.

5) How would the distribution change if the domain of X or Y were altered (e.g., X>0 but Y≥-1/x to satisfy positivity in the exponent)? Could that cause confusion?

Yes, boundary changes can drastically alter the shape of the joint pdf. Our problem statement explicitly says X>0 and Y>0, so that’s a strict domain. If someone misread the pdf as x e^(-x(y+1)) for all real y, they might get nonsensical integrals or domain definitions. The positivity of y is also crucial in ensuring the integrals converge the way we derived.

If we allowed Y to be negative, that changes the region of integration for x y ≤ z (especially if z could also be negative). That would lead to entirely different expressions and possibly break the standard approach or complicate it significantly. So it’s essential to keep the domain constraints in mind.

6) What if we were asked about higher moments, like E[Z^n] for n>1, or the moment generating function (MGF) or characteristic function of Z?

To find E[Z^n] directly from the pdf, we would compute:

E[Z^n] = ∫(0 to ∞) z^n (2 e^(-z)(1 - e^(-z))) dz.

We can break that into:

2 ∫(0 to ∞) z^n e^(-z) dz - 2 ∫(0 to ∞) z^n e^(-2z) dz.

Each integral can be computed using the gamma function:

∫(0 to ∞) z^n e^(-az) dz = n! / a^(n+1), provided n is a nonnegative integer and a>0.

Hence:

E[Z^n] = 2 [ n! - n!/2^(n+1) ] = 2 n! (1 - 1/2^(n+1)).

Be aware that for large n, n! can grow very quickly, so E[Z^n] can be quite large.

For the MGF M_Z(t) = E[e^(tZ)], we would calculate:

M_Z(t) = ∫(0 to ∞) e^(t z) [2 e^(-z)(1 - e^(-z))] dz,

which can be split similarly. However, the integrals become:

2 ∫(0 to ∞) e^((t-1)z) dz - 2 ∫(0 to ∞) e^((t-2)z) dz,

provided the real part of t<1 and t<2 for those integrals to converge. If we assume t<1:

• The first integral is 2 / (1 - t).

• The second integral requires t<2, giving 2 / (2 - t).

Hence:

M_Z(t) = (2 / (1 - t)) - (2 / (2 - t)), for t<1 (and also t<2). One would simplify that expression for a more direct formula, but it shows how the MGF can exist only in a bounded region of t.

7) Could f_Z(z) be used to model real-world processes and how to validate it with empirical data?

Yes. Although derived in a theoretical problem, there may be situations where a random product arises naturally—for instance, X might represent a rate parameter that is itself random, and Y is a waiting time that depends on that random rate. This situation can occur in hierarchical or multi-stage processes (e.g., a random ‘failure rate’ X chosen first, then an actual waiting time Y based on that rate).

To validate, one would:

1. Collect pairs (x_i, y_i) from the system or measure only z_i = x_i y_i if partial data is missing.

2. If both x_i and y_i are available, check whether the joint distribution assumption x e^(-x(y+1)) is consistent with the data. This can be done by maximum likelihood or a goodness-of-fit test in two dimensions.

3. If only z_i = x_i y_i is recorded, estimate the pdf of z_i via a histogram or kernel density estimator and see if it matches 2 e^(-z)(1 - e^(-z)) for z>0. For more robust inference, one can attempt parametric fitting and see if the derived distribution is a good fit.

Potential pitfalls:

• If the real process has a different correlation structure between X and Y, the above formula will not match.

• If the domain (x>0,y>0) is not appropriate for the real data (maybe y can be zero or negative), that mismatch must be accounted for.

• In practice, small-sample issues can make it difficult to differentiate between multiple candidate distributions. A thorough approach includes statistical tests, residual analysis, or information criteria to decide if this distribution is indeed appropriate.

8) In what scenario would the transformation technique (U=X, V=XY) become cumbersome, and how do we address that?

The transformation approach requires careful computation of the Jacobian and setting up integration limits in the (u, v) plane. Specifically:

• U = X, V = XY => Y = V / U.

• The support region for (U, V) can get complicated. Here, U>0 and V>0, but the condition that Y>0 means V = u * y>0 only if u>0 and y>0. That implies v>0 and also we have to carefully specify 0<u<∞, 0<v<∞.

• The Jacobian is |det| = |∂(x,y)/∂(u,v)| = 1/u.

We would get:

f_{U,V}(u, v) = f_{X,Y}(u, v/u) * (1 / u),

but we must remain cautious about domain constraints so that v/u>0 => v>0, etc. Once that is derived, to find f_V(v), we integrate over u>0:

f_V(v) = ∫(0 to ∞) f_{U,V}(u,v) du.

While this is perfectly correct, each step must be done carefully to avoid sign or domain mistakes. For complex problems or piecewise domains, the transformation can become cumbersome. Address it by double-checking domain boundaries, sketching (u,v) ranges, and verifying that for each (x,y), the transformation is invertible in the region. Also, watch out for singularities if u=0 or v=0 can appear in the domain.

In summary, while transformation is powerful, the region-of-integration complexities and potential for algebraic slips are typical pitfalls. Sometimes, the direct approach of computing P(Z ≤ z) = ∫∫ I(xy ≤ z) f(x,y) dx dy might be more straightforward, as was done here.

9) Could the random variable Z = XY exhibit any interesting shape features like a unimodal pdf or monotonic hazard function?

• Shape: The pdf f_Z(z)=2 e^(-z)(1 - e^(-z)) starts at 0 when z=0, initially rises to a peak, and then decays toward 0 as z→∞. That is a unimodal shape. In fact, the mode can be found by setting the derivative of 2 e^(-z)(1 - e^(-z)) to zero.

• Hazard function: The hazard rate h_Z(z) = f_Z(z) / (1 - F_Z(z)). With F_Z(z) = (1 - e^(-z))^2, 1 - F_Z(z) = 1 - (1 - e^(-z))^2 = 2 e^(-z) - e^(-2z). Substituting:

  h_Z(z) = [2 e^(-z)(1 - e^(-z))] / [2 e^(-z) - e^(-2z)].

  After simplifying, one could study how it behaves as z→0 or z→∞. This might be relevant in reliability theory if we interpret Z as a “time to failure.”

Such shape analyses can matter in real-world modeling if we suspect a single peak in probability density, or if we want to characterize how “failure rate” changes over time. The hazard might not be monotonic here, so it’s a subtle distribution from a reliability perspective.

10) What if we wanted the distribution of 1/Z instead of Z, or other functions like Z^(-1/2)?

Sometimes in Bayesian or reliability contexts, we might need to consider the reciprocal 1/Z. You can do so by the standard transformation approach for a single variable:

If W = 1/Z, then Z = 1/W. The pdf transformation formula for a one-to-one map W -> Z is:

f_W(w) = f_Z(z(w)) · |d z(w)/d w|,

with z(w)=1/w. Then dz/dw = -1/w^2, so |dz/dw|=1/w^2. So:

f_W(w) = f_Z(1/w) * (1 / w^2),

for w>0 (since Z>0 => W>0). Substituting f_Z(1/w)=2 e^(-(1/w))(1 - e^(-(1/w))) for w>0. You’d end up with:

f_W(w) = 2 (1/w^2) e^(-1/w) (1 - e^(-1/w)).

One must be careful with the domain and potential numerical stability if w is very large or very small. Also, integrals with e^(-1/w) can have tricky expansions near w=0. So, while straightforward in principle, it’s important to handle boundary behaviors carefully.

These details highlight that any function of Z can be studied, but the transformations require careful domain and Jacobian considerations.