🚀 go-pugleaf

   Let I (kappa) = Sum(r=0,oo) ((kappa/2.0)^(2*r+p))/((r!)*((r+p)!))
        p
be the modified Bessel function of order p, with power series as given above.

   The von Mises distribution on a circle {theta in [0, 2*pi]  or  theta in
[-pi, pi] }    is given by
vM(theta)  =  (1/C(kappa))*exp(kappa*(cos(theta-mu))) ,
                                           where C(kappa) = 2*pi*I (kappa) .
                                                                  0
The von Mises distribution around the circle has similarities to the
Gaussian distribution on the line.   For the von Mises, mu is the mean
direction, and kappa is a spread parameter (like the Gaussian 1/(sigma**2) ).

Let A(kappa)= E(cos(theta)) = I1(kappa)/I0(kappa)   (not lowering subscripts).

   Given the rapidly (faster than exp(.) ) converging power series for I0 and
I1 above,  A is easy to calculate.
   The FIRST (of the two) PURPOSE(s) of this posting is to find out whether
there is a good quick way to calculate A (without having to use I1 and I0)
or a good quick way to calculate A_-1, the inverse function of A.
(I am currently doing it by binary chop.  My problem is that, in calculating
I0 and I1 (vaguely like exp(.) ), my machine gets into trouble for large
inputs to A  (something which one of the estimators below tends to require)
and I get a core dump.)

   The mean, mu, is easy to estimate, and all reasonable estimators agree.
Let x = Sum(i=1,n) cos(theta_i)   and    y = Sum(i=1,n) sin(theta_i)    .
if (x .eq. 0.0) then mu_hat = pi/2.0 else mu_hat = atan(y/x) endif .

   The SECOND (of the two) PURPOSE(s) of this posting is to find some
estimators of kappa.  (We are working on one, and I (only) know of two others.)

This completes my requests.  The remainder of this posting gives some of the
estimators for kappa (one of which tends to require large inputs to A and
hence, by my current method, to I0 and I1, and highlights my machine overflow
problem).
   To get these estimators, we must firstly define r .
Let phi(i) = theta(i) - mu     for i = 1,...,n .
Then, r = sqrt(x^2 + y^2) = Sum(i=1,n) cos(phi(i))   {equivalent defn's}.

   The Maximum Likelihood Estimate (m.l.e) kap_ml of kappa is given by
kap_ml = A_-1(r/n),  or (equivalently)  A(kap_ml) = r/n .
{Reference: Upton+Fingleton,p229,(9.30) .}

   A second ("pendulum") estimator {we say "pendulum" because there is a
thermodynamic analogy between this estimator and a system of pendula in
equilibrium}, kappa_tilde (or kappa~ ) , satisfies
(Upton+Fingleton,p230,eqn(9.32)) :
  n*A(kappa_tilde) = r*A(r*kappa_tilde)  }     (r**2)/n >= 1
  kappa_tilde      = 0                   }     (r**2)/n <  1

(The r*kappa_tilde input to A (and to I1 and I0) above is what is currently
causing (my machine and) me problems.)

   For r/n > 0.9 (ish),  Upton+Fingleton gives an approximation to kappa~ .
(Upton+Fingleton,p230,(9.33))    (for r/n > 0.9 (ish) ),
ktaprx =   (1.0 - 1.0/n)/(2.0*(1.0-r/n))
         + (r/n - 1.0/(n**2))/(4.0*(r/n)*(1.0-1.0/n)) .

References: Upton, G.J.G. + B. Fingleton, "Spatial Data Analysis by Example
(Categorical and Directional data)", Volume 2.
The book by K.V. Mardia is also very helpful and relevant (and contains the
Ip power series), and a book by N.I. Fisher, T. Lewis and B.J.J. Embleton
(of which I am told there is to be a sequel) is also of some help.

   Re-capping, comments on how best to calculate A or (the inverse of A) A_-1 ,
with particular relevance to preventing machine overflow from the "pendulum"
estimator would be appreciated.  References to alternative estimators for
kappa would also be welcomed.

Thank you.         - David Dowe.

(Dr.) David Dowe, Dept of Computer Science, Monash University, Clayton,
Victoria 3168, Australia           dld@bruce.cs.monash.edu.au