<s>[INST] Derivation of integrating over many parameters in Neyman-Scott Problem? I am trying to follow the derivation for the variance estimator in the Neyman-Scott problem given in this article. However, I'm not sure how they go from the 2nd to the 3rd line of this derivation. Any help is appreciated, thanks!

 [/INST] Each of the integrals indexed by $i$ in the product has the form
$$\int_{\mathbb R} \exp\left[\phi(\mu_i,x_i,y_i,\sigma)\right]\,\mathrm{d}\mu_i.$$
When $\phi$ is linear or quadratic in $\mu_i,$ such integrals have elementary values.  The more difficult circumstance is the quadratic.  We can succeed with such integrations merely by knowing
$$\int_\mathbb{R} \exp\left[-\frac{1}{2}z^2\right]\,\mathrm{d}z = \sqrt{2\pi}.$$
Generally, let $\theta$ denote everything that isn't $\mu_i$ ($\theta=(x_i,y_i,\sigma)$ in this instance).  Writing $x=\mu_i,$ the quadratic case is when $\phi$ can be expressed in the form
$$\phi(x,\theta) = -A(\theta)^2 x^2 + B(\theta)x + C(\theta)$$
for arbitrary functions $A, B, C$ and $A(\theta)\ne 0.$  The reasons for expressing the coefficient of $x^2$ in this way are (i) to guarantee the integral exists and (ii) to avoid using square roots.
Thus, we are concerned with evaluating
$$f(\theta) = \int_\mathbb{R} \exp\left[\phi(x,\theta)\right]\,\mathrm{d}x=\int_\mathbb{R} \exp\left[-A(\theta)^2x^2 + B(\theta)x + C(\theta)\right]\,\mathrm{d}x.$$
This done by completing the square, in exactly the same way the quadratic formula is traditionally derived.  The result amounts to changing the variable of integration to $z$ where
$$z = A\sqrt{2}\,x - \frac{B}{A\sqrt{2}};\quad \mathrm{d}z = A\sqrt{2}\,\mathrm{d}x.$$
In terms of $z,$
$$\exp\left[\phi(x,\theta)\right] = \exp\left[-\frac{1}{2}z^2 + C(\theta) + \frac{B(\theta)^2}{4A(\theta)^2}\right]= \exp\left[-\frac{1}{2}z^2\right] \exp\left[C(\theta) + \frac{B(\theta)^2}{4A(\theta)^2}\right].$$
The integral thereby becomes
$$\begin{aligned}
f(\theta) &= \frac{1}{A\sqrt{2}} \int_\mathbb{R} \exp\left[-\frac{1}{2}z^2\right] \exp\left[C(\theta) + \frac{B(\theta)^2}{4A(\theta)^2}\right]\,\mathrm{d}z \\
&= \frac{1}{A\sqrt 2} \sqrt{2\pi}\exp\left[C(\theta) + \frac{B(\theta)^2}{4A(\theta)^2}\right] \\
&=\frac{\sqrt \pi}{A} \exp\left[C(\theta) + \frac{B(\theta)^2}{4A(\theta)^2}\right].
\end{aligned}$$
Let's memorialize this for future reference:
$$\boxed{\int_\mathbb{R} \exp\left[-A(\theta)^2x^2 + B(\theta)x + C(\theta)\right]\,\mathrm{d}x = \frac{\sqrt \pi}{A} \exp\left[C(\theta) + \frac{B(\theta)^2}{4A(\theta)^2}\right].}\tag{*}$$

To apply this to the derivation in the question, simply look at the argument of the exponential in the integrals and break it down into the form of $\phi;$ namely, as a linear combination of $\mu_i^2,$ $\mu_i,$ and a constant:
$$-\frac{(x_i-\mu_i)^2 + (y_i-\mu_i)^2}{2\sigma^2} = -\frac{1}{\sigma^2}\mu_i^2 + \frac{x_i+y_i}{\sigma^2} \mu_i + -\frac{x_i^2 + y_i^2}{2\sigma^2},$$
from which we read off
$$\cases{A = \frac{1}{\sigma} \\ B = \frac{x_i+y_i}{\sigma^2} \\ C =  -\frac{x_i^2 + y_i^2}{2\sigma^2},}$$
whence
$$C(\theta) + \frac{B(\theta)^2}{4A(\theta)^2} =  -\frac{x_i^2 + y_i^2}{2\sigma^2} + \left(\frac{x_i+y_i}{\sigma^2}\right)^2 \frac{\sigma^2}{4} = -\frac{(x_i-y_i)^2}{4\sigma^2}.$$
Plugging everything into the formula $(*)$ gives--also by visual inspection--
$$f(x_i,y_i,\sigma) = f(\theta) = \sigma\sqrt\pi \exp\left[ -\frac{(x_i-y_i)^2}{4\sigma^2}\right].$$
This takes us from the second line to the third line in the derivation.
</s>