Variance computed using Taylor series does not agree with numerical experiment

Question

I would like to estimate an angle $\theta\in\left(-\frac{\pi}{2},\frac{\pi}{2}\right)$ given the noisy observations of its sine and cosine (this is related to my earlier question). My estimator is the inverse tangent of the ratio of the means of the observations of sine and cosine. Let's assume that the noise is additive and Gaussian. I am having trouble showing that the mean square error (MSE) of this estimator is decreasing in $n$, though numerical experiments (and intuition) suggest that it does.

Formally, suppose independent sequences of i.i.d. zero-mean Gaussian random variables $(X_n)$ and $(Y_n)$ have variance that is inversely proportional to $n$. That is $X_n\sim\mathcal{N}\left(0,\frac{\sigma^2}{n}\right)$, $Y_n\sim\mathcal{N}\left(0,\frac{\sigma^2}{n}\right)$, and the joint p.d.f. is:

$$f_{X_n,Y_n}(x_n,y_n)=\frac{n}{2\pi\sigma^2}\exp\left[-\frac{(x_n^2+y_n^2)n}{2\sigma^2}\right].$$

I am interested in the mean squared error of the estimator $\hat{\theta}(X_n,Y_n)=\tan^{-1}\left(\frac{\sin(\theta)+X_n}{\cos(\theta)+Y_n}\right)$ as $n\rightarrow\infty$. Specifically, I would like to show that:

$$\lim_{n\rightarrow\infty}\mathbb{E}_{X_n,Y_n}[(\hat{\theta}(X_n,Y_n)-\theta)^2]=\lim_{n\rightarrow\infty}\mathbb{E}_{X_n,Y_n}\left[\left(\tan^{-1}\left(\frac{\sin(\theta)+X_n}{\cos(\theta)+Y_n}\right)-\theta\right)^2\right]=0.$$

Furthermore, I would like to show that $\mathbb{E}_{X_n,Y_n}[(\hat{\theta}(X_n,Y_n)-\theta)^2]=\mathcal{O}\left(\frac{1}{n}\right)$.

I performed a numerical experiment in MATLAB (see code below), which suggests the inverse linear scaling of the MSE:

