,B tI 1997 uKS.IAG IsRARY RESEARCH SUPPORT UNIVERSITY OF MICHIGAN BUSINESS SCHOOL DECEMBER 1996 A FULL BAYESIAN ANALYSIS OF CIRCULAR DATA USING VON MisEs DISTRIBUTION WORKING PAPER #9612-37 PAUL DAMIEN UNIVERSITY OF MICHIGAN BUSINESS SCHOOL AND STEPHEN WALKER IMPERIAL COLLEGE, LONDON, UK

I A Full Bayesian Analysis of Circular Data using von Mises Distribution Paul Damien1 and Stephen Walker2 Department of Statistics and Management Science, School of Business, University of Michigan, Ann Arbor, 48109-1234, USA. 2 Department of Mathematics, Imperial College, 180 Queen's Gate, London SW7 2BZ SUMMARY This paper presents a full Bayesian analysis of circular data paying special attention to the von Mies distribution. We obtain samples from the posterior distribution, given an independent sample from the von Mises distribution, using the Gibbs sampler after the introduction of some strategic latent variables which ensures all the full conditional distributions are of known type. Keywords: Circular data, latent variable, Gibbs sampler, von Mises distribution. I

1 Introduction The am of this paper is to provide the basis for a full Bayesian analysis for circular data. In particular we will be concentrating on the von Mines distribution (von Mises, 1918) though our method of analysis should also be applicable to spherical and cylindrical data. Circular data arise naturally in a number of areas and for a recent survey from a frequentist perspective the reader is referred to the book of Fisher et al. (1987). Frequentist approaches have also been Arnold (1941), Fisher (1953), Gumbel (1954), Mardia (1972, 1975) and Bingham (1964). Some theoretical results in this context also appears in Lvy (1939). The Bayesian literature is far less extensive and to some extent has not been very successful. This has been due to the difficulties in working with the parametric distributions commonly associated with circular data, in particular, the von Mises distribution. The von Mises distribution is the most important distribution in the statistics of circular data and is the 'natural' analogue on the circle of the normal distribution on the real line, The earliest attempt at Bayesian inference for the von Mises distribution was given by Mardia and El Atoum (1976). However these authors assume the concentration parameter to be known and only provide point eatimatos for the directional parameter. Guttorp and Lockhart (1988) consider the case when both parameters are unknown but are forced to use the posterior maximum likelihood (mode) estimate for the concentration parameter (Lenth, 1981). The Bayesian analysis of Bagchi (1987) and Bagchi and Guttman (1988) also restrict their attention to the von Mines distribution. More recently Bagchi and Kadane (1991) developed Laplace approximations to posterior distributions involving the von Mises distribution. They only provide Bayes estimates (posterior means) for the cosine of the directional parameter after taking the concentration parameter to be known. In this paper we implement a full Bayesian analysis involving the von Mises distribution in which both parameters are assumed unknown and assigned the conjugate prior distribution introduced by Guttorp and Lockhart (1988). We derive the posterior distribution and via the Introduction of strategic latent variables (Damien and Walker, 1996) are able to use a Gibbs sampler (Smith and Roberts, 1993) to simulatefrom this postcrior, providing the basis for a full Bayesian analysis. 2

I 2 The von Mises distribution The von Mises distribution is a symmetric unimodal distribution which is the most commontmodel for unimodal samples of circular data. The probability density function is given by f(lp, c) = [2tr1o()]-exp[cos(0 - )] 0 < 0 < 2, c 2 0, where Io() = (2vr)-1' exp[n co. J]d is the modified Bessel function of order zero. The mean direction is # and a is the concentration parameter, Prior distribution We use the conjugate prior suggested by Guttorp and Lockhart (1988). This is given, up to a constant of proportionality, by f(p, c) oc l "(,)ecxp[,Ro&os(e - po)]. We will only be considering the case when c is a non-negative integer which is in the spirit of the interpretation of the prior parameters. Esaentially the prior parameters can be thought of as representing c observations in the direction po and RA can be thought of as the component on the z axis (i.e., in the known direction) of the resultant of c observations. The uniform distribution on the circle is a limiting case of this prior distribution. Posterior distribution Let 0 = (0l... ^,) be a sample of size rn. The posterior distribution of (Is, t) is given, up to a constant of proportionality, by f(J, K0) oc ) e o xp(, () cos(S -,)1, where m = c + n and,, and R, are obtained from I, cos,, = R coS /o + Co O; and, sin /, = R sin Po + E sin,. In the next section we develop a Gibbs sampler for generating random variatee from the posterior. 3

