Circular Dot-plot.

The simplest representation of circular data is a circular raw data plot, in which each observation is plotted as a point on the unit circle.

roulette=c(43,45,52,61,75,88,88,279,357)
plot.circular(0,main="Circular Raw DotPlot of Roulette Data")
points.circular(rad(roulette))

Higher Moments.

Since \(e^{iα} = (\cos α + i \sinα)\) we note that,

\[ \bar{C}_n+i\bar{S}_n=\frac{1}{n}\sum_{j=1}^ne^{i\alpha_j} \]

Similarly we can define the higher order sample moments by taking higher powers of \(\left\{e^{i\alpha_j}\right\}\) and averaging these. We have,

\[ \frac{1}{n}\sum_{j=1}^n\left(e^{i\alpha_j}\right)^p=\frac{1}{n}\sum_{j=1}^n\cos p\alpha_j + i \frac{1}{n}\sum_{j=1}^n\sin p\alpha_j = \bar{C}_n(p) + i\bar{S}_n(p) \]

for \(p\in I^+\), where \(\left(\bar{C}_n(p),\bar{S}_n(p)\right)\) are called the \(p\)th order trigonometric moments based on the sample.

Circular PDF

A circular distribution is a probability distribution whose total probability is concentrated on the circumference of a unit circle. Since each point on the circumference represents a direction, such a distribution is a way of assigning probabilities to different directions or defining a directional distribution. \(\theta \in \left(0,2\pi\right]\) or \(\left(-\pi,\pi\right]\) ,

1. \(f\left(\theta\right)\geq0\)

2. \(\int_{0}^{2\pi}f\left(\theta\right)d\theta=1\)

3. \(f\left(\theta\right)=f\left(\theta+2k\pi\right)\) for any integer \(k\) (i.e. \(f\) is periodic).

Complex Fourier Series Representation of \(\boldsymbol{f\left(x\right)}\) with period \(\boldsymbol{2\pi}\).

Let \(f\left(x\right)\) be a periodic function of period \(2\pi\), then the complex fourier series representation \(f\left(x\right),x\in\left[-\pi,\pi\right]\) is,

\[ f\left(x\right)=\sum_{k=-\infty}^{+\infty}c_ke^{ikx} \]

where,

Characteristic Function.

As on the real line, such a distribution can also be described via its “characteristic function”. However, since \(\theta\) is a periodic r.v. having the same distribution as \(\left(\theta+2\pi\right)\), the characteristic function of such a r.v. has the property :

\[ \varphi_\theta\left(t\right)\overset{def}{=}\mathbb{E}\left(e^{it\theta}\right)=\mathbb{E}\left(e^{it\left(\theta+2\pi\right)}\right)=e^{it2\pi}\varphi_\theta\left(t\right) \]

Hence, either \(\varphi_{\theta}\left(t\right)=0\) or \(e^{it2\pi}=1\) i.e., \(t\) must be an integer. Thus, for circular r.v.s, the characteristic function needs to be defined only at integer values.

Trigonometric Moments of \(\boldsymbol{\theta}\)

The value of the characteristic function at an integer \(p\) is also called the \(p\)th trigonometric moment of \(\theta\).

\[ \varphi_\theta\left(t\right)=\mathbb{E}\left(e^{ip\theta}\right)=\mathbb{E}\left(\cos p\theta + i\sin p\theta\right) = \rho_pe^{i\mu_p}\,\,\,\, p\in\mathbb{Z} \]

Say,

\[ \alpha_p=\mathbb{E}\left(\cos p\theta\right),\beta_p=\mathbb{E}\left(\sin p\theta\right)\,\,\,p\in\mathbb{Z} \]

and these are related by the equation,

\[ \rho_p= \sqrt{\alpha^2_p+\beta^2_p}\,\,\,\mathrm{and}\,\,\mu_p=\arctan^*\left(\frac{\beta_p}{\alpha_p}\right) \]

Fourier Representation of \(\boldsymbol{f\left(\theta\right)}\) with period \(\boldsymbol{\theta}\).

Here the pdf \(f\left(\theta\right)\) is a periodic function of \(\theta\), so from the complex fourier series, we can write,

\[ c_p=\frac{1}{2\pi}\int_{-\pi}^{\pi}f\left(\theta\right)e^{-ip\theta}d\theta \]

\[ \iff2\pi c_p=\mathbb{E}\left(e^{-ip\theta}\right)=\varphi_{-p}\left(\theta\right) \]