I then tried to confirm it analytically using the Taylor series expansion of the function $\tan^{-1}\left(\frac{\sin(\theta)+X_n}{\cos(\theta)+Y_n}\right)$ at the origin $(x,y)=(0,0)$. However, utmost care in the control of the remainder must be exercised when taking the expectation of Taylor expansion--see this discussion (in particular, Mike McCoy's answer). With that in mind, I use Wolfram Mathematica (see code below) I obtain:

$$\begin{align}&\tan^{-1}\left(\frac{\sin(\theta)+x}{\cos(\theta)+y}\right) =\theta+x \cos (\theta )-y \sin (\theta )\\ &\qquad+\frac{1}{2} \left[\sin (2 \theta ) \left(y^2-x^2\right)-2 x y \cos (2 \theta )\right] \\ &\qquad+\frac{1}{3} \left[-y \sin (3 \theta ) \left(y^2-3 x^2\right)-x \cos (3 \theta ) \left(x^2-3 y^2\right)\right] \\ &\qquad+\frac{1}{4} \left[\sin (4 \theta ) \left(x^4-6 x^2 y^2+y^4\right)+4 x y \cos (4 \theta ) \left(x^2-y^2\right)\right]\\ &\qquad+\frac{1}{5} \left[-y \sin (5 \theta ) \left(5 x^4-10 x^2 y^2+y^4\right)+x \cos (5 \theta ) \left(x^4-10 x^2 y^2+5 y^4\right)\right]\\ &\qquad+\frac{1}{6} \left[-\sin (6 \theta ) \left(x^6-15 x^4 y^2+15 x^2 y^4-y^6\right)\right. \\ &\qquad\qquad\qquad\left.+\cos (6 \theta ) \left(-\left(6 x^5 y-20 x^3 y^3+6 x y^5\right)\right)\right] \\ &\qquad+\frac{1}{7}\left[-y \sin (7 \theta ) \left(-7 x^6+35 x^4 y^2-21 x^2 y^4+y^6\right)\right. \\ &\qquad\qquad\qquad\left.-x \cos(7 \theta ) \left(x^6-21 x^4 y^2+35 x^2 y^4-7 y^6\right)\right]\\ &\qquad+\ldots\end{align}$$

The zeroth-order term is of course $\theta$. If we just plug $X_n$ and $Y_n$ into the first order term $x \cos (\theta )-y \sin (\theta )$ and compute the variance, we get the desired inverse dependence on $n$. However, per Mike McCoy's answer to the aforementioned post, I need to ensure that the error from the remainder $R_2(X_n,Y_n;\theta)$ after the first order term goes to zero as well, where the remainder is defined as follows:

$$R_2(x,y)=\tan^{-1}\left(\frac{\sin(\theta)+x}{\cos(\theta)+y}\right)-\theta-x \cos (\theta )+y \sin (\theta ).$$

Remainder is an arithmetic mess, however, subtracting $\theta$ from the Taylor series expansion above, squaring, plugging in $X_n$ and $Y_n$, and finding expected value by integrating in Mathematica (see code below) yields:

$$\mathbb{E}_{X_n,Y_n}[(\hat{\theta}(X_n,Y_n)-\theta)^2]=\left(\frac{\sigma^2}{n}\right)+\left(\frac{\sigma^2}{n}\right)^2+\frac{8}{3} \left(\frac{\sigma^2}{n}\right)^3+12 \left(\frac{\sigma^2}{n}\right)^4+\frac{384}{5}\left(\frac{\sigma^2}{n}\right)^5+640 \left(\frac{\sigma^2}{n}\right)^6+\frac{46080}{7} \left(\frac{\sigma^2}{n}\right)^7+80640 \left(\frac{\sigma^2}{n}\right)^8+\ldots$$

From these first few terms one recognizes the following series:

$$\mathbb{E}_{X_n,Y_n}[(\hat{\theta}(X_n,Y_n)-\theta)^2]=\sum_{p=1}^\infty\frac{2^pp!\sigma^{2p}}{pn^p},$$

which diverges for any fixed $n$ (and $\sigma^2$), implying that the MSE is infinite.

What am I doing wrong? Is there any other way to show that the MSE converges to zero?

---

CODE

This is the MATLAB code used in generating the figure:

for n=1:100
  x=randn(1,1000)/sqrt(n*10);
  y=randn(1,1000)/sqrt(n*10);
  v(n)=var(atan((sin(theta)+x)./(cos(theta)+y)));
end

This is the Mathematica code to obtain the first seven terms of the Taylor series expansion:

Sum[FullSimplify[
 D[ArcTan[(Sin[\[Theta]] + t*x)/(Cos[\[Theta]] + t*y)], {t, i}] /. 
 t -> 0, Assumptions -> \[Theta] > -Pi/2 && \[Theta] < 
  Pi/2]/(i!), {i, 0, 7}]

and here is the code used to first the MSE to the first few terms:

Integrate[(Out[2] - \[Theta])^2/(2*Pi*s2)*
  Exp[-(x^2 + y^2)/(2*s2)], {x, -Infinity, Infinity}, {y, -Infinity, 
  Infinity}, 
 Assumptions -> \[Theta] > -Pi/2 && \[Theta] < Pi/2 && s2 > 0]

Answer 1

It does not matter that whether the mean squared error is decreasing with $n$. All that matters is that it can be made arbitrarily small. Here is a simple demonstration.

Fixing $\theta$, let $s=\sin(\theta)$ and $c=\cos(\theta)$. The estimator $\hat\theta$ is the slope of the ray through the origin and the point

$$(c + Y_n, s+X_n).$$

Let $\epsilon \gt 0$. Consider the disk of radius $\sin\delta$ around $(c,s)$. Euclidean geometry shows us this disk is contained within the wedge at the origin with angles $\theta-\delta$ to $\theta+\delta$. By choosing $n$ sufficiently large, the chance that $(c+Y_n, s+X_n)$ lies within that disk can be made to exceed $1-\alpha$ for any tiny $\alpha$. In fact, because $n(X_n^2 + Y_n^2)$ follows a $\chi^2_2$ distribution,

$$n = \lceil\frac{(\chi^2_2)^{-1}(1-\alpha)}{(\sin\delta)^2}\rceil$$

will work.

In this sketch, $(c,s)$ is at the red dot, the disk is drawn in black, the wedge in blue, and 10,000 simulated values of $(c+Y_n, s+X_n)$ are shown. $\epsilon$ was equal to $1/50$ here, corresponding to a root mean error in the slope estimate of no more than $\sqrt{1/50}\approx 0.14$. With $\alpha=0.05$, the value of $n$ was $300$, corresponding to a standard deviation of $\sqrt{1/300} \approx 0.058$ in the individual coordinates of the simulated points.

Now, within that disk the angular error does not exceed $\delta$ and outside that disk the angular error cannot exceed $\pi$ under any circumstance (since the inverse tangent always produces values in the range $[-\pi/2, \pi/2]$). This bounds the expected squared error:

$$\mathbb{E}([\hat\theta - \theta]^2) \le (1-\alpha)(\sin(\delta))^2 + \alpha \pi^2.$$

By choosing, say, $\delta = \arcsin{\sqrt{\epsilon}}$ and $\alpha=\epsilon/\pi^2$, the right hand side will be less than $2\epsilon$. Because $\epsilon \gt 0$ was arbitrary, the limiting mean square error is zero.

(In the figure, the mean square error was $0.0033 \ll 0.04 = 2\epsilon$.)

---

The argument concerning angular error in the disk assumed the inverse tangent was continuous within a neighborhood of $\theta$. That will not be the case for $\theta=\pm \pi/2$, but a slight change in the definition of the inverse tangent in cases where $\theta$ is close to these values will fix the problem.