Confidence interval for sum of random variables Let $X_1, X_2$ be independent (not necessarily identically distributed) random variables. Assume we have estimates for the means $m_1, m_2$ and confidence internals $\left[a_1, b_1\right],\left[a_2, b_2\right]$
i.e. $\mathbb{E}\left[X_1\right]=m_1, \mathbb{E}\left[X_2\right]=m_2, \mathbb{P}\left(X_1 \in\left[a_1, b_1\right]\right) \geq 95 \%$ $\mathbb{P}\left(X_2 \in\left[a_2, b_2\right]\right) \geq 95 \%$. We are interested in the sum $Y:=X_1+X_2$. We know that the mean behaves nicely i. e. $\mathbb{E}[Y]=m_1+m_2$. What would be the best way to derive confidence intervals for $Y$ by computational means? What would be the statistical assumptions for each approach be?
On a side note. In my setting, I have many $X_1, \ldots, X_n$. The confidence intervals are symmetric around the mean.
I have thought about assuming gaussian $X_1, \ldots, X_n$ and then reverse-engineering the standard deviation for bootstrap or simply sum of standard deviation or something like that, but did not make it work.

We can't obtain similar intervals for $Y$ without additional information (or making additional assumptions) about the underlying distributions of the $X_i$. For example, consider the case where $\mathbb{P}(X_1\in[-1,1])=95\%$ and $\mathbb{P}(X_2 \in[-2,2])=95\%$. If we assume $X_1\sim N(0,\sigma_1)$ and $X_2\sim N(0,\sigma_2)$, then $X_1+X_2\sim N\Big(0,\sqrt{\sigma_1^2+\sigma_2^2}\Big)$. We can solve for $a_{12}$ and $b_{12}$ where $\mathbb{P}(X_1+X_2\in[a_{12},b_{12}])=95\%$:
ci1 <- c(-1, 1)
ci2 <- c(-2, 2)
hl1 <- diff(ci1)/2
hl2 <- diff(ci2)/2
mu1 <- mean(ci1)
mu2 <- mean(ci2)
# gaussian model
sigma1 <- uniroot(function(sigma) qnorm(0.025, 0, sigma) + hl1, c(0.5, 0.6)*hl1, tol = sqrt(.Machine$double.eps))$root
sigma2 <- sigma1*hl2/hl1
setNames(qnorm(c(0.025, 0.975), mu1 + mu2, sqrt(sigma1^2 + sigma2^2)), c("a12", "b12"))
#>       a12       b12 
#> -2.236068  2.236068

If we instead assume $X_1\sim\text{Cauchy}(0,\gamma_1)$ and $X_2\sim\text{Cauchy}(0,\gamma_2)$, then $X_1+X_2\sim\text{Cauchy}(0,\gamma_1+\gamma_2)$ and $a_{12}=a_1+a_2=-3$ and $b_{12}=b_1+b_2=3$, which are quite different from what was obtained using the Gaussian model.