That is,

\[ 2\pi f\left(\theta \right)=\sum_{k=-\infty}^{+\infty}\varphi_{-p}\left(\theta\right)e^{ip\theta}=1+\sum_{p=1}^{\infty}\left\{\varphi_{p}\left(\theta\right)e^{-ip\theta}+\varphi_{-p}\left(\theta\right)e^{ip\theta}\right\} \]

which gives,

\[ 2\pi f\left(\theta\right)=\left[1+2\sum_{p=1}^{\infty}\left(\alpha_p\cos p\theta + \beta_p\sin p\theta\right)\right] \]

Indeed, one can recover the pdf from either of these sequences. These trigonometric moments of \(\theta\) are the same as the Fourier coefficients in the Fourier series expansion of the pdf \(f\left(\theta\right)\). This is analogous to the “inversion formula” for characteristic functions of real-valued r.v.s.

Circular Uniform Distribution.

If the total probability is spread out uniformly on the circumference of a circle, we get the Circular Uniform (CU) distribution with the constant density,

\[ f\left(\theta\right)=\frac{1}{2\pi}\boldsymbol{1}\left[0\leq\theta\leq2\pi\right] \]

Circular Uniform Distribution.

If the total probability is spread out uniformly on the circumference of a circle, we get the Circular Uniform (CU) distribution with the constant density,

\[ f\left(\theta\right)=\frac{1}{2\pi}\boldsymbol{1}\left[0\leq\theta\leq2\pi\right] \]

Properties.

1. \(\theta\sim CU\left(0,2\pi\right)\implies\varphi_{p}\left(\theta\right)= \begin{cases}1 & \mathrm{for}\:p=0\\0 & \mathrm{for}\:p\in\mathbb{Z}\setminus\left\{0\right\} \end{cases}\)

2. \(\mu=\arctan^*\left(\frac{0}{0}\right)\) which is undefined. For Circular uniform distribution, there is no preferred direction.

3. \(\rho=\sqrt{\alpha^{2}+\beta^{2}}=0\)

4. \(\theta_{1},\theta_{2},\ldots,\theta_{n}\overset{iid}{\sim}CU\left(0,2\pi\right)\implies\sum_{i=1}^{n}\theta_i\sim CU\left(0,2\pi\right)\)

Circular Representation of \(\boldsymbol{CU\left(0,2\pi\right)}\)

cdensity(dcircularuniform,main="",nlabels=0)
segments(x0=0,y0=0,x1=1/sqrt(pi),y1=0) 
box(lwd=3)

Some Useful Results.

For \(\kappa>0,p\in\mathbb{Z}\),

1. \(I=\int_{0}^{2\pi}\left(\sin p\theta\right)e^{\kappa\cos\theta}d\theta=-I\implies I=0\)
2. \(I_{p}\left(\kappa\right)=\frac{1}{2\pi}\int_{0}^{2\pi}\left(\cos p\theta\right)e^{\kappa\cos\theta}d\theta.\, p^{th}\) order modified Bessel’s function of first kind.
3. \(I_{0}\left(\kappa\right)=\sum_{r=0}^{\infty}\left(\frac{\kappa}{2}\right)^{2r}\left(\frac{1}{r!}\right)^{2}\)
4. \(I=\frac{1}{2\pi}\int_{0}^{2\pi}e^{a\cos\theta+b\sin\theta}d\theta=\sum_{r=0}^{\infty}\left(\frac{1}{2}\sqrt{a^{2}+b^{2}}\right)^{2r}\left(\frac{1}{r!}\right)^{2}\)
5. \(\frac{\partial}{\partial\kappa}I_{0}\left(\kappa\right)=I_{1}\left(\kappa\right)\)

Circular Normal \(\boldsymbol{\left(CN\right)}\) / Von-Mises \(\boldsymbol{\left(\nu M\right)}\) Distribution.

A circular r.v. \(\theta\) is said to have a von Mises or a Circular Normal (CN) distribution if it has the density function :

\[ f\left(\theta;\mu,\kappa\right)=\frac{1}{2\pi I_0\left(\kappa\right)}e^{\kappa\cos\left(\theta-\mu\right)}\boldsymbol{1}\left[0\leq\theta<2\pi\right] \]

where \(0\leq\mu<2\pi\) and \(\kappa>0\) are parameters. Here \(I_{0}\left(\kappa\right)\) in the normalizing constant is the modified Bessel function of the first kind and order zero , and is given by,

