From: Ravi Varadhan <rvaradha_at_jhsph.edu>

Date: Wed 29 Jun 2005 - 05:36:22 EST

Ravi Varadhan, Ph.D.

Assistant Professor, The Center on Aging and Health Division of Geriatric Medicine and Gerontology Johns Hopkins University

Ph: (410) 502-2619

Fax: (410) 614-9625

Email: rvaradhan@jhmi.edu

R-help@stat.math.ethz.ch mailing list

https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html Received on Wed Jun 29 05:50:15 2005

Date: Wed 29 Jun 2005 - 05:36:22 EST

Hi Mike,

You should think of "vector" sums rather than arithmetic sum. So the mean is really the direction of the resultant vector of all the unit vectors. When two unit vectors are exactly opposing each other (as in your example with angles 45 and 225 degrees), they cancel each other out and the resultant vector has length zero, and hence the direction is arbitrary and undefined. However, there are two issues which give rise to seemingly anomalous results: (1) due to finite numerical precision in computing the sines and cosines, exact cancellation does not happen, and (2) discontinuity in atan function due to branch-point (0, 0) and a branch cut along negative real axis.

For example,

> deg(circ.mean(c(rad(30),rad(210))))

[1] -56.30993247402022

*> sx <- sin(rad(210)) + sin(rad(30))
**> cx <- cos(rad(210)) + cos(rad(30))
**> sx
*

[1] -1.665334536937735e-16

*> cx
*

[1] 1.110223024625157e-16

*> atan(sx,cx) # this should really be 0
*

[1] -0.9827937232473291

*>
*

If you now look at the resultant vector it has components (sin(-0.9828),
cos(-0.9828)) = (-0.83,0.55), which is completely wrong.

However, the above example is only a minor perturbation away from the
following (where the resultant vector is not zero):

*> deg(circ.mean(c(rad(211),rad(30)))) # note 211 = -149 (mod 360)
*

[1] -59.49999999999996

which is perfectly alright.

Hope this helps,

Ravi.

Ravi Varadhan, Ph.D.

Assistant Professor, The Center on Aging and Health Division of Geriatric Medicine and Gerontology Johns Hopkins University

Ph: (410) 502-2619

Fax: (410) 614-9625

Email: rvaradhan@jhmi.edu

> -----Original Message-----

*> From: r-help-bounces@stat.math.ethz.ch [mailto:r-help-
**> bounces@stat.math.ethz.ch] On Behalf Of Michael Peckford
**> Sent: Tuesday, June 28, 2005 1:35 PM
**> To: r-help@stat.math.ethz.ch
**> Subject: [R] Circular Mean Question
**>
**>
**> Hi, a question about the circular mean function in the package
**> CircStats:
**>
**> Can anyone shed some light on why the circ mean function seems to make
**> sense for the first 2 set of bearings and then the mean of 225 and 45
**> degrees gives an unexpected 180 deg.
**>
**> > deg(circ.mean(c(rad(222),rad(45))))%%360
**> [1] 133.5
**> > deg(circ.mean(c(rad(224),rad(45))))%%360
**> [1] 134.5
**> > deg(circ.mean(c(rad(225),rad(45))))%%360
**> [1] 180
**> > deg(circ.mean(c(rad(226),rad(45))))%%360
**> [1] 315.5
**>
**> Can anyone explain this???
**>
**> This problem was first detected when I was trying to take the circ
**> weighted means of my data:
**>
**> With 2 groups of bearings:
**> x <- c(270,180)
**> y <- c(45,270)
**>
**> the circular mean of these bearings gives:
**> > deg(circ.mean(c(rad(x),rad(y))))%%360
**> [1] 257.2356
**>
**> When finding the weighted means I get this:
**> > meany <- circ.mean(rad(y))
**> > meanx <- circ.mean(rad(x))
**>
**> > deg(circ.weighted.mean(c(meanx,meany),c(2,2)))%%360
**> [1] 281.25
**>
**> The function for weighted mean I am using:
**>
**> circ.weighted.mean <- function (x,w)
**> {
**> sinr <- sum(w*sin(x))
**> cosr <- sum(w*cos(x))
**> circmean <- atan(sinr, cosr)
**> circmean
**> }
**>
**> I am assuming that the problem that mention above is the cause of the
**> different mean bearings.
**>
**> Am I missing something fundamental here?
**>
**> Thanks,
**> Mike
**>
**> ______________________________________________
**> R-help@stat.math.ethz.ch mailing list
**> https://stat.ethz.ch/mailman/listinfo/r-help
**> PLEASE do read the posting guide! http://www.R-project.org/posting-
**> guide.html
*

R-help@stat.math.ethz.ch mailing list

https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html Received on Wed Jun 29 05:50:15 2005

*
This archive was generated by hypermail 2.1.8
: Fri 03 Mar 2006 - 03:33:06 EST
*