Asymptotic behaviour of a multiple integral on the unit hypercube

Let $X_j$ be a sequence of i.i.d. uniform(0,1) random variables, and set $Y_n=\sum_{j=1}^n X^2_j/\sum_{j=1}^n X_j$. The integral in the problem can be expressed as $I_n=\mathbb{E}(Y_n)$.

Note that $0\leq Y_n\leq 1$ almost surely, and that by the strong law of large numbers $$Y_n={\sum_{j=1}^n X^2_j \over \sum_{j=1}^n X_j}= {\sum_{j=1}^n X^2_j \over n}\times {n\over \sum_{j=1}^n X_j}\to {1\over 3}\times{2\over 1},$$ almost surely. By the dominated convergence theorem, we conclude that $I_n\to 2/3$.

In fact I fully agree with Byron Schmuland. You can see easily that:

$$ \frac{2}{3}=\frac{\int_{(0,1)^n} \sum_{i=1}^{n}x_i^2}{\int_{(0,1)^n}\sum_{i=1}^{n}x_i}\leq\int_{(0,1)^n} \frac{\sum_{i=1}^n x_i^2}{\sum_{i=1}^n x_i} d\mu $$

[Warning: after an hint from @Biron, I edit this. From Chebyshev we can just obtain a lower bound of $\frac{1}{2}$]

$$ \int_{(0,1)^n} \frac{\sum_{i=1}^n x_i^2}{\sum_{i=1}^n x_i} d\mu \leq\int_{(0,1)^n} \sum_{i=1}^n x_i^2d\mu \int_{(0,1)^n}\frac{1}{\sum_{i=1}^n x_i}d\mu = \frac{1}{3} \int_{(0,1)^n}\frac{n \space d\mu}{\sum_{i=1}^n x_i} $$

These are consequences of Chebyshev's integral inequality. But in the last integral you can use the new variable: $$s=\frac{1}{n}\sum_{i=1}^{n}x_i$$ and express the measure: $$d \mu^{(n)} = \frac{d \mu^{(n)}}{ds}ds=p^{(n)}(s)ds$$ where you know $p^{(n)}(s)$ is the multiple convolution of $$p^{(1)}(x)=\chi_{(0,1)}(x)$$ and using tauberian calculus and Laplace transform you can discover that this converges to: $$p^{(n)}(s)=\frac{\sqrt{n}}{\sqrt{2\pi}\sigma}e^{-\frac{n(s-\bar{x})}{2\sigma^2}^2}$$ where: $$\sigma^2=\mathbb{E}[(x-\bar{x})^2]=\frac{1}{12}$$ using the Chebyshev's inequalities you can now easily obtain a good inequality for an upper limit of the integral:

$$\int_{(0,1)^n}\frac{n \space d\mu}{\sum_{i=1}^n x_i} \leq 2 + O\left(\frac{1}{n}\right)$$

so that integral converges to a limit not greater than: $$\frac{2}{3}$$

[After I wrote this I obtained a proof that the limit is $2/3$ following a different approach]

In order to gain a finer knowledge of the asymptotic behavior of $I_n$ you can write down: $$K^{(n)}=\int_{(0,1)^n}\frac{n \space d\mu}{\sum_{i=1}^n x_i}=\int_0^1 \frac{p^{(n)}(s)}{s}ds=2\int_{-1/2}^{1/2}\frac{\bar{p}^{(n)}(x)}{1-2x}dx $$ where $s=1/2+x$ and: $$\bar{p}^{(n)}(x)=p^{(n)}(1/2+x)$$ so that: $$K^{(n)} = 2 \sum_{i=0}^{\infty}\mathbb{E_{p^{(n)}(x)}}[(2x)^{2i}]\simeq 2 + 2\sum_{i=1}^{\infty} \left(\frac{1}{3n}\right)^i(2i-1)!!$$

where the geometric series has been used, and where, because the distribution is pair, the odd moments are missing. In the last expression you can see the gaussian extimation of the moments, however this expression doesn't converge. In fact, increasing the exponent of $x^{2i}$ the peaks of function $x^{2i}\frac{1}{\sqrt{2\pi\sigma}}e^{-\frac{x^2}{2\sigma^2}}$ are $\pm\sigma\sqrt{2i}$, so the gaussian extimation is accurate only if the peaks are in the range of variation of our variable x: $\frac{\sqrt{2i}}{\sqrt{12n}} \ll \frac{1}{2} $ or, in other terms: $i \ll 6n$. In order to obtain a non diverging extimation of moments for bigger value of i a way is to consider the partial moments of gaussian approximation in the restricted range $(-\frac{1}{2}+\frac{1}{2(n-1)},\frac{1}{2}-\frac{1}{2(n-1)})$ and to consider the right volume in the range $(-\frac{1}{2},-\frac{1}{2}-\frac{1}{2n})$ which amount exactly to $\frac{(x+1/2)^{n-1}}{(n-1)!}dx$. The asymptotic behavior of the partial moments results now upper bounded in fact, if we assume the meausure concentrate in extremum of interval, we obtain: $$\mathbb{E}_{p^{(n)}(x)}[x^{2i}] \leq 2\left(\frac{1}{2}-\frac{1}{2n}\right)^{2i}$$ and the contribute of the two symmetric symplectic range in the hypercube is bounded from: $$ \frac{1}{n!(2n)^n}+\frac{1}{2(n-1)^{n-1}(n-1)(n-1)!} $$

My question now is: what is the asymptotic behaviour of $I_n$?

Being unsure about what @Tetis' three (so far) answers achieve exactly, let me post the asymptotic behaviour @Byron's approach yields, pushing things one step further.

Using the notations in @Byron's post, one sees that $I_n$ is $n$ times the expectation of $X_1^2/S_n$, where $S_n=\sum\limits_{k=1}^nX_k$. Define $Z_n$ by the identity $$S_n=nE[X]+\sqrt{n}\sigma(X)Z_n,$$ where $\sigma^2(X)$ denotes the variance of every $X_k$. Then $Z_n$ converges in distribution to a centered normal random variable. Using the expansion $1/(1+t)=1-t+t^2+o(t^2))$ when $t\to0$, one gets $$ n\frac{X_1^2}{S_n}=\frac{X_1^2}{E[X]}\left(1-\frac{\sigma(X)Z_n}{\sqrt{n}E[X]}+\frac{\sigma^2(X)Z_n^2}{nE[X]^2}+o\left(\frac1n\right)\right), $$ hence $$ I_n=nE\left[\frac{X_1^2}{S_n}\right]=\frac1{E[X]}\left(E[X^2]-\frac{\sigma(X)E[X_1^2Z_n]}{\sqrt{n}E[X]}+\frac{\sigma^2(X)E[X_1^2Z_n^2]}{nE[X]^2}+o\left(\frac1n\right)\right). $$ Expanding $Z_n$, one sees that $$ \sigma(X)\sqrt{n}E[X_1^2Z_n]=E[X^3]-E[X^2]E[X], \qquad E[X_1^2Z_n^2]=E[X^2]+o(1), $$ hence, $$ \lim_{n\to\infty}n\cdot\left(I_n-\ell_X\right)=\kappa_X, $$ where $$ \ell_X=\frac{E[X^2]}{E[X]},\qquad\kappa_X=\frac{E[X^2]^2-E[X^3]E[X]}{E[X]^3}. $$ In the case at hand, $E[X^k]=\frac1{k+1}$ for every positive $k$, hence $$ \lim_{n\to\infty}n\cdot\left(I_n-\frac23\right)=-\frac19. $$