\[ I_0\left(\kappa\right)=\frac{1}{2\pi}\int_0^{2\pi}\exp\left(\kappa\cos\theta\right)d\theta=\sum_{r=0}^{\infty}\left(\frac{\kappa}{2}\right)^{2r}\left(\frac{1}{r!}\right)^2 \]

Properties of \(\boldsymbol{CN\left(\mu,\kappa\right)}.\)

Properties of \(\boldsymbol{CN\left(\mu,\kappa\right)}.\)

Circular Representation of \(\boldsymbol{CN\left(0,\kappa\right)}\)

#---Circular Representation of CN densities----
#--mu=0, kappa= 0.5,1,2,4---
cdensity(function(x){dvonmises(x,mu=0,kappa=4)},nlabels=0,
main="Circular Representation of CN Densities",col="navy",xlim=c(-0.8,1.5),ylim = c(-0.8,0.8)) 
cdensity(function(x){dvonmises(x,mu=0,kappa=2)},nlabels=0,col="brown",add=T)
cdensity(function(x){dvonmises(x,mu=0,kappa=1)},nlabels=0,col="darkgreen",add=T)
cdensity(function(x){dvonmises(x,mu=0,kappa=0.5)},nlabels=0,col="orange",add=T)
legend("bottomright",legend=c("0.5","1","2","4"),col=c("orange","darkgreen","brown","navy")
,title=expression(kappa),lwd=2)
segments(x0=0,y0=0,x1=1/sqrt(pi),y1=0) 
box(lwd=3)

Properties of \(\boldsymbol{CN\left(\mu,\kappa\right)}.\)

Linear Representation of \(\boldsymbol{CN\left(0,\kappa\right)}\)

#---Linear Representation of CN densities----
#--mu=0, kappa= 0.5,1,2,4---
cols=c("0.5"="orange","1"="darkgreen","2"="brown","4"="navy") 
ggplot()+xlim(rad(c(-180,180)))+ylim(c(0,0.8))+geom_vline(xintercept=0,lwd=0.8)+
labs(y="Density",x=expression(theta),title = "Linear Representation of CN Densities",
caption=expression(paste(mu,"= 0")))+  
geom_function(fun=function(x){dvonmises(x,mu=0,kappa=0.5)},aes(col="0.5"),lwd=1)+
geom_function(fun=function(x){dvonmises(x,mu=0,kappa=1)},aes(col="1"),lwd=1)+
  geom_function(fun=function(x){dvonmises(x,mu=0,kappa=2)},aes(col="2"),lwd=1)+
geom_function(fun=function(x){dvonmises(x,mu=0,kappa=4)},aes(col="4"),lwd=1)+ 
 scale_color_manual(expression(kappa),values = cols) 

Variations of \(\boldsymbol{CN\left(\mu,\kappa\right)}.\)

CN Distribution on part of the Circumference.

One can define a CN distribution restricted to any arc of arbitrary length say \(\frac{2\pi}{\ell}\) for an integer \(\ell\), rather than the full circle, with the resulting density,

\[ g\left(\theta\right)=\frac{\ell}{2\pi I_0\left(\kappa\right)}\exp\left[\kappa\cos\ell\left(\theta-\mu\right)\right]\boldsymbol{1}\left[0\leq\theta<\frac{2\pi}{\ell}\right] \]

where \(0\leq\mu<\frac{2\pi}{\ell}\) .

• When \(\ell=2\), and the circular r.v. is restricted to a semi-circular arc, i.e., \(\theta\in\left[0,\pi\right)\), it is called an “axial distribution”. Occasionally, the measurements result in “axial data” as for instance when one observes the axes of pebbles rather than directions, and such distributions are relevant in modeling axial data.

Variations of \(\boldsymbol{CN\left(\mu,\kappa\right)}.\)

\(\boldsymbol{\ell}\)-modal CN Distribution.

A multi-modal CN density, say with \(\ell\) equidistant modes, is obtained through the density,

\[ g\left(\theta\right)=\frac{1}{2\pi I_0\left(\kappa\right)}\exp\left[\kappa\cos\ell\left(\theta-\mu\right)\right]\boldsymbol{1}\left[0\leq\theta<{2\pi}\right];\,\,\,0\leq\mu<2\pi \]

Variations of \(\boldsymbol{CN\left(\mu,\kappa\right)}.\)

Mixtures.

