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!


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.
