4
$\begingroup$

I am having this issue with directional statistics. I have a probability distribution on a circle (in this case, an Orientation Distribution Function). This is for measuring the directionality of collagen fibers in soft tissue and therefore, the directionality doesn't matter as collagen facing one direction is the same as collagen facing that direction + 180 degrees.

I want to find the mean direction/orientation using directional statistics.


My progress

I've tried the moment formula available on the directional statistics Wikipedia page.

Directional Statistics Info - Courtesy of Wikipedia

However, the vector I have found is not the correct preferred direction vector. I was wondering if any one had any ideas on what I should be doing here.

[Top: The ODF plotted along an x-axis of degrees. Bottom: The ODF wrapped around a circle with my attempt at finding the preferred direction vector

Edit - Added more detail per the request of Whuber

As stated above, I am trying to use the moment equation provided by wikipedia. This is the first moment (the mean vector) equation in integral form. In the first step, I take my ODF (which I am using a Fourier series to represent, though this is irrelevant for the question (I think)) and wrap it around the circle with the Euler's identity.In the second step, I find the mean vector by taking the integral of the circular data from 0 to 2pi.

The math that I am doing right now

This is my current code (in MATLAB).

% Wrap the data around the circle

circ = (exp(1i*examplePoint.theta));

% Theta is angles starting at Pi/180 through 2pi

vals = abs(examplePoint.odf).*circ;

% examplePoint.odf is the ODF values. 
% The absolute value is just to make sure all ~0 values are 
% on the positive side.

% Find the angle of the preferred direction vector

m1 = sum(vals)/(examplePoint.theta(end)-examplePoint.theta(1));

% My integral

PrefD = angle(m1);

% Calculates the angle of the vector the integral found

Then I plot the ODF and a vector of angle PrefD on the polar plot (shown).

The problem is that I can rotate the data any direction and the angle is always within 10 degrees but never correct. It is always slightly off to the right or left, seemingly randomly.

I have tried different ways of calculating the integral and I always get more or less the same results.

$\endgroup$
9
  • 2
    $\begingroup$ Could you please show exactly what vector you did find, explain how you found it, and describe how you know it is incorrect? $\endgroup$
    – whuber
    Feb 7 '17 at 18:12
  • $\begingroup$ Sure! I added some more detail to the question. $\endgroup$
    – user148316
    Feb 7 '17 at 21:40
  • $\begingroup$ I don't see where you actually followed out the Wikipedia prescription: how does your calculation of m1 actually integrate anything? It seems to be a sum divided by a different of angles--but that won't even approximate the integral in question. $\endgroup$
    – whuber
    Feb 7 '17 at 22:41
  • $\begingroup$ I've tried it this way and using simpson's rule to approximate the integral and they both get me around the same results. I guess simpson's rule is the more accurate way to do it as it is a more accurate integral approximation but they both get about the same thing. $\endgroup$
    – user148316
    Feb 8 '17 at 1:21
  • $\begingroup$ I don't see how Simpson's rule applies to data. Could you show, perhaps with a tiny dataset (two or three angles will do) how you are attempting to calculate the mean direction? $\endgroup$
    – whuber
    Feb 8 '17 at 1:54
2
$\begingroup$

You don't have an ODF: you have empirical data. Computing the Fourier Transform will work, if carried out properly, but it's redundant and inefficient. It amounts to finding the mean direction by separately averaging the x- and y-components of the directions and then converting that average to a direction.

For orientation data, all you need to do is double all angles. This identifies every direction with its opposite direction. After computing the mean direction of the doubled-angle data, halve its angle (so that the result will be between $-\pi/2$ and $\pi/2$). This designates the mean orientation.

Figure

This figure shows ten random orientations and their mean orientation in red.

The following R code shows how the figure was created. The mean.orientation function performs the calculation. It allows for weighted data via its weights argument as shown in the commented-out line for computing alpha.bar.

mean.orientation <- function(alpha, weights) {
  theta <- 2*alpha
  if (missing(weights)) weights <- rep(1, length(alpha))
  atan2(mean(weights*sin(theta)), mean(weights*cos(theta))) / 2
}

n <- 10
alpha <- rnorm(n, pi/7, pi/15) + pi*sample.int(2, n, replace=TRUE)
weights <- rexp(n)
#alpha.bar <- mean.orientation(alpha, weights)
alpha.bar <- mean.orientation(alpha) # For equally-weighted data

u <- seq(0, 2*pi, length.out=360)
w <- max(weights)
circle <- w*cbind(cos(u), sin(u))
plot(circle, type="l", col="#e0e0e0", asp=1, bty="n", xlab="", ylab="")
arrows(-weights*cos(alpha), -weights*sin(alpha), 
       weights*cos(alpha), weights*sin(alpha), length=0)
arrows(-w*cos(alpha.bar), -w*sin(alpha.bar), 
       w*cos(alpha.bar), w*sin(alpha.bar), length=0, lwd=2, col="Red")
$\endgroup$
3
  • $\begingroup$ The difference is that you're assuming all angles have the same intensities. In my problem, the angles have different intensities (or, more precisely, relative probabilities of the collagen fiber at this point being aligned in that direction). Can you show how you would solve the problem with different intensities at different angles? $\endgroup$
    – user148316
    Feb 8 '17 at 22:38
  • 1
    $\begingroup$ I wondered whether that was the case. The fix is simple: multiply each (cos, sin) pair in your data by any weight you like. (The weights may even be negative--only their absolute values matter.) Then proceed as before. There's no need to divide by the sum of weights: any such division would cancel in the computation of the angle at the end. (An equivalent description is to note that weighted directional data are just complex numbers. Their mean direction is the argument (angle) of their arithmetic mean.) I have modified the code to accommodate arbitrary weights. $\endgroup$
    – whuber
    Feb 8 '17 at 23:03
  • $\begingroup$ Hi w huber, Follow up question here, now that I have a mean vector, what are your thoughts about finding the standard deviation? Everything that I have read has been an example of taking the standard deviation from unit circle data and not weighted data like I have. $\endgroup$
    – user148316
    Mar 1 '17 at 19:53

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.