Estimating direction from a distribution on a circle - MathOverflow most recent 30 from http://mathoverflow.net 2013-05-26T06:43:09Z http://mathoverflow.net/feeds/question/42139 http://www.creativecommons.org/licenses/by-nc/2.5/rdf http://mathoverflow.net/questions/42139/estimating-direction-from-a-distribution-on-a-circle Estimating direction from a distribution on a circle Andrei Kolin 2010-10-14T12:12:40Z 2011-06-15T11:15:46Z <p>Let there be $n$ points on a unit circle. It is known they come from "normal" distribution around particular unknown direction (i.e. sum of 2 "normal" distributions on circle - one centered at point $p$ and the other at its opposite $-p$). What is the best way to estimate this direction? By best I mean an algorithm that is a. analytical, b. efficient and c. simple.</p> http://mathoverflow.net/questions/42139/estimating-direction-from-a-distribution-on-a-circle/42145#42145 Answer by Robby McKilliam for Estimating direction from a distribution on a circle Robby McKilliam 2010-10-14T12:59:13Z 2010-10-14T12:59:13Z <p>The standard way to solve this is to just consider each of your data points as unit vectors, then take the average of those unit vectors. The direction of this averaged vector is the estimated direction.</p> <p>There is a large literature on this topic which generally goes by the name of <a href="http://en.wikipedia.org/wiki/Directional_statistics" rel="nofollow">directional statistics</a>. The seminal text on is Mardia and Jupp's book <a href="http://www.amazon.com/Directional-Statistics-Kanti-V-Mardia/dp/0471953334" rel="nofollow">Directional Statistics</a>. This field has a huge number of applications in astronomy, biology, meteorology, engineering etc. </p> http://mathoverflow.net/questions/42139/estimating-direction-from-a-distribution-on-a-circle/42229#42229 Answer by Robby McKilliam for Estimating direction from a distribution on a circle Robby McKilliam 2010-10-14T22:59:11Z 2011-06-15T11:15:46Z <p>I see now that Andrei would like to know what to do when the distribution has 2 modes and is <em>symmetric</em> about these modes. It seems better to just give a second (more detailed) answer rather than complicate the simple answer I gave above (basically I think the idea in gowers comment above is sound, but it's a bit tricky to actually implement).</p> <p>So, how do we deal with estimating the 'mean direction' of a distribution that looks something like:</p> <p><img src="http://robbymckilliam.github.com/pictures/mathoverflow/bimod.png" alt="alt text"></p> <p>Good questions at this point are ''what <em>is</em> mean direction anyway?'' and specifically for the distribution above ''does a mean direction even <em>exist</em>?''</p> <p>This has been a question I have been looking at a few months now. I'm wary of blowing my own horn a bit here, but I am going to attach a part of my <a href="http://robbymckilliam.github.com/papers/thesis.pdf" rel="nofollow">thesis</a> which I think gives satisfactory answers to these questions (I would love to give you the whole thesis, but it's not quite ready for the public to see). I suggest that there are (atleast) two different, but equally reasonable and intutive definitions of mean direction. I argue that the distribution above <em>has no mean</em> in a rigorously definable sense for both of these definitions. </p> <p>Given $N$ data points $\Theta_1,\dots, \Theta_N$ on a circle there exist very accurate and efficient O(N)-time algorithms to estimate both of these means <em>if they exist</em>. Neither algorithm will converge if used on circular data drawn from the bimodal distribution above as (according to my definition) the means <em>do not exist</em>. </p> <p>Still, given $N$ data points $\Theta_1,\dots, \Theta_N$ drawn from the bimodal distribution above, if what you want to do estimate one of the ''modes'' rather than the mean direction, then my gut tells me that there probably are efficient and accurate algorithms to do this, although I don't know if they exist in the literature. You could try Fishers book <a href="http://books.google.com.au/books?id=wGPj3EoFdJwC&amp;printsec=frontcover&amp;dq=Fisher+statistical+analysis&amp;source=bl&amp;ots=Ph2uvBACR4&amp;sig=P9t7LCZBqlnD77cxCsxIda2P364&amp;hl=en&amp;ei=YIe3TOCAGYuPcezI5LYG&amp;sa=X&amp;oi=book_result&amp;ct=result&amp;resnum=4&amp;ved=0CCgQ6AEwAw#v=onepage&amp;q&amp;f=false" rel="nofollow">The Statistical Analysis of Circular Data</a>.</p> http://mathoverflow.net/questions/42139/estimating-direction-from-a-distribution-on-a-circle/42240#42240 Answer by Niels Diepeveen for Estimating direction from a distribution on a circle Niels Diepeveen 2010-10-15T00:39:01Z 2010-10-15T00:39:01Z <p>This seems too simple to be true, but, combining some of the ideas posted earlier, I think you could just interpret the vectors as complex numbers and take the RMS.</p> <p>Squaring will turn the bimodal distribution into a unimodal one. Then the square roots of the mean should give a good estimate of the modes of the original distribution.</p> http://mathoverflow.net/questions/42139/estimating-direction-from-a-distribution-on-a-circle/42342#42342 Answer by Robby McKilliam for Estimating direction from a distribution on a circle Robby McKilliam 2010-10-15T23:15:02Z 2011-06-15T11:13:22Z <p>Ok, so now I will describe why Niels's estimator works so well. Take a bimodal and symmetric circular density function $f$ with modes $p$ and $-p$ (we will assume that $p$ is positive) such as the one plotted in my previous answer. Let $\Theta_1, \Theta_2, \dots, \Theta_N$ be $N$ observations drawn from $f$.</p> <p>Niels's estimator first computes the complex numbers $e^{i 2 \Theta_n}$ and takes their average $$\bar{C} = \sum_{n=1}^{N} e^{i 2 \Theta_n} .$$ The estimate, denoted $\hat{p}$, is given by taking the complex argument of $\bar{C}$ and dividing by 2, that is $$\hat{p} = \frac{\angle{\bar{C}}}{2}$$ where $\angle{\bar{C}} \in [0,2\pi)$ denotes the complex argument. The next theorem describes the asymptotic properties of this estimator. I use the notation $\langle x \rangle_{\pi}$ do denote $x$ taken to its representative inside $[-\pi, \pi)$. So, for example, $\langle 2\pi \rangle_{\pi} = 0$ and $\langle \pi + 0.1 \rangle_{\pi} = -\pi + 0.1$.</p> <blockquote> <p><strong>Theorem:</strong> Let $\lambda$ denote the difference $\lambda = \tfrac{1}{2}\langle 2\hat{p} - 2p \rangle_{\pi}.$ Then $\lambda$ converges almost surely to zero as $N \rightarrow \infty$ and the distribution of the normalised difference $\sqrt{N}\lambda$ converges to the zero mean normal with variance $$\frac{\sigma_s^2}{c}$$ where $$\sigma_s^2 = \int_{-\pi/2}^{\pi/2}\sin^2(\theta) f(\langle \theta + p \rangle_\pi) d\theta \qquad \text{and} \qquad c = \int_{-\pi/2}^{\pi/2}\cos(\theta) f(\langle \theta + p \rangle_\pi) d\theta.$$</p> </blockquote> <p>The definition of the difference $\lambda$ might seem a little strange at first, but it is actually very natural. To see why note that $p$ and the estimate $\hat{p}$ are both in $[0,\pi)$ but, for example, if $p = 0$ and $\hat{p} = \pi - 0.01$ then the difference between these is <em>not</em> $\pi - 0.01$, because the two modes are actually very close to aligned in this case. The correct difference is $\lambda = \tfrac{1}{2}\langle 2(\pi-0.01) - 2 \times 0 \rangle_{\pi} = 0.01$. </p> <p>The proof of this theorem follows from a very similar argument to Theorem 6.1 (page 87) from <a href="http://robbymckilliam.github.com/papers/thesis.pdf" rel="nofollow">my thesis</a>. The original argument is due to <a href="http://www.stat.mq.edu.au/our_staff/staff_-_alphabetical/staff1/barry_quinn/" rel="nofollow">Barry Quinn</a>. Rather than restate the proof I'll just give you some convincing numerical evidence.</p> <p>I've run some simulations for the case when the noise is a sum of two weighted von Mises circular distributions with <em>concentration parameter</em> $\kappa$. So, when $\kappa$ is large the distribution is concetrated and looks something like the picture on the left below ($\kappa = 20$ in this case) and when $\kappa$ is small the distribution is quite spread out and looks something like the picture on the right below ($\kappa = 0.5$). We obviously expect the estimator to perform better when the distribution is quite concentrated ($\kappa$ is large). </p> <p><img src="http://robbymckilliam.github.com/pictures/mathoverflow/bimodthin.png" alt="alt text"> <img src="http://robbymckilliam.github.com/pictures/mathoverflow/bimodthick.png" alt="alt text"></p> <p>Here are the results. The plot below show the simulated variance of $\lambda$ after 5000 trials (the dots) versus the variance predicted in the theorem above for a range of values of $\kappa$ and number of observations $N$. You can see that the theorem does a very good job of accurately predicting the perfomance if $\kappa$ isn't too small.</p> <p><img src="http://robbymckilliam.github.com/pictures/mathoverflow/direst.png" alt="alt text"></p> <p>There is still an open question as to whether this is the <em>best</em> estimator (in the sense of maximally reducing the variance of $\lambda$). It would be possible to derive a <a href="http://en.wikipedia.org/wiki/Cram%25C3%25A9r%25E2%2580%2593Rao_bound" rel="nofollow">Cramer-Rao bound</a> for this estimation problem to give an idea of the best possible performance of an unbiased estimator. I suspect that this estimator performs very near the Cramer-Rao bound. So, in that sense it is close to best possible.</p>