# Marginal density of $X_1$ given that $X_1 + X_2 = d$ where $X_1$ and $X_2$ are iid Weibull?

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):

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).

library(RColorBrewer)

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")
dev.off()


• 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)$? – Zhanxiong Jun 5 '18 at 13:27
• 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. – 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.

Code in .m:

k = 0.5;
lambda = 1;
d = 10

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

vK = 0.4:0.4:2;

figure(1)
for i=1:length(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);
plot(vX, p_x1_x2/aux)
hold on;
end;


For more accurate graph results, reduce the value of eps

• 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. – Yves Jun 5 '18 at 14:14
• @Yves, right, corrected code and figure – Carlos Campos Jun 5 '18 at 14:15
• Thanks a lot for your answer (+1). Could you explain where the term $f(X_2=d-X_1)$ comes from? – COOLSerdash Jun 6 '18 at 6:56
• @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) – Glen_b Jun 6 '18 at 7:15