I 3 Inference via the Gibbs sampler Our approach will depend on the introduction of latent variables to define a joint distribution with (p/, ). This joint distribution will be constructed to ensure that all full conditional distributions (required for the Gibbs sampler) are of known type and can be sampled directly. Damien and Walker (1996) exemplify this concept in a variety of applications after developing suitable theory, and Walker and Damien (1996) extend the ideaa to the analysis of data in the context of neutral to the right processes. Consider firt the latent variablec w, with w defined on (0,oa), and v, also defined on (0,oo), and define their joint distribution with (p, P) by f(F,, xw, u) ac &l1(v < R^^[l+ac(P-n)) (Wm-c-wo(n)} Clearly the marginal distribution for (,/, n) is as required. Our next step involves writing Io(x) in the form 00 1o(n)-, where Ak = (k!)~O.52k. Therefore the joint distribution becomes f(i,,w,ii) c e^I(v <c eRfttl+c0(4->M)]){wm-1 ft c-^ }. kwo Next we introduce the latent variable u = (ul, u2,...) and r and define the joint distribution with (ii, w, v) by f(P,.,w,tWIu, ) ac eAI l~(y < eRnnh +o#-?C)J, r < wml)(e i I(u < e-cG k)}. trAgain it is clear that the marginal for distribution for (#j, x) is u required. To implement the Gibbs sampler we need the full conditional distributions (denoted by a *). Also the symbol U denotes the uniform distribution. The full conditional for z is given by r(r) = U(O, w-1) 4

I and for v is given by V= U(o, R^"o(L-( nl)J The full conditional for p is given by f*() = U(A), where A is given by the set pn +cosl-'[(( ^)~ log u - 1]. The full conditional for w is given by f'(w) c c-W('l1( -') < w < mink{-(Akt)-c loguk}). and the full conditional for K is given by f'(c) oc e-"" I(maxf{O, } < n < minkm{-(wA)-llogtJ U1f(2'k)), where v, - log v/[RE +, cofs( - /)]. Lastly, the full conditional for ech of the uk is given by fi (U)= U (0e=g- ). It is clear that we cannot sample all the us since there are an infinite number of them. However we do not need to sample all of them in order to implement the Gibbs sampler (recall it is the (p, r) samples we are interested in). In fact it it easy to see from the full conditionals of w and c that we only need to obtain mink{-(Akrck)-t 1ogui} and min^{[-(wAh)-l loguJ]t/(k^). We will only concern ourselves with the former of these as the latter will follow by essentially the same method. Since the full conditional for ui is the uniform distribution on (0, e""L) we can take ui - r=e- T X^^ where 1, are iid U(O, 1) and t is the 'current' w. Therefore the upper bound for the (exponential) full conditional distribution of w is given by mini{w - aA(c) log T}, where ack(c) = (Ah,92)-1, and note that ae -+ oo. Essentially then we only need to consider {/a(cx) = -ak(K) log rT}. In figure 1 we illustrate ((Kr) (on log scale) for s = 1, 5, 10,20,40 and 80 and for each of these it is clear that obtaining the minimum value is very tractable. From figure 1 we observe that the minimum value seems to occur when k is approximately K/2 (certainly for the larger values of x.) This can be 5

I explained theoretically. Our interest in the following discussion is in understanding _ = k: Mn = — a(&c)logmr where Mx = mink{-ac(c) log n}). Recall that the rks are iid U(0,1) and so clearly the 'expected' value of k will be that k minimising (k!)2(n/2)'23. Taking logarithms we see that this reduces to finding that k which minimisea log k! - klog(i/2). If we allow k to take non-integral values and replace k! by r(k + 1), where r represents the gamma function, then we want that k such that b(k + 1) = log(c/2) (here i is the digamnma function). According to Abramowitz and Segun (1968) the asymptotic formula for,k is given by #(k) = logk - 1/(2k) - 1/(12k)2 ^ - and so a good approximation for the required k (in particular for large c) is provided by t/2. To investigate this point further and to see what the variance of k looks like we simulated 1000 Mso random variables and the resulting histogram representation of the sample is given in figure 2. The mean value of the sample is 39.8 and the standard deviation is 4.6. 4 Numerical examples Example 1. Here we consider the von Mises distribution and the roulette data with sample size 9 (see Mardia 1972, p.2). We simulated from the joint posterior distribution for (p, t) using the algorithm described in Section 3. We take c - 0, 0o s 0 and Re = 0 and the resulting marginal posterior distributions for cosp and t are presented in figure 3. Example 2. We also consider the data set presented in Guttorp and Lockhart (1988, Table 1.) Seven signals from an unknown location (the target) are picked up at observed compass readings by receiver stations at known locations. If A is the true direction from station i to the target then f/(4 p, Pc) - [2TIo(.)]l expf[ Coa( - A)] Therefore only minor modifications arc required to the algorithm outlined in the previous section. We take a flat prior for the j and take f(c) o J7 (Pc) exp(Ro-c) as the prior for i with (Ro, c) = (5,5) (see Guttorp and Lockhart for this prior and alternatives). The posterior plot of the location of the target is given in figure 4. 6

I 5 Discussion In this paper we have provided a full Bayesian analysis of circular data modelled via the von Mises distribution. It is quite straightforward to extend the idea of latent variable substitution resulting in full conditional distributions of known type in the Gibbs sampler to observations on the surface of a p-dimensional hypersphere S, of unit radius and centered on the origin. This would make use of the von Mises-Fisher distribution (see, for example, Mardia and El-Atoum, 1976) with density function given by f(1,rl P) = c,() exp[cp'0, where as before p is the mean direction, s the concentration parameter and cp(r) a= i(-t)//[(21r)/2Ip,_l)/] is the constant factor where I,(n) is the modified Bessel function of the first kind given (which is important to us) as the series ) = (r + k + l)r(k + 1)' where r is the gamma function. Therefore we can proce'd in the same way as for the circulax data using the conjugate or noninformative priors proposed by Kardia and El-Atoum (1976). 6 References Abramowitz,M. and Segun,I.A. (1968). Handbook of Mathematical TFnctions, Dover Publications Inc., New York, Arnold.K.J. (1941). On spherical probability distributions, Ph.D. Thesis, Massachusetts Institute of Technology. Bagchi,P. (1987). Bayesian analysis of directional data. Ph.D. Thesis, University of Toronto. Bagchi,P. and Guttman,I. (1988). Theoretical considerations of the von Mises-Fisher distribution. Journal of Applied Statistics 15, 149-169. 7

Bagchi,P. and Kadane,J.B. (1991). Laplace approximations to posterior moments and marginal distributions on circles, spheres, and cylinders. The Canadian Journal of Statistics 19, 67-77. Bingham,C. (1964). Distributions on the sphere and on the projective plane. Ph.D. Thesis, Yale University. Damien,P. and Walker,S.G. (1996). Sampling probability densities via uniform random variables and a Gibbs sampler. Submitted for publication. Fisher,N.I.,Lewis,T. and Embleton,B.J.J. Statistical analysis of spherical data. Cambridge University Press. Fisher,R.A. (1953). Dispersion on a sphere. Proceedings of the Royal Society A217, 295-305. Gunbel,E.J. (1954). Applications of the circular normal distribution. Journal of the Amer;can Statistical Association 48, 233-246. Guttorp,P. and Lockhart,R.A. (1988). Finding the location of a signal: A Bayesian analysis. Journal of the American Statistical Association 83, 322 -329. Lenth,R.V. (1981). On finding the source of a signal. Technometrics 23, 149-154. LevyP, (1939). L'addition des variables aldatoires defines sur une circonferences. Bull. Soc. Math. funce 67, 1-41. MardiaK.V. (1972). Statistics of directional data. London: Academic Press. Mardia,K.V. (1975). Statistics of directional data. Journal of the Royal Statistical Society Series B 37, 97-133. Mardia,K.V. and El-Atoum,S.A.M. (1976). Bayesian inference for the von Mises-Fisher distribution. Biometrika 63, 203-205. 8

Von Mises,R. (1918). Uber die "GAnzzahligk-eit" der Atomgewicht und Verwandte Fragen. Physikal. g. 19, 490-500. Smith,A.F.M. and Roberts,G.O. (1993). B&yesian computations via the Gibbs sampler and related Markov chain Monte Carlo methods. Journal of the Royal Statistical Society Series B 55, 3-23. WalkerS.G. and Damien,P. (1996). A full Bayeaian analysis involving a neutral to the right process. Submitted for publication. Q

1 I I I II kLd 3 I a S I I a4 -j r, j I tI 1 1 i # I'mbN I I ir II q It 7I 1 - I I _ # l at 1 I ti at I " It es * q # w1 n Figure 1: Plots of log beta (k,x) versus k for x=1,5,10,20,40 and 80 10

8. 0 d q W................ 0 i- -— I! i I I — I — i ii - I 25 30 35 40 45 50 55 M_80 Figure 2: Histogram approximation of the density for Mso I 11

I T l -"*-*@ W & -- ""'""'^"^i*l *'n.......... —. —1 —1 - to C 7t ~ ~ ~ ~ ~ ~~bpp 0 Figure 3: Hiltogram approximation of the marginal posterior densities for coBp (top) andt P (bottom) 0~~~~~~~~~~~~o~u 0' COBp (top) and,c (bottom) 12

m -M * I a, a I 0* ci, C I 1 I *39 4%4 I II h of fie 0 g4 s0 * Pg I L Ia A 4 CI0 a I I.. 1 ' I 0 1 w ci 0 v # 3 a qUI I. * *4 *, * * C. I 1 C U,4 0R PI I -- - - T- r i- I 5.5 6.0 OS5 7.0 7.5 6.0 x (meat vakn 73) Figure 4: Positerior plot of location of the target asE 13