Report_img.ipynb

<p><strong>Von Mises Distribution Implementation in Python</strong></p>
<p>Vedant Mehta, Kanchan Satpute</p>
<p><em>Texas A&amp;M University, College Station, Texas</em></p>
<blockquote>
<p><strong>Abstract:</strong> This article tries to bridge the gap between the quality of package available for R and Python for circular statistics and von Mises distribution in specific. There is a function available in Python to generate random deviates from vonMises distribution. But there are no functions available to calculate the probability density, cumulative distribution, quantiles, etc.</p>
</blockquote>
<h1 id="introduction">1. Introduction</h1>
<p>Directional statistics or circular statistics is a sub-discipline of statistics that deals with directions, axes and rotation. Think of it as a regular linear data converted into a circular data by giving it attributes like rotation, angle, etc. Circular statistics is a lot different than linear statistics. Firstly, there is no true zero. Namely, 0 and 360 degrees are equal. So, labeling a value as high or low is arbitrary. Due to these characteristics, method of analysis of this kind of data changes completely. The kind of data that has angles, or periodicity, or does not have a true zero can be labeled as directional data. Some of the examples include temporal periods (e.g. time of day, month, hour, week, etc.), compass directions, daily wind directions, ocean currents, etc.<sup>[5]</sup></p>
<p>Calculation of mean, median and variance of a circular data is quite different from that in linear statistics. If given a data of angles, it cannot be simply averaged like it is done in linear statistics.</p>
<p>Method for mean calculation:</p>
<p><strong>Example:</strong> Given angular data <span class="math inline">(<em>α</em><sub>1</sub>,  <em>α</em><sub>2</sub>, ….,<em>α</em><sub><em>n</em></sub>)</span>. Calculate the sine and cosine of all the angles.</p>
<p>Further, <span class="math inline">$X = \ \frac{\sum_{i = 1}^{n}{\cos\alpha_{i}}}{N}$</span> and <span class="math inline">$Y = \frac{\sum_{i = 1}^{n}{\sin\alpha_{i}}}{N}$</span>. Also, <span class="math inline">$\overset{\overline{}}{r} = \sqrt{X^{2} + Y^{2}}$</span>. So, mean cosine will be</p>
<p><span class="math inline">$\cos\overset{\overline{}}{\alpha} = \frac{X}{r}$</span> and mean sine will be <span class="math inline">$\sin\overset{\overline{}}{\alpha} = \frac{Y}{r}$</span>. Finally, mean angle will be</p>
<p><br /><span class="math display">$$\theta_{r} = \arctan\left( \frac{\sin\overset{\overline{}}{\alpha}}{\cos\overset{\overline{}}{\alpha}} \right)$$</span><br /></p>
<p>Decide the resultant quadrant in following way: (Figure given for reference)</p>
<p>Sin +, cos + : mean angle computed directly</p>
<p>Sin +, cos - : mean angle = <span class="math inline">180  −  <em>θ</em><sub><em>r</em></sub></span></p>
<p>Sin -, cos - : mean angle = <span class="math inline">180 + <em>θ</em><sub><em>r</em></sub></span></p>
<p>Sin -, cos +: mean angle = <span class="math inline">360 − <em>θ</em><sub><em>r</em></sub></span></p>
<p>Circular variance measures variation in the angles about the mean direction.</p>
<p>Variance <span class="math inline">$V = 1 - \overset{\overline{}}{r}$</span>. So it ranges from 0-1. When the variance is 1, it means the vectors are concentrated in one direction. Value of 0 means the vectors are equally dispersed around the circle.</p>
<p>There is another kind of data known as <em>bimodal data.</em> When data have opposite angles they are said to have diametrically bimodal circular distributions. The mean angle of diametrically bimodal data is orthogonal (at right angle) to the true mean. There is a procedure called <em>angle doubling</em> to deal with the diametrically bimodal data. But this article won’t be discussing on that topic.</p>
<p>There are different types of distributions defined. Generally speaking, any kind of probability density function can be wrapped around the circumference of a circle. Von Mises distribution is one of the circular distributions that are defined in circular statistics and can be considered as analogous to normal distribution in linear statistics. Also, it is a close approximation to “wrapped normal” distribution.</p>
<p>Probability density function is given by:</p>
<p><br /><span class="math display">$$f\left( x \middle| \mu,\kappa \right) = \frac{1}{2\pi I_{o}(\kappa)}\exp\left\lbrack \kappa\cos{(x - \mu)} \right\rbrack$$</span><br /></p>
<p>where, <span class="math inline"><em>I</em><sub><em>o</em></sub>(<em>κ</em>)</span> – Modified Bessel function of zero<sup>th</sup> order.</p>
<p><span class="math inline"><em>μ</em></span> – measure of location (similar to mean in Normal distribution)</p>
<blockquote>
<p><span class="math inline"><em>κ</em></span> – measure of concentration ( <span class="math inline">1/<em>κ</em></span> is analogous to <span class="math inline"><em>σ</em><sup>2</sup></span> )</p>
</blockquote>
<p><span class="math inline"><em>I</em><sub><em>o</em></sub>(<em>κ</em>)</span> is defined as:</p>
<p><br /><span class="math display">$$I_{0}\left( x \right) = \sum_{r = 0}^{\infty}\frac{1}{r!r!}\left( \frac{x}{2} \right)^{2r}$$</span><br /></p>
<p><br /><span class="math display">$$\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \  = \frac{1}{2\pi}\int_{0}^{2\pi}e^{x\cos\theta}\text{dθ}$$</span><br /></p>
<p>The above equation is the zero<sup>th</sup> order modified Bessel function.</p>
<p>When <span class="math inline"><em>κ</em></span> = 0, Von Mises distribution reduces to the uniform distribution. As <span class="math inline"><em>κ</em></span> increases, the von Mises distribution approaches normal distribution.</p>
<p><img src="./Figures/media/image1.png" width="202" height="151" /> <img src="./Figures/media/image2.png" width="201" height="151" /></p>
<p>Left panel depicts the PDF and right panel depicts the CDF of the von Mises distribution.</p>
<p>In the right panel, as <span class="math inline"><em>κ</em></span> increases, the S curve will gradually become a straight line. That makes sense in a way, that as the concentration increases, the probability is more cumulated near the mean, i.e. zero in our case.</p>
<h1 id="methods">2. Methods</h1>
<p>We have used the vonMises function in the circular package of R as a reference for generating the algorithms for each method.</p>
<p><strong>rvonmises(n, mu, kappa)</strong></p>
<p>Description – A method for generating random numbers for a von Mises circular distribution.</p>
<p>Arguments –</p>
<p>n – number of observations</p>
<p>Examples -</p>
<p><strong>dvonmises(x, mu, kappa) </strong></p>
<p>Description – A method for calculating the probability density at the given points for a von Mises circular distribution.</p>
<p>Arguments –</p>
<p>x – A vector containing the points at which the density is to be calculated. The object is from class ‘circular’</p>
<p>log – logical; if True, probabilities p is given as log(p). The default value for log is given as False.</p>
<p>Examples –</p>
<p><strong>pvonmises(q, mu, kappa)</strong></p>
<p>Description – Method used to calculate the cumulative distribution at the given points for a von Mises distribution.</p>
<p>Arguments –</p>
<p>q – A vector containing the points at which the distribution is to be calculated. The object is from class ‘circular’</p>
<p>tol – the precision in evaluating the distribution function. Default value = 1e-20</p>
<p>Examples –</p>
<p><strong>qvonmises(p, mu, kappa)</strong></p>
<p>Description – A method used to calculate the quantiles for the given probabilities for a von Mises distribution.</p>
<p>Arguments –</p>
<p>p – A vector containing the probabilities at points at which the quantiles are to be calculated. The object is from class ‘circular’</p>
<p>from_ - a value used for evaluating pvonmises and qvonmises. Default = None</p>
<p>tol – machine epsilon value raised to 0.6</p>
<p>Examples –</p>
<p>Common arguments for all the methods:</p>
<p>mu – The mean direction of the distribution. This object is from class ‘circular’</p>
<p>kappa – non-negative value for the concentration of the distribution</p>
<h1 id="results-and-discussion">3. Results and Discussion</h1>
<p>We run the functions pvonmises, qvonmises, dvonmises with various values of parameters mu and kappa. Below shown is the table that shows the comparison of the values obtained in R and values obtained by the package we built in Python.</p>
<table>
<thead>
<tr class="header">
<th>Method</th>
<th>R</th>
<th>Python</th>
</tr>
</thead>
<tbody>
<tr class="odd">
<td><p>pvonmises(2, 1, 6)</p>
<p>pvonmises([2, 0.8], 2, 6)</p></td>
<td><p>[0.9888944]</p>
<p>[0.5 , 0.003595458]</p></td>
<td><p>[0.988894]</p>
<p>[0.5 , 0.00359546]</p></td>
</tr>
<tr class="even">
<td><p>dvonmises(0.5, 1, 6)</p>
<p>dvonmises([1, 3], 3, 6)</p></td>
<td><p>[0.4581463]</p>
<p>[1.949157e-04, 9.54982e-01]</p></td>
<td><p>[0.45814625]</p>
<p>[1.949157e-04, 9.54982e-01]</p></td>
</tr>
<tr class="odd">
<td><p>qvonmises(0.5, 1, 6)</p>
<p>qvonmises([0.2, 0.6], 2, 7)</p></td>
<td><p>[1]</p>
<p>[1.67413597, 2.09767203]</p></td>
<td><p>[1]</p>
<p>[1.67413597, 2.09767203]</p></td>
</tr>
</tbody>
</table>
<p>Now, we will plot some graphs to demonstrate how precise our values are when compared to those in R</p>
<p>When we run the function rvonmises(n=1000, mu=1, kappa=1), it generates following output in R and Python respectively.</p>
<p><img src="./Figures/media/image3.png" width="234" height="211" /> <img src="./Figures/media/image4.png" width="242" height="215" /></p>
<p>Figure 1: rvonmises in R (left panel) and Python (right panel)</p>
<p>When we run the function dvonmises(x = np.linspace(-pi, pi, 1000), mu=1, kappa=6), it generates following output in R and Python respectively.</p>
<p><img src="./Figures/media/image5.png" width="220" height="177" /> <img src="./Figures/media/image6.png" width="229" height="172" /></p>
<p>Figure 2: dvonmises in R (left panel) and Python (right panel)</p>
<h1 id="future-scope">4. Future Scope</h1>
<p>We need to make the package more robust so that the function can accept different kind of inputs. When we ran the benchmarking tests, we saw that our code took longer time to execute as compared to that in R. So we need to optimize the code in order to decrease the execution time. We can include other functions from the ‘circular’ package of R into Python.</p>
<h1 id="reference">5. Reference</h1>
<p>[1] <a href="https://www.researchgate.net/figure/Wind-data-for-KRDM-the-nearest-FAA-weather-reporting-station-at-the-Redomond-OR_fig5_261417337" class="uri">https://www.researchgate.net/figure/Wind-data-for-KRDM-the-nearest-FAA-weather-reporting-station-at-the-Redomond-OR_fig5_261417337</a></p>
<p>[2] <a href="https://ncss-wpengine.netdna-ssl.com/wp-content/uploads/2013/01/Rose-Plot.png" class="uri">https://ncss-wpengine.netdna-ssl.com/wp-content/uploads/2013/01/Rose-Plot.png</a></p>
<p>[3] <a href="http://webspace.ship.edu/pgmarr/geo441/lectures/lec%2016%20-%20directional%20statistics.pdf">http://webspace.ship.edu/pgmarr/geo441/lectures/lec%2016%20-%20directional%20statistics.pdf</a></p>
<p>[4] <a href="https://ncss-wpengine.netdna-ssl.com/wp-content/themes/ncss/pdf/Procedures/NCSS/Circular_Data_Analysis.pdf" class="uri">https://ncss-wpengine.netdna-ssl.com/wp-content/themes/ncss/pdf/Procedures/NCSS/Circular_Data_Analysis.pdf</a></p>
<p>[5] <a href="https://en.wikipedia.org/wiki/Von_Mises_distribution" class="uri">https://en.wikipedia.org/wiki/Von_Mises_distribution</a></p>
<p>[6] <a href="https://packaging.python.org/tutorials/distributing-packages/#your-package" class="uri">https://packaging.python.org/tutorials/distributing-packages/#your-package</a></p>
<p>[7] <a href="https://r-forge.r-project.org/scm/viewvc.php/pkg/R/vonmises.R?view=markup&amp;root=circular" class="uri">https://r-forge.r-project.org/scm/viewvc.php/pkg/R/vonmises.R?view=markup&amp;root=circular</a></p>
<p>[8] <a href="https://cran.r-project.org/web/packages/circular/circular.pdf" class="uri">https://cran.r-project.org/web/packages/circular/circular.pdf</a></p>