Expectation and distribution of ratio of correlated Gamma/Chi-square random variables This question is very similar to: Distribution of the ratio of dependent chi-square random variables
But the big difference is what happens when we don't have standard normal variables.

I want to know the expectation (and potentially other moments/distribution) of:
$$ W = \frac{X_1^2}{X_1^2 + \sum_{i=2}^{n} X_i^2} $$
where $X_i \sim N(0, \sigma_i^2)$ for $i = 1, ..., n$.
If all the variances of the $X$-s are equal, then this reduces to a ratio of chi-square random variables and the distribution of $W$ is a beta distribution. However, when the variances are not the same, I can't see a way to use that result.
Is the answer a simpe textbook result, or does it need to be worked out from scratch? If it's the latter, what's the best way to go about it? I haven't done this before.
Another way to phrasing it is what is the distribution of $W = \frac{U}{U+ V}$ where $U \sim \Gamma(\frac{N -1}{2}, 2\sigma^2)^{\text{is this correct?}}$, but then what is $V$? I guess that's another question: if $V = \sum_i \sigma_i ^2 Z_i^2$ is the sum of a bunch of scaled standard normal variables, what is its distribution? Is it some generalised form of the Gamma distribution?

Given some commenters suggested that this may not have a closed form solution, here is an alternative question. With large $n$, $\sum_{i=2}^{n} X_i^2$ tends towards normally distributed. So the alternative question is this:
If we have
$$ W_2 = \frac{Y^2}{Y^2 + a + bZ}$$
where $Y \sim N(0, 1)$ and $Z \sim N(0, 1)$ are independent and standard normal, and $a$ and $b$ are positive constants, does this have a closed form solution?
 A: I don't see where you've described any correlation structure (despite the title using the term "correlated").  So maybe starting small with $n=2$ and assuming independence will be a start.
Using Mathematica (or probably just basic algebra) one can find the distribution of
$$W_2=X_1^2/(X_1^2+X_2^2)$$
where $X_i \sim N(0,\sigma_i^2)$.
dist = TransformedDistribution[x^2, x \[Distributed] NormalDistribution[0, s]];
distw2 = TransformedDistribution[x1/(x1 + x2), 
  {x1 \[Distributed] dist /. s -> \[Sigma][1], 
   x2 \[Distributed] dist /. s -> \[Sigma][2]}];
Mean[distw2]

with the result being
$$\frac{\sigma_1}{\sigma_1+\sigma_2}$$
That suggests that maybe for general $n$, there might be a simple formula for the mean.  Either further manipulations or simulations could give support or squash that hope.
The variance of $W_2$ is $\frac{\sigma_1 \sigma_2}{2 (\sigma_1+\sigma_2)^2}$.
Addition:
The density of $W_3=X_1^2/(X_1^2+X_2^2+X_3^2)$ is
$$-\frac{\sigma_1^2 \sigma_3 E\left(-\frac{(w-1) \sigma_1^2
   \left(\sigma_2^2-\sigma_3^2\right)}{\left((w-1) \sigma_1^2-w
   \sigma_2^2\right) \sigma_3^2}\right)}{\pi  \sqrt{w \left(\sigma_2^2
   w-\sigma_1^2 (w-1)\right)} \left(\sigma_1^2 (w-1)-\sigma_3^2 w\right)}$$
where $\sigma_2 \leq \sigma_3$, $E(.)$ is the complete elliptic integral, and found with the following Mathematica code:
dist2 = TransformedDistribution[x2 + x3, 
  {x2 \[Distributed] GammaDistribution[1/2, 2 s[2]^2], 
   x3 \[Distributed] GammaDistribution[1/2, 2 s[3]^2]},
   Assumptions -> 0 < s[2] <= s[3]];

dist123 = TransformedDistribution[x1/(x1 + x23), 
  {x1 \[Distributed] GammaDistribution[1/2, 2 s[1]^2],
   x23 \[Distributed] dist2}];
Simplify[PDF[dist123, w], Assumptions -> s[1] > 0]

I have not found a closed-form representation for the mean when $n=3$ but when the variances are known, then numerical integration should easily find the mean and other desired moments.
Special case
Suppose $\sigma_1^2=\sigma^2_U$ and $\sigma_i^2=\sigma^2_V$ for $i=2,\ldots,n$.  Then $W=U/(U+V)$.  The mean of that distribution can be found with
dist = TransformedDistribution[
  u/(u + v), {u \[Distributed] GammaDistribution[1/2, 2 \[Sigma]u^2],
   v \[Distributed] GammaDistribution[(n - 1)/2, 2 \[Sigma]v^2]}]

mean = FullSimplify[Mean[dist], Assumptions -> {\[Sigma]u > 0, \[Sigma]v > 0}]

$$\frac{(n+1) \sigma_u^2 \, _2F_1\left(-\frac{1}{2},\frac{n}{2};\frac{n+2}{2};1-\frac{\sigma_v^2}{\sigma_u^2}\right)-\left((n-1) \sigma_v^2+\sigma_u^2\right) \, _2F_1\left(\frac{1}{2},\frac{n}{2};\frac{n+2}{2};1-\frac{\sigma_v^2}{\sigma_u^2}\right)}{n \sigma_u \sigma_v}$$
where $_2 F_1$ is the hypergeometric function.
Specifying values of $n$ results in great simplifications:
FullSimplify[mean /. n -> 2, Assumptions -> {\[Sigma]u > 0, \[Sigma]v > 0}]

$$\frac{\sigma_U}{\sigma_U+\sigma_V}$$
FullSimplify[mean /. n -> 6, Assumptions -> {\[Sigma]u > 0, \[Sigma]v > 0}]

$$\frac{\sigma_U^2 (3 \sigma_U+\sigma_V)}{3 (\sigma_U+\sigma_V)^3}$$
