In their tutorial (page 23) on heavy-tailed distributions, Nair et al. present the following graph (taken from a pre-publication chapter from a book by the same authors):

Conditional distribution of Weibull

Pictured are the marginal densities of $X_1$ given that $X_1 + X_2 = 10$, where $X_1$ and $X_2$ are iid Weibull r.v.s (either with shape 2 or 0.5, the scale is not provided by the authors).

Unfortunately, I can't reproduce how the authors calculated the marginal density. According to my internet search, the sum of Weibull variables is generally difficult to calculate (cf. Nadarajah S (2008)).

My question is therefore: Can someone show how the authors calculated the marginal density of $X_1$, given that $X_1 + X_2 = d$?

Edit: Many thanks to @CarlosCampos for his solution and illustration. I ported his Matlab code to R (added below).


k <- 0.5
lambda <- 1
d <- 10

eps <- 0.001
vX <- seq(eps, d-eps, eps)
vK <- seq(0.4, 2, 0.4)

cols <- brewer.pal(length(vK), name = "Set1")

png("weibull_marginal.png", width = 7*1.1, height = 5*1.1, units = "in", res = 600, type = "cairo-png")
par(cex = 1.2, mar=c(2,2.1,0,0)+0.1)
plot(x = vX, y = rep(1, length(vX)), las = 1, type = "n", xlab = "", ylab = "", xlim = c(0, d), ylim = c(0, 1))
abline(h = 0)
for(i in seq_along(vK)) {

  k <- vK[i]
  p_x1_x2 <- (k/lambda)^2*(vX*(d - vX)/(lambda^2))^(k - 1)*exp(-(vX^k + (d - vX)^k)/lambda^k)
  aux <- sum(p_x1_x2*eps)
  p_x1_x2_norm <- p_x1_x2/aux
  lines(p_x1_x2_norm~vX, col = cols[i], lwd = 2)

par(cex = 1)
legend(x=8, y=1.03, legend = vK, col = cols, lwd = 2, box.lwd = 0, box.col = "white", bg = "white")

Reproduced Graph using R

  • $\begingroup$ Have you considered making the 1-1 transformation $Y_1 = X_1, Y_2 = X_1 + X_2$, then find the joint density of $(Y_1, Y_2)$ from that of $(X_1, X_2)$? $\endgroup$ – Zhanxiong Jun 5 '18 at 13:27
  • 2
    $\begingroup$ Nice reference. You can simply evaluate the density of the $X_i$ on a regular grid and use a numeric convolution to find the density of $X_1 + X_2$; Then the conditional density results. The precision is limited, but may be acceptable for an illustration. $\endgroup$ – Yves Jun 5 '18 at 13:45

Let's apply Bayes theorem: $$f(X_1 \vert X_1+X_2 = d) = \frac{f(X_1+X_2 = d \vert X_1)f(X_1)}{f(X_1+X_2 = d)} = cte \cdot f(X_2=d-X_1)f(X_1)$$ Substituting expressions for Weibull distributions: $$f(X_1 \vert X_1+X_2 = d) = cte \cdot \left(\frac{k}{\lambda}\right) \left(\frac{d-x_1}{\lambda}\right)^{k-1} e^{((d-x_1)/\lambda)^k}\left(\frac{k}{\lambda}\right) \left(\frac{x_1}{\lambda}\right)^{k-1} e^{(x_1/\lambda)^k}$$ for $0<x_1<d$, otherwise $f(x)=0$ Where $cte$ can be easily evaluated. I provide a .m code to generate that figure.

enter image description here

Code in .m:

k = 0.5;
lambda = 1;
d = 10 

eps = 0.1;
vX = eps:eps:d-eps;

vK = 0.4:0.4:2;

for i=1:length(vK)
  p_x1_x2 = (k/lambda)^2*(vX.*(d-vX)/(lambda^2)).^(k-1).*exp(-(vX.^k+(d-vX).^k)/lambda^k);
  aux = sum(p_x1_x2*eps);
  plot(vX, p_x1_x2/aux)
  hold on;

For more accurate graph results, reduce the value of eps

  • 2
    $\begingroup$ Minor pb: you should divide by aux * eps instead of aux to get the right normalisation. The densities on the original plot are not either suitably normalised. $\endgroup$ – Yves Jun 5 '18 at 14:14
  • $\begingroup$ @Yves, right, corrected code and figure $\endgroup$ – Carlos Campos Jun 5 '18 at 14:15
  • $\begingroup$ Thanks a lot for your answer (+1). Could you explain where the term $f(X_2=d-X_1)$ comes from? $\endgroup$ – COOLSerdash Jun 6 '18 at 6:56
  • 2
    $\begingroup$ @COOLSerdash The product of $f$'s because of independence. Then you just use the fact that if $X_1+X_2 = d$ then $X_2=d-X_1$, so once $d$ is fixed, the integral is a function of only $X_1$; this is just a completely standard convolution integral en.wikipedia.org/wiki/Convolution_of_probability_distributions (however note we're abusing notation here by conflating the random variables with the values they take) $\endgroup$ – Glen_b -Reinstate Monica Jun 6 '18 at 7:15

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.