/*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
cyclic.c - circular averages
by Andy Allinger, 2015, released to the public domain
This program may be used by any person for any purpose.
cymean() cyclic mean by a method due to David Radcliffe, 2005
cymed() cyclic median by an original variation of Radcliffe's algorithm
These functions are based on the FORTRAN subroutines CYMEAN and CYMED.
However, they use double precision and return PHASE ONLY!
For a full explanation, refer to the article:
"Circular Values Math and Statistics with FORTRAN" at www.codeproject.com
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
#include
#include
#include
#include
/* argument to qsort() */
static int cmpdbl (const void *a, const void *b)
{
double dif;
dif = *(double *)a - *(double *)b;
if (dif < 0.) {
return -1;
}
else if (dif > 0.) {
return 1;
}
else {
return 0;
}
} /* end function cmpdbl */
/*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
cymean() - cyclic mean
origin: David Radcliffe, 2005 (or earlier)
He proposed the algorithm:
1. Find the average and standard deviation of the given numbers.
2. Increase the smallest number by 360.
3. Repeat steps 1 and 2 until all numbers are greater than 360.
4. Choose the average that yields the smallest standard deviation.
5. If the average is greater than 360 then subtract 360.
Notice that unlike the linear mean, the cyclic mean does not
always have a distinct value. Instead, the array ave[] contains in
its beginning entries the various averages, and the integer ties denotes
how many averages were found.
WARNING: X is returned shifted and sorted
__Name____Type________________In/Out___Description______________________
*x double array In data
n int In length of x
t double In period
*ave double array Out cyclic mean
*ties pointer to int Out how many values in ave[]
*dev pointer to double Out standard deviation
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
void cymean (double *x, int n, double t, double *ave, int *ties, double *dev)
{
int i;
double r, s, v, s2, t2, nn, tol, vmin;
if (1 == n) {
*ties = 1;
ave[0] = x[0];
*dev = 0.;
return;
}
nn = (double) (n);
tol = sqrt(t * DBL_EPSILON);
t2 = t * t;
s = 0.;
s2 = 0.;
/* move all data to the range 0...t, sort */
for (i = 0; i < n; ++i) {
r = fmod(x[i], t);
if (r < 0.) r += t;
s += r;
s2 += r * r;
x[i] = r;
}
(void) qsort(x, n, sizeof(double), cmpdbl);
/* find phase of seam of lowest deviation */
vmin = DBL_MAX;
*ties = 0;
for (i = 0; i < n; ++i) {
v = s2 - s * s / nn; // n @ variance
if (vmin - v > tol) { // new low deviation
*ties = 0;
vmin = v;
}
if (abs(v - vmin) < tol) { // tie for low deviation
++(*ties);
ave[*ties-1] = s / nn;
}
/* apply update formula to sum and sum of squares */
s += t;
s2 = s2 + 2. * x[i] * t + t2;
}
for (i = 0; i < *ties; ++i) {
if (ave[i] >= t) ave[i] -= t; // subtract if out of 0...t
}
*dev = sqrt(vmin / nn); // (n@variance/n)^(1/2)
return;
} /* end function cymean */
/*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
cymed() - the cyclic median
This variation of Radcliffe's algorithm computes the median instead.
Notice that unlike the linear median, the cyclic median does not
always have a distinct value. Instead, the array ave[] contains in
its beginning entries the various averages, and the integer ties denotes
how many averages were found.
WARNING: X is returned shifted and sorted
__Name_______Type________________In/Out____Description____________________
*x double array In data
n int In length of x
t double In period
*ave double array Out cyclic medians
*ties pointer to int Out how many values in ave[]
*dev pointer to double Out mean absolute deviation
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
void cymed (double *x, int n, double t, double *ave, int *ties, double *dev)
{
int i, j0, j1;
double r, s1, smin, tol;
bool odd;
if (1 == n) {
*ties = 1;
ave[0] = x[0];
*dev = 0.;
return;
}
j1 = n / 2;
odd = (1 == (n % 2)); // is n an odd number?
tol = t * DBL_EPSILON;
/* move all data to the range 0...t, sort */
for (i = 0; i < n; ++i) {
r = fmod(x[i], t);
if (r < 0.) r += t;
x[i] = r;
}
qsort(x, n, sizeof(double), cmpdbl);
/* initial value of cyclic median */
if (odd) {
r = x[j1];
} else {
j0 = j1 - 1;
r = (x[j0] + x[j1]) / 2.; // split the difference
}
/* sum of absolute deviations */
s1 = 0.;
for (i = 0; i < n; ++i) {
s1 += abs(x[i] - r);
}
/* find phase of seam of lowest deviation */
smin = DBL_MAX;
*ties = 0;
for (i = 0; i < n; ++i) {
if (smin - s1 > tol) { // new low deviation
*ties = 0;
smin = s1;
}
if (abs(s1 - smin) < tol) { // tie for low deviation
++(*ties);
if (odd) {
ave[*ties-1] = x[j1];
} else {
j0 = j1 - 1; // j0 is the prior point
if (j0 < 0) j0 = n-1;
ave[*ties-1] = (x[j0] + x[j1]) / 2.;
}
}
/* apply update formula */
r = x[i];
if (odd) {
j0 = j1 + 1; // j0 is the new median
if (j0 >= n) j0 = 0;
s1 = s1 + 2. * r + t - x[j0] - x[j1];
} else {
s1 = s1 + 2. * r + t - 2. * x[j1];
}
x[i] = r + t; // advance seam
++j1;
if (j1 >= n) j1 = 0;
}
for (i = 0; i < *ties; ++i) {
if (ave[i] >= t) ave[i] -= t; // subtract if out of 0...t
}
*dev = smin / (double) (n); // mean absolute deviation
} /* end function cymed */
/*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
test program to replicate examples from Lior Kogan's codeproject.com article
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
#include
int main()
{
double x1[4] = {0., 30., 60., 90.}; // mean = 45
double x2[4] = {0., 90., 180., 270.}; // mean = {45, 135, 225, 315}
double x3[4] = {30., 130., 230., 330.}; // mean = 0
double x4[3] = {100., 200., 300.}; // median = 200
double x5[4] = {40., 100., 200., 300.}; // median = 70
double x6[4] = {0., 180., 270., 271.}; // median = 270.5 (some will argue)
double t, ave[4], dev;
int i, n, ties;
t = 360.;
n = 4;
printf("The set: {");
for (i=0; i