Neyman–Scott problems, by hand

Companion to the notes on Chapter 9 of Amari's Information Geometry and Its Applications. Back to the notes. A model $p(x;u,v)$ has a parameter of interest $u$ and a nuisance parameter $v$; in a Neyman–Scott problem every observation has its own $v_i$, so the nuisance parameter grows with the data and the maximum likelihood estimate can settle on the wrong value. The cure of the chapter is an estimating function $f(x,u)$ whose expectation vanishes for every $v$. All widgets use the book's examples (the coefficient of proportionality, the variance of a scale, the gamma shape of a spike train) or the exact finite model of the notes; the numbers match code/neyman_scott.py up to Monte Carlo error, and simulations use a seeded generator.

1. What an unknown nuisance parameter costs

$N$ pairs $(x,y)$ share one specimen weight $v$: $x\sim N(v,1)$, $y\sim N(uv,1)$. The scores $e_u$ and $e_v$ are two vectors in a plane with the Fisher inner product $g$; knowing $v$ you could use all of $e_u$, but when $v$ is unknown only the part orthogonal to $e_v$ is usable, the efficient score $\bar e_u=e_u-g_{uv}g_{vv}^{-1}e_v$ (9.14), whose squared length is $\bar g=g_{uu}-g_{uv}^2/g_{vv}$ (9.10). Here $\bar g=v^2/(1+u^2)$: the variance of $\hat u$ is multiplied by $1+u^2$. Move $u$ and watch the angle open (the loss is $\sin^2$ of it); the button simulates the MLE $\hat u=\bar y/\bar x$ at $N=2000$. The curves $v'=v\sqrt{1+u^2}=\text{const}$ make an orthogonal nuisance parameter: in $(u,v')$ the loss vanishes.

2. No orthogonal nuisance parameter for a vector $u$

For a scalar $u$ the condition $g_{u\lambda}=0$ is an ordinary differential equation $dv/du=-g_{uv}/g_{vv}$, always solvable. For two parameters of interest it asks that the "horizontal" directions (orthogonal to the nuisance direction) be level sets of a function $v'$, which needs a loop in the $u$-plane to lift to a closed curve. Take $x\sim N(\mu,I_3)$ with $\mu=(u_1,u_2,0)+v\,d(u)$, $d=(-u_2,u_1,1)/\sqrt{1+|u|^2}$. Lifting the circle $|u|=\rho$ horizontally, $v$ comes back changed by $-2\pi\rho^2/\sqrt{1+\rho^2}$: the nuisance parameter that would make $g_{a\lambda}=0$ everywhere does not exist (the remark after (9.19)). Orthogonality at one point is not obstructed.

3. The MLE settles on the wrong value

Two Neyman–Scott problems with one nuisance parameter per box. Gaussian variance: $N$ boxes of $m$ measurements, each box with its own mean; the joint MLE of $\sigma^2$ (here $\sigma^2=2$) is $(m-1)/m$ times the within-box estimator (9.114), which is unbiased. Matched pairs: $N$ pairs of binary outcomes with a free log-odds per pair and odds ratio $u=2$ between the two members; the joint MLE of $u$ is $(n_{01}/n_{10})^2\to u^2$, the ratio of discordant counts $n_{01}/n_{10}\to u$. Add boxes and the histograms narrow, but around the wrong centre. (The sums of squares are drawn from their exact laws; 300 data sets per histogram.)

4. Five solutions of the proportionality problem

Each specimen $i$ has its own weight $v_i$ (drawn here from $N(\mu_v,\tau^2)$, the semiparametric reading of the book): $x_i=v_i+\varepsilon_i$, $y_i=uv_i+\varepsilon_i'$ with unit noise. Fitted slopes through the origin: least squares (9.23) shrinks towards $0$ (limit $u\,\mathbb E v^2/(\mathbb E v^2+1)$); the average of ratios (9.24) is wild, since $y/x$ has no mean; the gross average (9.25) and total least squares (9.26), which is also the MLE, are consistent. TLS has a second root $-1/\hat u$, the line orthogonal to the fit. Try a small $\mu_v$ (mean zero is where the gross average dies) and increase $n$.

5. The averaging method has no law of large numbers

$\mathbb E|y/x|=\infty$ because $x$ has a positive density at $0$. The average of $N$ ratios does not converge in probability: it converges in law to a Cauchy distribution with scale $\pi a$, $a=\mathbb E_v[\varphi(v)\,\mathbb E|y|]$, centred at the principal-value mean $u\,\mathbb E_v[v\sqrt2 F(v/\sqrt2)]$ ($F$ is Dawson's integral), which is not $u$. Six runs of the running average (orange) and one of the gross average (green); the histogram is 400 independent averages of $N=1000$ ratios against the predicted Cauchy curve. Shrink $\tau$ and raise $\mu_v$: $a$ becomes exponentially small and the effect hides until $N$ is astronomical.

6. The class (9.97) and the formula (9.103)

Every estimating function $f=(y-ux)\,h(x+uy)$ gives a consistent estimator. For $h(s)=s+c$ ($c=0$ is TLS, $c\to\infty$ the gross average) the asymptotic variance per specimen is $(1+u^2)\{(c+(1+u^2)\bar v)^2+(1+u^2)^2\operatorname{Var}v+(1+u^2)\}/(c\bar v+(1+u^2)\overline{v^2})^2$ with the whole brace multiplied by $1+u^2$; (9.103) as printed (with its evidently missing $1/N$ and the square of the denominator supplied) multiplies only the first term, which makes the curve (dashed) dip below the efficiency bound for Gaussian $k$. The minimum is at $c=\bar v/\operatorname{Var}v$ (9.105), and for Gaussian $k$ it equals the bound $q/\mathbb E[(\mathbb E[v|s])^2]$, $q=1+u^2$. Move the spread $\tau$: below the crossover TLS loses to the gross average, above it wins.

7. The geometry in an exact finite model

$x,y\in\{0,\dots,m\}$, odds ratio $u$ between $y$ and $x$, nuisance $v$ = log-odds of $x$ (the natural parameter of $s=x+y$). The mixing law $k$ lives on the nine-point grid $v=-2,-1.5,\dots,2$. The nuisance tangent space is $\{h(s)\}$, the efficient score is $\dot l_E=(y-E[y|s])/u$ with $E[y|s]$ from a non-central hypergeometric law, and it is the same for every $k$, while the efficient information $\bar g=E\dot l_E^2$ depends on $k$ and sits below the information with $k$ known, which sits below the information with $v$ known. Any other estimating function is $\dot l_E+a$, $a\perp\dot l_E,T_K$: its variance is $(1+t^2)/\bar g$ when $\|a\|^2=\bar g$ and $a$ enters with weight $t$ (Theorem 9.6).

8. A spike train: four estimators of the gamma shape

Intervals are gamma with shape $\kappa$ (regularity) and a rate that changes from box to box (two intervals per box, $n$ boxes, rate $e^{\sigma_v Z}$). The efficient estimating equation (9.131) sets $S=-\tfrac12\overline{\log(4D(1-D))}$, $D=T_1/(T_1+T_2)$, equal to $\phi(2\kappa)-\phi(\kappa)-\log2$; $L_V$ (9.133) is consistent and less efficient; the joint MLE (a rate per box) solves the wrong equation $S=\log\kappa-\phi(\kappa)$ and overshoots by nearly a factor 2; the naive estimate pools all intervals and reads the rate changes as irregularity. $D$ does not depend on the rate at all. Raise $\varepsilon$, the fraction of artefact boxes with $D=10^{-6}$, to see why $L_V$ may be preferred.