# Variance computed using Taylor series does not agree with numerical experiment

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]

• Your messy remainder starts out as an expression that depends on $\theta$. When you integrate that dependency out, you are including the area near $\theta=\frac{\pi}{2}$. In that region, the remainder term grows rapidly as $\theta \to \frac{\pi}{2}$ so you would not expect the remainder to go to zero. If your region of integration were restricted to eliminate this, then the power series yo obtained would be messier but would assumedly converge. – Mark Fischler Jun 24 '16 at 21:50
• Are you sure the single quadrant atan(z) and ArcTan[z] are adequate as opposed to atan2 and ArcTan[x,y]? – horchler Jun 24 '16 at 22:57
• @Mark Could you please clarify your comment? I'm not integrating over $\theta$. I think sines and cosines get squared and cancel out when one calculates variance. Then again, in my numerical experiments, as $\theta$ gets closer to $\pm\pi/2$, the $n$ required to converge to small MSE gets larger. So, perhaps there is a problem there... However, I am not sure how to exclude $\theta$ from the regions around $\pm\pi/2$... – M.B.M. Jun 25 '16 at 1:33
• @horchler Using ArcTan[x,y] yields the same result... – M.B.M. Jun 25 '16 at 4:33
• The Taylor series approach has little hope of succeeding, because the series converges only for arguments in absolute value less than $1$, whereas adding the Normal error assures there is positive probability of arguments greater than this. – whuber Jun 25 '16 at 18:31

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.