For better illustration we generate \(1000\) random samples from Mixed \(\nu M\) distribution with \(\left(\mu_{1},\mu_{2},\kappa_{1},\kappa_{2}\right)'=\left(\frac{3\pi}{2},\frac{\pi}{2},2,4\right)\) with proportion \(\frac{1}{2}\) and draw a circular histogram with the density curve.

sam=rmixedvonmises(1000,mu1=1.5*pi,mu2=pi/2,kappa1=2,kappa2=4,prop=0.5) 
chist(sam,nlabels=0,main="Circular Histogram for the Mixed Sample",border="orange") 
cdensity(function(x){dmixedvonmises(x,mu1=1.5*pi,mu2=pi/2,kappa1=2,kappa2=4,prop=0.5)},add=T)
segments(x0=0,y0=0,x1=1/sqrt(pi),y1=0) 
box(lwd=3)

Maximum Likelihood Estimation .

Multivariate \(\boldsymbol{\nu M}\) Distribution.

A unit random vector \(\boldsymbol{I}\) is said to have a \(p\)-variate von Mises-Fisher distribution if its p.e. is,

\[ C_p\left(\kappa\right)e^{\kappa\boldsymbol{\mu}'\boldsymbol{I}}dS_p\,\,\,I\in S_p \]

where \(\kappa\geq0\) and \(\boldsymbol{\mu}'\boldsymbol{\mu}=1\) . Here \(\kappa\) is called the concentration parameter and \(C_{p}\left(\kappa\right)\) is the normalizing constant. If a random vector \(\boldsymbol{I}\) has its p.e. of the form, we’ll say it’s distributed as \(\nu M_p\left(\boldsymbol{\mu},\kappa\right)\) . For the cases of \(p=2\) and \(p=3\), the distributions are called von Mises and Fisher distributions after the names of the originators ¹, and these particular cases explain the common nomenclature for \(p>3\).

Normalizing Constant.

It can be shown that,

\[ C_p\left(\kappa\right)=\frac{\kappa^{\left(\frac{p-1}{2}\right)}}{\left\{\left(2\pi\right)^{\frac{p}{2}}\boldsymbol{I}_{\left(\frac{p-1}{2}\right)}\left(\kappa\right)\right\}} \]

where \(I_r\left(\kappa\right)\) denotes the modified Bessel Function of the first kind and order \(r\), which can be defined by the formula,

\[ I_r\left(\kappa\right)=\frac{\left(\frac{\kappa}{2}\right)^r}{\Gamma\left(r+\frac{1}{2}\right)\Gamma\left(\frac{1}{2}\right)}\int_0^\pi e^{\pm\kappa\cos\theta}\sin^{2r}\theta d\theta \]

Random Sample Generation from \(\nu M_3\left(\boldsymbol{\mu},\kappa\right)\).

Simulation and density evaluation for \(p=3\). Taking \(\boldsymbol{\mu}=(0,0,1)',\kappa=0.5\)

mu <- c(0, 0, 1);kappa <- 0.5;n = 1e3;col=viridis(n)
x <- r_vMF(n = n, mu = mu, kappa = kappa)
rglobj=plot3d(x,box=F, col = col[rank(d_vMF(x = x, mu = mu, kappa = kappa))], size = 5, xlab = "x",ylab="y",zlab="z")
segments3d(x=c(0,0),y=c(0,0),z=c(0,1),lwd=2,col="red")
rglwidget()

Random Sample Generation from \(\nu M_3\left(\boldsymbol{\mu},\kappa\right)\).

Taking \(\boldsymbol{\mu}=(0,0,1)',\kappa=2\)

mu <- c(0, 0, 1);kappa <- 2;n = 1e3;col=viridis(n)
x <- r_vMF(n = n, mu = mu, kappa = kappa) 
rglobj=plot3d(x,box=F ,col = col[rank(d_vMF(x = x, mu = mu, kappa = kappa))], size = 5, xlab = "x",ylab="y",zlab="z") 
segments3d(x=c(0,0),y=c(0,0),z=c(0,1),lwd=2,col="red")
rglwidget()

Random Sample Generation from \(\nu M_3\left(\boldsymbol{\mu},\kappa\right)\).

Taking \(\boldsymbol{\mu}=(0,0,1)',\kappa=5\)

mu <- c(0, 0, 1);kappa <- 5;n = 1e3;col=viridis(n) 
x <- r_vMF(n = n, mu = mu, kappa = kappa)
rglobj=plot3d(x,box=F ,col = col[rank(d_vMF(x = x, mu = mu, kappa = kappa))], size = 5, xlab = "x",ylab="y",zlab="z",zlim=c(-1,1))
segments3d(x=c(0,0),y=c(0,0),z=c(0,1),lwd=2,col="red")
rglwidget()

Rayleigh Test for Uniformity.

Let \(I_{1},I_{2},\dots,I_{n}\) be a random sample of \(\nu M\left(\boldsymbol{\mu},\kappa\right)\). Consider,

\[ H_0:\kappa=0\,\,\mathrm{against}\,\,H_1:\kappa\neq0 \]

where \(\boldsymbol{\mu}\) is unknown , let \(\lambda\) be the LR statistics. Then it can be shown that,

\[ \log\lambda=-n\left\{\log\left(C_p\left(\hat{\kappa}\right)\right)+\hat{\kappa}A\left(\hat{\kappa}\right)\right\} \]

Note that \(\lambda\) is a function of \(\bar{R}\) alone. It can be shown that \(\lambda\) is a monotonically decreasing function of \(\hat{\kappa}\). Finally , on differentiating with respect to \(\hat{\kappa}\) and using results we get , it is found that \(\hat{\kappa}\) is a monotonically increasing function of \(\bar{R}\). Therefore the critical region of the Rayleigh test reduces to \(\bar{R}>k\).

For large \(n\) it can be shown that,

\[ pn\bar{R}^2\overset{H_0}{\sim}\chi^2_p \]

About the Dataset.

The circular package in R contains a dataset named pigeons. This has \(108\) rows and \(2\) columns. The observations are the vanishing bearings of homing pigeons displaced and released at two unfamiliar locations. The data are pooled with respect to the home direction (home direction set in \(360\) grades). This contains the following column,

• treatment : A factor with levels:

  • c : Control pigeon (unmanipulated).

  • v1 : Pigeons subjected to bilateral section of the ophthalmic branch of the trigeminal nerve.

  • on : Pigeons subjected to bilateral section of the olfactory nerve.

• bearing : Vanishing bearing of each bird in degrees.

The bearings are in degrees, we need to change them into radians.

pigeons$bearing=rad(pigeons$bearing)

Exploratory Data Analysis.

With respect to the factors, we first make a multiple stacked histogram of bearing, to get an idea about the data.

with(pigeons,cmhist((bearing),treatment,nlabels=0,main="Stacked Histogram of Bearing"))
segments(x0=0,y0=0,x1=1/sqrt(pi),y1=0) 
box(lwd=3,col="navy")

Exploratory Data Analysis Contd.

Let’s obtain the summary statistics.

with(pigeons,circ.mean((bearing))) 

[1] 0.1742687

with(pigeons,circ.disp((bearing))) 

    n        r      rbar       var
1 108 61.81342 0.5723464 0.4276536

Here the variability is much high. A natural question may be,

• Does the pigeons have any preferred direction ?.

Before going into that, we need to look for appropriate distribution of bearing. Will \(CN\left(\mu,\kappa\right)\) give a better fit for this data for some suitable \(\mu,\kappa\) ?

Testing for Uniformity.

The proposed model gives a nice fit. It’s plausible to work with CN distribution in this case. Now, we check whether the pigeons have any preferred direction by rayleigh test.

(rayleigh.test(circular(pigeons$bearing)))

       Rayleigh Test of Uniformity 
       General Unimodal Alternative 

Test Statistic:  0.5723 
P-value:  0 

This test strictly rejects \(H_{0}\) at level \(\alpha=0.05\). The pigeons have a preferred direction and an estimate of the direction is \(0.1743\) radians towards north-east.

References.

• Directional Statistics - Kanti V. Mardia, Peter E. Jupp - John Wiley & Sons.

• Multivariate Analysis - KV Mardia, JT Kent, JM Bibby - Academic Press

• Topics in Circular Statistics - Sreenivasa Rao Jammalamadaka, A. Sengupta - World Scientific.

• Circular Statistics in R - Arthur Pewsey, Markus Neuhauser, Graeme D. Ruxton - Oxford University Press.

• FISHER, N. I. (1993), Statistical Analysis of Circular Data, Cambridge University Press, Cambridge.

• CSORGO, M. AND HORVATH, L. (1996), A note on the change-point problem for angular data, Statist. Probab. Lett, 27, 61—65.

• Wikipedia