En tu pregunta impones la condición de independencia por pares, pero no es posible encontrar la distribución de una combinación lineal de más de dos variables aleatorias con esta especificación. Por lo tanto, supondré que estás dispuesto a imponer la condición más fuerte de independencia mutua para que la distribución de la combinación lineal esté totalmente determinada.
Sobre esta base, lo que se busca es la familia de distribuciones estables que se puede obtener considerando las propiedades de funciones características (CF) (es decir, las transformaciones de Fourier de las funciones de densidad). Si $X_1, ..., X_n$ son variables aleatorias mutuamente independientes con funciones características denotadas por $\varphi_k(t) = \mathbb{E}(\exp(it X_k))$ entonces cualquier combinación lineal de estas variables aleatorias tiene una CF que es un producto simple de las CF de las variables aleatorias subyacentes:
$$Y = \sum_{k=1}^n \alpha_k X_k \quad \quad \implies \quad \quad \varphi_Y(t) = \prod_{k=1}^n \varphi_k(\alpha_k t).$$
Lo que se busca son clases de distribuciones en las que la forma de la distribución es cerrada bajo combinaciones lineales. Esta propiedad se da en las clases de distribuciones en las que la forma de la función característica se cierra bajo la multiplicación. Esto ocurre para la distribución normal como caso especial, pero también ocurre para una familia más amplia de distribuciones llamada distribución alfa-estable de Lévy (ver más abajo).
Distribución normal: Esto tiene la función característica $\varphi_k(t) = \exp(i \mu_k t - \tfrac{1}{2} \sigma_k^2 t^2)$ . Podemos ver que esta forma funcional es cerrada bajo la multiplicación:
$$\begin{equation} \begin{aligned} \varphi_Y(t) = \prod_{k=1}^n \varphi_k(\alpha_k t) &= \prod_{k=1}^n \exp(i \alpha_k \mu_k t - \tfrac{1}{2} \alpha_k^2 \sigma_k^2 t^2) \\[8pt] &= \exp \Big( \sum_{k=1}^n (i \alpha_k \mu_k t - \tfrac{1}{2} \alpha_k^2 \sigma_k^2 t^2) \Big) \\[8pt] &= \exp \Big( i \mu_Y t - \tfrac{1}{2} \sigma_Y^2 t^2 \Big), \\[8pt] \end{aligned} \end{equation}$$
donde $\mu_Y = \sum_{k=1}^n \alpha_k \mu_k$ y $\sigma_Y^2 = \sum_{k=1}^n \alpha_k \sigma_k^2$ . A partir de esta función característica vemos que la combinación lineal de variables aleatorias también es normal, con media y varianza actualizadas que son las correspondientes combinaciones lineales de las medias y varianzas subyacentes.
Distribución alfa-estable de Lévy: Se trata de una generalización de la distribución normal, que permite una función característica generalizada que sigue siendo cerrada bajo la multiplicación. Tiene una función característica definida por:
$$\varphi_k(t) = \exp(i \mu_k t - \tfrac{1}{\gamma} |\sigma_k t|^{\gamma} (1-ib_k \cdot \text{sgn}(t) \Phi_{\gamma}(t))),$$
donde $0 < \gamma \leqslant 2$ es un parámetro fijo (común a cada una de las distribuciones subyacentes) y:
$$\Phi_\gamma(t) = \begin{cases} \tan (\tfrac{\pi \gamma}{2}) & & \text{for } \gamma \neq 1, \\[4pt] - \tfrac{2}{\pi} \ln |t| & & \text{for } \gamma = 1. \end{cases}$$
(También tenemos la restricción del parámetro $|b_k| \leqslant 1$ .) Esta función característica también es cerrada bajo la multiplicación (cuando el parámetro $\gamma$ es común a todas las variables aleatorias subyacentes), por lo que define una clase de distribuciones más amplia que la distribución normal. La distribución normal se da en el caso especial de que $\gamma = 2$ y $b_k = 0$ . La función de densidad correspondiente a esta distribución no existe en forma cerrada, pero es posible derivar muchas de sus propiedades (que dejaré fuera del alcance de la presente respuesta).