Taking The Average Of A Dihedral Angle
4
1
Entering edit mode
12.7 years ago
dimkal ▴ 730

I have a time series a dihedral angle measured from a Molecular dynamics simulation. This dihedral angle fluctuates around 0 degrees. How do i properly take an average of something that's either +1 and +359 degrees, and take the periodicity of the system in to account to get zero???

I guess i have to use polar coordinates, but entirely sure on how to do this.

• 6.0k views
ADD COMMENT
2
Entering edit mode
12.7 years ago

I do not know of a way to do this with a single formula, but you can solve it as an optimization problem.

Let us call the unknown average x and the measured angles x_1, x_2, x_3, etc. For a given value of x, you can define a sum of squared errors, which is easy enough to calculate in polar coordinates:

error = sum_i (diff(x, x_i)**2), where diff(x, x_i) = min(abs(x-x_i), 360-abs(x-x_i))

You now simply find the value of x that minimizes this error function using whatever numeric optimization algorithm you like.

ADD COMMENT
1
Entering edit mode

yes that seems like it would work, and this is what i'm sort of doing at the moment. However I'm looking to see if there's a much simpler approach.

ADD REPLY
2
Entering edit mode
12.7 years ago
Ahill ★ 2.0k

Editing this answer, I'll suggest a solution to the optimization problem posed by Lars. The "average" angle phi that minimizes the cartesian distance to the observed angle unit vectors on the unit circle would be:

phi = atan(sum_i(sin(theta_i))/sum_i(cos(theta_i)))    [if sum_i(cos(theta_i))>=0]
      atan(sum_i(sin(theta_i))/sum_i(cos(theta_i)))+pi [otherwise]

where theta_i are the observed angles (in radians).

In R:

mean.angle <- function(theta.i) {
  sum.sin <- sum(sin(theta.i))
  sum.cos <- sum(cos(theta.i))
  mean.angle <- atan(sum.sin/sum.cos)
  phi <- ifelse(sum.cos>=0, mean.angle, mean.angle+pi)
}
ADD COMMENT
2
Entering edit mode

I don't think that will work. If you tage the average of 1 and 359 degrees you will get 180, but the real answer would be 0.

ADD REPLY
1
Entering edit mode

Oops i made a mistake in my question, it should of been +359 and not -359. in which case your algorithm would average it as ~180, but it should be zero.

ADD REPLY
0
Entering edit mode

but the correct answer should be zero.

ADD REPLY
0
Entering edit mode

I see I've misunderstood the question, agree original answer will not work. See edit to answer above.

ADD REPLY
1
Entering edit mode
12.7 years ago

If you convert to the interval [0,2pi], you should get a mean of pi. Mean-centering that data by subtracting pi should give you data on the interval [-pi,pi], where the mean of a uniform distribution is 0. You could also try that with +/-180, if the conversion is too time consuming. The difference in sign corresponds to mechanical interpretations of opposite forces. I would also try to convert to flory convention, depending on the polymer you're working with. Hope I didn't misunderstand your question, I'm interested a more thorough explanation.

ADD COMMENT
0
Entering edit mode

Sorry, but I think you have misunderstood the problem. Suppose I start out with two angles of 1 and 359 degrees. After your transformations, the average becomes 0, again corresponding to 180 degrees whereas the right answer is 0 degrees.

ADD REPLY
0
Entering edit mode
12.7 years ago
Kris • 0

What if you try the following transformation of your angles (in degrees):

[0,360) ---> (-180, 180]

newangle = oldangle, if oldangle is between 0 and 180 degrees

newangle = oldangle - 360, if oldangle is between 180 and 360 degrees

(In other words, if the point of the unit circle lies in quadrants III and IV (and only then), you choose a different angle to represent this point, an angle that differs from your angle by 360 degrees and is negative)

With this transformation, angle +1 stays angle +1, and angle +359 becomes angle -1, so the two angles average to 0. The uniform distribution on [0,360) gets mapped to the uniform distribution on (-180,180] and the new average is 0. Does that work?

ADD COMMENT

Login before adding your answer.

Traffic: 2006 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6