## 1. Introduction

Many scientific and engineering problems involve circular values. Most of the times these values are orientation (angles) or cyclical timing (time-of-day), but it may be any other circular quantity. Such problems are very common in physics, geodesics and navigation, but also appear in other fields such as psychology, criminology (time-of-day statistics), bird-watching and biology (directional statistics and time-of-year statistics).

In contrast to linear values (values on the real line), circular values have a limited range, and a wrap-around property.

Mathematics and Statistics of circular values is very tricky and error prone, both for simple operations such as addition and subtraction, and for more complex operators such as calculating average and median, estimations and interpolations.

For angles, it is even more error prone, since sometimes the range is used, sometimes . Sometimes it is and sometimes - even in the same application. For time-of-day computations, the range may be or something else. We need a scheme for using several representations in the same application, and to safely perform operations between them.

## 2. Features

In this article I am going to present:

- Definition of circular value type and its associated operators
- A C++11 infrastructure for doing mathematics with circular values
- C++11 statistical distribution classes for circular values: Wrapped Normal and Wrapped Truncated Normal
- Average, Weighted average and Median of circular values
- Circular parameter estimation based on noisy independent measurements
- Interpolation and average estimation of sampled continuous-time circular signal

This is an original work. To the best of my knowledge, some of the ideas mentioned here were not published before.

## 3. Requirements

To compile the code, you’ll need Visual C++ 2012.

The given code relies heavily of C++11 features.

The code may be easily converted to any other C++11 compiler.

The code may be converted to plain old C++. However, it will lose some of its charm.

## 4. Some helper functions

We’ll define the following helper functions:

- Square function: Sqr(r1). This is simple:
template <typename T>
T Sqr(const T& x)
{
return x*x;
}

- Floating-point modulo function:
*reminder* = *mod*(*dividend*, *divisor*)The *dividend*, *divisor* and the resulting *reminder* are all floating-point values.

The *quotient* and the *reminder* must satisfy:

- *quotient* is integer

- |*reminder*| < |*divisor*|

- *dividend* = *quotient* ⋅ *divisor* + *reminder*

Still, many definitions are possible:

A. The *reminder* has the same sign as the *dividend*

- Implies that the *quotient *rounds towards zero

B. The *reminder *has the same sign as the *divisor*

- Implies that the *quotient *rounds toward negative infinity

C. The *reminder *has a different sign than the *dividend*

- Implies that the *quotient *rounds away from zero

D. The *reminder *has a different sign than the *divisor*

- Implies that the *quotient *rounds toward positive infinity

E. The sign of the *reminder *is always non-negative

F. The sign of the *reminder *is always non-positive

G. The reminder is closest to zero

- Implies that the *quotient *is the integer nearest to the exact value of division's result

IEEE 754 names this *REM* operator, and disambiguate it for a boundary case.

When the fractional part of the division's result is exactly 0.5, the *quotient* is even.

H. The *reminder* is farthest from zero

Values for (*quotient*, *reminder*) according to the various implementations:

In addition, we should define the behavior when the *divisor *is 0:

X. When the *divisor *is 0, the *reminder *is undefined

Y. When the *divisor *is 0, the *reminder *is defined as 0

Z. When the *divisor *is 0, the *reminder *is defined as the *dividend*

Most programming languages, numerical computing environment, and spreadsheets, use definitions A, B or G. For example:

- C (ISO/IEC 9899) and C++ (ISO/IEC 14882) *fmod*() function implements the AX or AY definition

(the standard allow either)

- C (starting from ISO/IEC 9899-1999) and C++ (starting from ISO/IEC 14882-2011,

but not yet supported by many compiles) defines an additional function - *reminder*(),

which implements definition GX or GY (the standard allow either)

- Microsoft Excel's *MOD*() function implements the BX definition.

- MATLAB's *mod*() and *rem*() functions implement the BZ and AX definitions respectively.

####

We will use the following definition:

####

#### Definition 1

This definition was given by Knuth in *The Art of Computer Programming*, Vol.1 (p.39 of the 3rd edition), and is equivalent to definition BZ.

A straightforward implementation would be:

template<typename T>
T Mod(T x, T y)
{
static_assert(!std::numeric_limits<T>::is_exact , "Mod: floating-point type expected");
if (0 == y)
return x;
return x - y * floor(x/y);
}

However, such implementation is not resilient to boundary cases resulting from the noncontinuity of the modulo function and the floating-point limited accuracy. Here are two examples for double-precision values:

- Mod(-1e-16, 360.) = 360. (should be 0)

The value 360,-1e-16 = 359.9999999999999999 cannot be represented by double-precision

- Mod(106.81415022205296, _TWO_PI) = -1.421e-14 (should be _TWO_PI - 1.01e-14)

Therefore, we'll use the following implementation:

template<typename T>
T Mod(T x, T y)
{
static_assert(!std::numeric_limits<T>::is_exact , "Mod: floating-point type expected");
if (0 == y)
return x;
double m= x - y * floor(x/y);
if (y > 0)
{
if (m>=y)
return 0;
if (m<0 )
{
if (y+m == y)
return 0 ;
else
return y+m;
}
}
else
{
if (m<=y)
return 0;
if (m>0 )
{
if (y+m == y)
return 0 ;
else
return y+m;
}
}
return m;
}

#### Proof

- Almost-equal function: IsAlmostEqual(r1, r2)
When implementing unit-testing for floating-point values computations, it is a common need to verify that 2 computations give equal results. Due to the finite-accuracy of floating-point types, many times the values won’t be equal, but ‘almost equal’. For example `double d= 300.;`

is only almost-equal to `double e= exp(log(300.));`

```
```

To test for almost-equality, we need to consider not the absolute magnitude of the error, but its magnitude relative to the result. Moreover, we need to take care of marginal situations such as INF/NAN values.

Google Test is a framework for writing C++ tests. I’ve extracted the relevant functionality for testing almost equality into a single file - FPCompare.h, which I tweaked a little. Based on it, I implemented the following functions:

template<typename T>
static bool IsAlmostEq(T x, T y)
{
static_assert(!std::numeric_limits<T>::is_exact , "IsAlmostEq: floating-point type expected");
FloatingPoint<T> f(x);
FloatingPoint<T> g(y);
return f.AlmostEquals(g);
}
static void AssertAlmostEq(const double f, const double g)
{
assert(IsAlmostEq(f, g));
}

For more information, look at Google test at http://code.google.com/p/googletest/

## 5. A circular value type

#### Definition 2

A Circular value type is defined by three constants:

- The range
- The zero-value .

Note the use of a right-open interval range.

For example, for angles in the range , it is natural to define the zero-value in the middle, while for angles in the range it is natural to define zero-value on the edge. Generally, the zero-value may be any value in the range – not necessarily 0.

The zero-value has some special properties which will be defined below.

We will use the following macro to define a circular value type:

#define CircValTypeDef(_Name, _L, _H, _Z) \
struct _Name \
{ \
static const double L ; \
static const double H ; \
static const double Z ; \
static const double R ; \
static const double R_2; \
\
static_assert((_H>_L) && (_Z>=_L) && (_Z<_H), \
#_Name##": Range not valid"); \
}; \
\
const double _Name::L = (_L) ; \
const double _Name::H = (_H) ; \
const double _Name::Z = (_Z) ; \
const double _Name::R = ((_H)-(_L)) ; \
const double _Name::R_2= ((_H)-(_L))/2.;

And use it to define the following circular value types:

CircValTypeDef(SignedDegRange , -180., 180., 0. )
CircValTypeDef(UnsignedDegRange, 0., 360., 0. )
CircValTypeDef(SignedRadRange , -M_PI, M_PI, 0. )
CircValTypeDef(UnsignedRadRange, 0., 2*M_PI, 0. )

And some additional circular value types - for testing:

CircValTypeDef(TestRange0 , 3., 10., 5.3)
CircValTypeDef(TestRange1 , -3., 10., -3.0)
CircValTypeDef(TestRange2 , -3., 10., 9.9)
CircValTypeDef(TestRange3 , -13., -3., -5.3)

## 6. Defining a circular value class with a given type

The templated class `CircVal`

is used for storing a single *circular* value.

The template parameter should be defined using the `CircValTypeDef`

macro.

template <typename Type>
class CircVal
{
…
}

Here is an example of defining a circular value variable with a given type:

CircVal<UnsignedDegRange> c1;

## 7. Range checking and wrapping

`CircVal::IsInRange`

static function is used for testing that a given

*linear* value is within the range of a

`CircVal`

class.

static bool IsInRange(double r)
{
return (r>=Type::L && r<Type::H);
}

The *wrap*(*r*) function is used to ‘wrap around’ a *linear *value to the range of a given CircVal class.
For example: for the range , the value 360 would be wrapper to 0, and 370 would be wrapped to 10. -350 would be wrapped to 10 as well.

#### Definition 3

We will use this for faster implementation.

The `CircVal::Wrap`

static function calculates *wrap(r)*:

static double Wrap(double r)
{
if (r>=Type::L)
{
if (r< Type::H ) return r ;
else if (r< Type::H+Type::R) return r-Type::R;
}
else
if (r>=Type::L-Type::R) return r+Type::R;
return Mod(r - Type::L, Type::R) + Type::L;
}

## 8. A Walk

#### Definition 4

For *linear* values, a walk is defined by:

- A start point
- A directed length

A walk is a movement along the *linear* axis, from the start point, with the given directed length.

When the length is positive, we'll call it ‘an increasing walk’ (toward positive infinity).

When the length is negative, we'll call it ‘a decreasing walk’ (toward negative infinity).

The end point is the point reached at the end of the walk.

We’ll call a walk with a start point *s* and an end point *e* “A walk from *s* to *e*”.

There is only one possible walk from *s* to *e*: <*s*, *e*-*s* >

For any *circular* value type, a walk is defined by:

- A start point
- A directed length

A walk is a movement along the *circular* axis, from the start point, with the given directed length.

When the length is positive, we'll call it ‘an increasing walk’ (depicted clockwise in the picture above).

When the length is negative, we'll call it ‘a decreasing walk’ (depicted counterclockwise in the picture above).

The end point is the point reached at the end of the walk.

We’ll call a walk with a start point *s* and an end point *e *“A walk from *s* to *e*”.

Since the walk is wrapped around the circle:

There are infinitely many walks from :

The shortest walk from s to *e* is the one that minimizes |d|. The walk may be increasing or decreasing.

## 9. Distance between two values

#### Definition 5

The *directed distance* from *linear* value :

,

is the directed length of the walk from .

For example,

Similarly, the directed distance from *circular *value :

is the directed length of the shortest walk from :

It can easily be deduced that:

Since

And since

We can write:

which we will use for faster implementation.

The `CircVal::Sdist`

static function calculates *sdist(c1,c2)*:

static double Sdist(const CircVal& c1, const CircVal& c2)
{
double d= c2.val-c1.val;
if (d < -Type::R_2) return d + Type::R;
if (d >= Type::R_2) return d - Type::R;
return d ;
}

Note: To define a metric space, we can use |*sdist(c1,c2)*| as our metric (we won't prove this here).

For *circular* values, we'll define an additional distance function:

, The* increasing distance* from *circular *value , is the length of the shortest increasing-walk (depicted clockwise in the picture below) from .

It is always in the range .

There is no similar distance function for *linear *values.

#### Definition 6

which may also be expressed as:

It can easily be deduced that:

Hence,
it cannot be used as a metric

The `CircVal::Pdist`

static function calculates *pdist(c1,c2)*:

static double Pdist(const CircVal& c1, const CircVal& c2)
{
return c2.val>=c1.val ? c2.val-c1.val : Type::R-c1.val+c2.val;
}

#### Theorem 2

#### Proof

Rephrasing the definitions:

#### Theorem 3

(Very intuitive)

#### Proof

According to the definition above:

And also

, which are equal in all cases.

#### Theorem 4

For any two *circular* values:

#### Proof

#### Theorem 5

For any two circular values:

#### Proof

#### Theorem 6

For any two circular values:

#### Proof

Trivial, based on theorems 4, 5

## 10. Conversion between different types of circular values

We want to be able to perform operations between different types of *circular* values.

For example:

CircVal<UnsignedDegRange> d1= 10.;
CircVal<SignedDegRange > d2= -10.;
CircVal<SignedRadRange > d3 ;
d3= d1+d2;

To convert a circular value from circular type to circular type , we’ll retain the arc-length between and . For that, we can use either of these formulae:

#### Theorem 7

I and II are equivalent

#### Proof

Therefore, in both cases I =II

The conversion can be accomplished easily be implementing appropriate constructor and assignment operator:

template<typename CircVal2>
CircVal(const CircVal2& c2)
{
double val2= c2.Pdist(c2.GetZ(), c2);
val= Wrap(val2*Type::R/c2.GetR() + Type::Z);
}
template<typename CircVal2>
CircVal& operator= (const CircVal2& c2)
{
double val2= c2.Pdist(c2.GetZ(), c2);
val= Wrap(val2*Type::R/c2.GetR() + Type::Z);
return *this;
}

## 11. Constructing circular value / Fetching circular value

Construction a `CircVal`

objects based on *linear* value-representation is implemented by the appropriate constructor and assignment operator. Construction of a *linear* value based on a `CircVal`

object is implemented with `operator double()`

const overloading.

CircVal(double r)
{
val= Wrap(r);
}
CircVal& operator= (double r)
{
val= Wrap(r);
return *this;
}
operator double() const
{
return val;
}

## 12. Circular value operators – Intuitive description

First, here is an intuitive description for each operator (formal definition will be given later):

- Negation operator
*(-c)* is defined for any *circular* value. The result is a *circular* value.

Negation of a *circular* value is the symmetric value relative to *z* (the zero-value).

- Opposite operator
*(~c)* is defined for any *circular* value. The result is a *circular* value.

Opposite of a *circular* value is the *circular* value pointing to the opposite direction.

- Addition and subtraction operators are defined for any two
*circular* values.

The result is a *circular* value - the point reached and the end of the walk:

Addition: Starting from *c1*, walking with length *pdist*(*z,c*2).

Subtraction: Starting from *c1*, walking with length -*pdist*(*z,c*2).

- Multiplication operator
*(c·r)* is defined for any *circular* value *c* and any *linear* value *r*.

The result is a *circular* value - the point reached and the end of the walk:
Starting from *z*, walking distance *r*·*pdist*(*c,z*).

- Division operator (c/r) is defined for any
*circular* value *c* and any non-zero *linear* value *r*.

The result is a *circular* value - the point reached and the end of the walk:

Starting from *z*, walking times distance *pdist*(*c,z*).

- Trigonometric functions
*sin, cos, tan* are defined for any *circular* value *c*, such that the result will be a *linear* value, in accordance with the common functions usage.
- Inverse trigonometric functions
*asin, acose, atan, atan2* shall be defined for any *linear* number (-1 <= *r *< 1, for *asin, acos*), such that the result will be a *circular* value of our class in accordance with the common functions usage.

## 13. Conversion between circular values and linear values

We want a conversion function *ToR*(*c*), such that any *circular* value *c* in our range [*l,h*) can be mapped to a *linear* value, and a conversion function *ToR*(*r*), such that any *linear* value *r* can be mapped to a *circular* value in our in our range [*l,h*). Note that these functions are different from the *circular* value construction/fetching described above.

Requirements to these conversion functions:

We’ll use these conversion functions to simplify the requirements and the definitions of several *circular *value operators.

## 14. Requirements for circular value operators

We want to define the *circular* value operators
such that the following conditions are satisfied:

The
set of real numbers in the range [l,h) with the + binary
operator is a bounded group (*z* is
the identity element, and *-c* is
the inverse element of *c*).

## 15. Comparison operators

Defining operators = and

is straightforward. The other 4 comparison operators {<, >, =<, >=} can be defined in several different ways. Requirements to these comparison operators (For any

*circular* values

*c1*,

*c2*,

*c3*) should satisfy:

All these requirements, however, tell us nothing about the meaning of the > operator.

For linear values we can define:

Based on this, we can define the circular values operator in different ways:

<>Hence, for two circular values

*c*1,

*c*2, for each comparison operator

, we can use:

Though all definitions satisfy all the requirements, they mean different things. The 1^{st} compares the *pdist* from *l*, the 2^{nd} compares the *pdisht* from *z*, and the 3^{rd} compares the *sdist* from *z*.

In many practical situations *z* = 0.

For *C* = < -180,180,0 > it seems more appropriate to use (i) or (iii) which are equivalent, while for *C* = < 0,360,0 >, (i) and (ii) which are equivalent fits better. This is why our selected definition is (i), but it may be different for specific use cases.

bool operator==(const CircVal& c) const { return val == c.val; }
bool operator!=(const CircVal& c) const { return val != c.val; }
bool operator> (const CircVal& c) const { return val > c.val; }
bool operator>=(const CircVal& c) const { return val >= c.val; }
bool operator< (const CircVal& c) const { return val < c.val; }
bool operator<=(const CircVal& c) const { return val <= c.val; }

## 16. Circular value operators – Formal definition

Now, here is the definition that satisfies all the requirements

#### Definition 7

#### Theorem 8

The - unary operator can be equivalently defined by each of these two equivalent formulae:

The + and - binary operators can be equivalently defined by each of these two equivalent formulae:

#### Proof

Here is the implementation of these operators:

template <typename Type>
class CircVal
{
double val;
public:
friend double ToR(const CircVal& c) { return c.val - Type::Z; }
const CircVal operator+ ( ) const { return val; }
const CircVal operator- ( ) const { return Wrap(Type::Z-Sdist(Type::Z,val)); }
const CircVal operator~ ( ) const { return Wrap(val+Type::R_2 ); }
const CircVal operator+ (const CircVal& c) const { return Wrap(val+c.val - Type::Z); }
const CircVal operator- (const CircVal& c) const { return Wrap(val-c.val + Type::Z); }
const CircVal operator* (const double& r) const { return Wrap((val-Type::Z)*r + Type::Z); }
const CircVal operator/ (const double& r) const { return Wrap((val-Type::Z)/r + Type::Z); }
CircVal& operator+=(const CircVal& c) { val= Wrap(val+c.val - Type::Z); return *this; }
CircVal& operator-=(const CircVal& c) { val= Wrap(val-c.val + Type::Z); return *this; }
CircVal& operator*=(const double& r) { val= Wrap((val-Type::Z)*r + Type::Z); return *this; }
CircVal& operator/=(const double& r) { val= Wrap((val-Type::Z)/r + Type::Z); return *this; }
CircVal& operator =(const CircVal& c) { val= c.val ; return *this; }
…
}
template <typename Type> static double sin (const CircVal<Type>& c) { return std::sin(ToR(CircVal<SignedRadRange>(c))); }
template <typename Type> static double cos (const CircVal<Type>& c) { return std::cos(ToR(CircVal<SignedRadRange>(c))); }
template <typename Type> static double tan (const CircVal<Type>& c) { return std::tan(ToR(CircVal<SignedRadRange>(c))); }
template <typename Type> static CircVal<Type> asin (double r ) { return CircVal<SignedRadRange>(std::asin (r )); }
template <typename Type> static CircVal<Type> acos (double r ) { return CircVal<SignedRadRange>(std::acos (r )); }
template <typename Type> static CircVal<Type> atan (double r ) { return CircVal<SignedRadRange>(std::atan (r )); }
template <typename Type> static CircVal<Type> atan2(double r1, double r2 ) { return CircVal<SignedRadRange>(std::atan2(r1,r2)); }
template <typename Type> static CircVal<Type> ToC (double r ) { return CircVal<Type>::Wrap(r + Type::Z); }

Note that for `asin`

, `acos`

, `atan`

and `atan2`

functions a template parameter should be used. For example:

`CircVal<SignedDegRange> d1= asin<SignedDegRange>(0.5);`

`d1= asin(0.5)`

won’t give us the expected result unless our range is `SignedRadRange`

.

## 17. Testing for near-equality of circular values

Based on the previously describe test for near-equality for two

*linear* values, we will construct a near-equality test for two

*circular* values:

template <typename Type>
class CircValTester
{
static bool IsCircAlmostEq(const CircVal<Type>& _f, const CircVal<Type>& _g)
{
double f= _f;
double g= _g;
if (::IsAlmostEq(f, g))
return true;
if (f < g)
return IsAlmostEq(f, g - Type::R);
else
return IsAlmostEq(f, g + Type::R);
}
static void AssertCircAlmostEq(const CircVal<Type>& f, const CircVal<Type>& g)
{
assert(IsCircAlmostEq(f, g));
}
…
}

## 18. Testing correctness of circular value class implementation

Using test-driven design, we’re check that our implementation fulfills our requirements.

The Test() function generates 10,000 random test-cases, and verify that all the requirements holds.

We will use C++11 random generation, which is a very useful addition to the language:

template <typename Type>
class CircValTester
{
static void Test()
{
CircVal<Type> ZeroVal= Type::Z;
AssertCircAlmostEq(ZeroVal , -ZeroVal);
AssertAlmostEq (sin(ZeroVal) , 0. );
AssertAlmostEq (cos(ZeroVal) , 1. );
AssertAlmostEq (tan(ZeroVal) , 0. );
AssertCircAlmostEq(asin<Type>(0.), ZeroVal );
AssertCircAlmostEq(acos<Type>(1.), ZeroVal );
AssertCircAlmostEq(atan<Type>(0.), ZeroVal );
AssertCircAlmostEq(ToC<Type>(0) , ZeroVal );
AssertAlmostEq (ToR(ZeroVal) , 0. );
std::default_random_engine rand_engine ;
std::uniform_real_distribution<double> c_uni_dist(Type::L, Type::H);
std::uniform_real_distribution<double> r_uni_dist(0. , 1000. );
std::uniform_real_distribution<double> t_uni_dist(-1. , 1. );
std::random_device rnd_device;
rand_engine.seed(rnd_device());
for (unsigned i= 10000; i--;)
{
CircVal<Type> c1(c_uni_dist(rand_engine));
CircVal<Type> c2(c_uni_dist(rand_engine));
CircVal<Type> c3(c_uni_dist(rand_engine));
double r (r_uni_dist(rand_engine));
double a1(t_uni_dist(rand_engine));
double a2(t_uni_dist(rand_engine));
assert (c1 == CircVal<Type>((double)c1) );
AssertCircAlmostEq(+c1 , c1 );
AssertCircAlmostEq(-(-c1) , c1 );
AssertCircAlmostEq(c1 + c2 , c2 + c1 );
AssertCircAlmostEq(c1 + (c2 +c3) , (c1 + c2) + c3 );
AssertCircAlmostEq(c1 + -c1 , ZeroVal );
AssertCircAlmostEq(c1 + ZeroVal , c1 );
AssertCircAlmostEq(c1 - c1 , ZeroVal );
AssertCircAlmostEq(c1 - ZeroVal , c1 );
AssertCircAlmostEq(ZeroVal - c1 , -c1 );
AssertCircAlmostEq(c1 - c2 , -(c2 - c1) );
AssertCircAlmostEq(c1 * 0. , ZeroVal );
AssertCircAlmostEq(c1 * 1. , c1 );
AssertCircAlmostEq(c1 / 1. , c1 );
AssertCircAlmostEq((c1 * (1./(r+1.))) / (1./(r+1.)) , c1 );
AssertCircAlmostEq((c1 / ( r+1.) ) * ( r+1. ) , c1 );
AssertCircAlmostEq(~(~c1) , c1 );
AssertCircAlmostEq(c1 - (~c1) , ToC<Type>(Type::R/2.) );
AssertAlmostEq (sin(ToR(CircVal<SignedRadRange>(c1))), sin(c1) );
AssertAlmostEq (cos(ToR(CircVal<SignedRadRange>(c1))), cos(c1) );
AssertAlmostEq (tan(ToR(CircVal<SignedRadRange>(c1))), tan(c1) );
AssertAlmostEq (sin(-c1) , -sin(c1) );
AssertAlmostEq (cos(-c1) , cos(c1) );
AssertAlmostEq (tan(-c1) , -tan(c1) );
AssertAlmostEq (sin(c1+ToC<Type>(Type::R/4.)) , cos(c1) );
AssertAlmostEq (cos(c1+ToC<Type>(Type::R/4.)) , -sin(c1) );
AssertAlmostEq (sin(c1+ToC<Type>(Type::R/2.)) , -sin(c1) );
AssertAlmostEq (cos(c1+ToC<Type>(Type::R/2.)) , -cos(c1) );
AssertAlmostEq (Sqr(sin(c1))+Sqr(cos(c1)) , 1. );
AssertAlmostEq (sin(c1)/cos(c1) , tan(c1) );
AssertCircAlmostEq(asin<Type>(a1) , CircVal<SignedRadRange>(asin(a1)));
AssertCircAlmostEq(acos<Type>(a1) , CircVal<SignedRadRange>(acos(a1)));
AssertCircAlmostEq(atan<Type>(a2) , CircVal<SignedRadRange>(atan(a2)));
AssertCircAlmostEq(asin<Type>(a1) + asin<Type>(-a1) , ZeroVal );
AssertCircAlmostEq(acos<Type>(a1) + acos<Type>(-a1) , ToC<Type>(Type::R/2.) );
AssertCircAlmostEq(asin<Type>(a1) + acos<Type>( a1) , ToC<Type>(Type::R/4.) );
AssertCircAlmostEq(atan<Type>(a2) + atan<Type>(-a2) , ZeroVal );
assert (c1 > c2 == (c2 < c1) );
assert (c1 >= c2 == (c2 <= c1) );
assert (c1 >= c2 == ( (c1 > c2) || (c1 == c2)) );
assert (c1 <= c2 == ( (c1 < c2) || (c1 == c2)) );
assert (c1 > c2 == (!(c1 == c2) && !(c1 < c2)) );
assert (c1 == c2 == (!(c1 > c2) && !(c1 < c2)) );
assert (c1 < c2 == (!(c1 == c2) && !(c1 > c2)) );
assert (!(c1>c2) || !(c2>c3) || (c1>c3) );
AssertCircAlmostEq(c1 , ToC<Type>(ToR( c1) ) );
AssertCircAlmostEq(-c1 , ToC<Type>(ToR(-c1) ) );
AssertCircAlmostEq(c1 + c2 , ToC<Type>(ToR(c1)+ToR(c2)) );
AssertCircAlmostEq(c1 - c2 , ToC<Type>(ToR(c1)-ToR(c2)) );
AssertCircAlmostEq(c1 * r , ToC<Type>(ToR(c1)*r ) );
AssertCircAlmostEq(c1 / r , ToC<Type>(ToR(c1)/r ) );
}
}
public:
CircValTester()
{
Test();
}
};

The code below executes these tests for all 8 circular-value defined types:

CircValTester<SignedDegRange > testA;
CircValTester<UnsignedDegRange> testB;
CircValTester<SignedRadRange > testC;
CircValTester<UnsignedRadRange> testD;
CircValTester<TestRange0 > test0;
CircValTester<TestRange1 > test1;
CircValTester<TestRange2 > test2;
CircValTester<TestRange3 > test3;

## 19. Using the CircVal class

This sample code below speaks for itself:

CircVal<SignedDegRange > d1= 10. ;
CircVal<UnsignedRadRange> d2= 0.2;
CircVal<SignedDegRange > d3= d1+d2;
d1+= 355.;
double d= d1;
d = sin(d1) / cos(d2) + tan(d3);
d1= asin<SignedDegRange>(0.5);

The code is simple, yet powerful.

## 20. C++11 Distribution classes for circular values

One of the additions to the C++ standard is the <random> header. It is now possible to generate pseudo-random values with a given distribution. Many distributions classes are provided: Bernoulli, Binomial, Cauchy, Chi-Square, Exponential, Extreme Value, Fisher F, Gamma, Geometric, Log-Normal, Negative-Binomial, Normal, Poisson, Student T, Uniform, Weibull, and some interval-related distributions.

Circular distribution classes, however, are not implemented. Some Circular distributions are von Mises, Wrapped Normal and Wrapped Cauchy.

Another distribution, which is very useful when working with circular values, is the truncated normal distribution.

We have implemented three distribution classes: wrapped normal, truncated normal and wrapped truncated normal.

## 21. Wrapped Normal Distribution class

A wrapped normal distribution results from the "wrapping" of the normal distribution around the range .
Therefore, the probability density function of the wrapped normal distribution is:

The PDF is described in the following image, for same , but different values of :

This is the distribution of the modulo of a normally-distributed value.

Example: Assume that we have clocks (date + time-of-day), which have a normal-distributed error. It’s time-of-day [00:00, 24:00) distribution would be wrapped-normal.

Wrapped normal distributed values can be generated as follows:

#include "WrappedNormalDist.h"
std::default_random_engine rand_engine;
std::random_device rnd_device ;
rand_engine.seed(rnd_device());
double fAvrg = 0.;
double fSigma= 45.;
double fL = -180.;
double fH = 180.;
wrapped_normal_distribution<double> r_wrp(fAvrg, fSigma, fL, fH);
double r1= r_wrp(rand_engine);

The implementation is based on VC 2012 implementation of the normal distribution, added wrapping.

## 22. Truncated Normal Distribution class

This class has nothing to do with circular values; however it is a base for the wrapped truncated normal distribution class described in the next section.

A truncated normal distribution is the probability distribution of a normally distributed random variable whose value is bounded to the range . The PDF is described in the following image, for different values of :

Many phenomena behave, or modeled by a truncated normal distribution. One type of examples is the estimation of quantities by “ordinary” humans, such as today’s temperature: if it is 25° Celsius, no “ordinary” human will estimate it below 10°, or above 40°. Another example is the life span of many species (based on data from Shiro Horiuchi - Interspecies Differences in the Life Span Distribution: Humans versus Invertebrates, The Population Council 2003).

Truncated normal distribution is also useful for generating pseudo-random values such as the described above.

Truncated normal distributed values can be generated as follows:

#include "TruncNormalDist.h"
std::default_random_engine rand_engine;
std::random_device rnd_device ;
rand_engine.seed(rnd_device());
double fAvrg = 0.;
double fSigma= 45.;
double fA = -40.;
double fB = 40.;
truncated_normal_distribution<double> r_trn(fAvrg, fSigma, fA, fB);
double r2= r_trn(rand_engine);

The implementation is based on VC 2012 implementation of the normal distribution as a skeleton, and on C. H. Jackson's R's implementation of the following paper: “Robert, C. P. - Simulation of truncated normal variables. Statistics and Computing (1995) 5, 121–125”

## 23. Wrapped Truncated Normal Distribution class

Same as for linear values, circular values wrapped normal distribution may be truncated as well.
Note that the normal distribution is first truncated to the range , and then wrapped around the range

Many phenomena behave, or modeled by a wrapped truncated normal distribution. One type of examples is the estimation of circular quantities by “ordinary” humans, such as the time of sunrise tomorrow, or the azimuth between two drawn points on a piece of paper: If it is 45°, no “ordinary” human will estimate it below 10° or above 80°. Another example may be the exact time-of-day at which John Doe wakes up on an “ordinary” working day.

Wrapped truncated normal distribution is also useful for generating pseudo-random values such as the described above.

Usually, the truncation range is such that *b-a < h-l*, however, this is not a necessity. For example, we may model the time-of-day error distribution of a given clock that was set 1 year ago, and its absolute error now is modeled by a truncated normal distribution of *[-50,50)* hours.

After wrapping, it is possible that *wrap(b) < wrap(a)* - See picture below, which may depict human estimation of an azimuth between two points, where the real azimuth is 0°.

Wrapped truncated normal distributed values can be generated as follows:

#include "WrappedTruncNormalDist.h"
std::default_random_engine rand_engine;
std::random_device rnd_device ;
rand_engine.seed(rnd_device());
double fAvrg = 0.;
double fSigma= 100.;
double fA = -500.;
double fB = 500.;
double fL = 0.;
double fH = 360.;
wrapped_truncated_normal_distribution<double> r_ wrp_trn(fAvrg, fSigma, fA, fB, fL, fH);
double d= r_wrp_trn(rand_engine);

## 24. 3 types of statistical problems

There are 3 distinct types of problems, which tends to be confused one with the other, especially when it comes to circular values. Although these problems may seem similar, the mathematical methods required to solve each of them should be considered separately.

For each problem type, we’ll give a description, and sample problem for *linear* values and for *circular* values.

### Problem Type 1

Given [circular] values *x*_{1}..*x*_{N} – calculate their mean

-or-

Given a random sample *x*_{1}..*x*_{N} from a [circular] random-variable *X* - calculate the sample mean

Usage: Averaging values

Examples:

*Linear*: Given the number of births occurred in the US for each day in the year 2000 -

Calculate the mean number of births per day
* Circular*: Given time-of-day [00:00-24:00) for each birth occurred in US in the year 2000 -

Calculate the mean time-of-day

### Problem Type 2

Given a multiset of measurements/observations with a [wrapped-/wrapped-truncated-]normal-distributed error of a [circular] parameter – Estimate the parameter

Usage: Estimation of an unknown constant value, based on a multiset of noisy measurements, where the noise has a [wrapped-/wrapped-truncated-]normal-distribution

Examples:

*Linear*: Given a multiset of distance measurements from a stationary transmitter to a stationary receiver, using a measurement technique with a normal distributed error – Estimate the distance.
* Circular*: Given a multiset of direction measurements from a stationary transmitter to a stationary receiver, using a measurement technique with a wrapped normal distributed error – Estimate the direction.
* Linear*: Given a multiset of distance estimates between two points, made by “ordinary” humans (assuming to subject to truncated normal distributed error) - Estimate the distance.
*Circular*: Given a multiset of azimuth estimates between two points, made by “ordinary” humans (assuming to subject to a wrapped truncated normal distributed error) – Estimate the direction.

### Problem Type 3

Given a multiset of samples x_{1}..x_{N }of a continuous-time [circular] signal, sampled at times t_{1}..t_{N} respectively, where the measurement technique/instrument is accurate - Estimate the mean value of the signal at period [t_{1}, t_{N}]

Usage: Estimating the average of a continuous-signal, in a given period, based on a set of time-tagged samples

Examples:

* Linear*: Given airplane velocity measured each 1 second [MPH] over a period of one hour – Estimate the mean velocity.
* Circular*: Given airplane heading measured each 1 second over a period of one hour - Estimate the mean heading.

## 25. Averaging 2 circular values

To average two *linear* values, , we can start at and take a walk with length . The end-point of the walk is the average. We can, by symmetry, start at and take a walk with length .

For *circular* values, given the range [0,360), the average of 330 and 30 is {0} and not {(330+30)/2} = {180}.
The average is in the middle of the walk with the lowest absolute-value length, from 330 to 30 (or from 30 to 330).

The average of *circular* values is a not a single *circular* value, but a set of *circular* values.

*average*(0,180) = {90, 270} - There is no reason to prefer one of these values over the other.

We will use the following formula:

## 26. Averaging n circular values

*Circular *values average calculation is really tricky. Observing algorithms described in many textbooks, and code written by many programmers, I've almost never seen a correct method used.

Let us start with *linear* values. Here is our reference problem:

- Given the number of births occurred in the US for each day in the year 2000 -

Calculate the mean number of births per day

The solution here is quite simple: Just sum the number of births, and divide it by the number of days.

Why does it work?

The general definition of the mean (arithmetic-average) of values for any mathematical field *F*, is the value of x∈*F* that minimizes . Formally: .

For *linear *values, we can calculate the average as follow:

Let *y* denote the expression we want to minimize.

To find the value of *x* that minimizes *y*, we need the derivative to be 0:

We got the well-known arithmetic mean formula for *linear *values.

Now let’s talk about *circular *values. Here is our reference problem:

- Given time-of-day [00:00-24:00) for each birth occurred in US in the year 2000 -
Calculate the mean time-of-day

For making the following explanation simple, let us consider the range [0,360).

Note that the true average of 0, 0, and 90 should be 30, and not atan , as suggested in most references.

Based on theorem 3, the absolute value of the distance between 2 *circular *values in the range [0,360) is

This formula is suitable, since we want to minimize the difference between a value and the average, regardless of which ‘direction’ of the average the value is.

Let *y* denote the expression we want to minimize.

The values of *x* that minimizes *y* is our average.

Here are some examples:

Example 1: For the input {0,30,60,90}, we’ll plot as a function of *x*:

We can see that the average is {45} which matches out intuition.

Example 2: For the input {0,90,180,270}, we’ll plot as a function of *x*:

In this case the average is {45.135, 225, 315}

(For all these values, is equal).

Example 3: For the input {30,130,230,330}, we’ll plot as a function of *x*:

In this case, the average is {0}. Again, it matches our intuition.

The following code was used to collect the data for the above graph:

vector<CircVal<UnsignedDegRange>> Angles2;
Angles2.push_back(CircVal<UnsignedDegRange>( 30.));
Angles2.push_back(CircVal<UnsignedDegRange>(130.));
Angles2.push_back(CircVal<UnsignedDegRange>(230.));
Angles2.push_back(CircVal<UnsignedDegRange>(330.));
auto y= CircAverage(Angles2);
ofstream f0("log0.txt");
for (double x= 0.; x<=360.; x+= 0.1)
{
double fSum= 0;
for (const auto& a : Angles2)
fSum+= Sqr(__min(abs(x-a), 360.-abs(x)));
f0 << x << "\t" << fSum << endl;
}

To find the values of that *x* minimizes *y*, we need the derivative to be 0.

However, due to the *min* and *abs* operators, the expression is not derivable.

Let us try to rephrase the expression.

We’ll divide the multiset into 3 distinct multi-subsets:

We’ll denote with *a, b, c, d* – the number of elements in *A, B, C, D* respectively.

Of course

The minimum of *y* may be where , or where *y* is not derivable.

Since in each sector *y*(*x*) is a 2^{nd} order polynomial with a positive coefficient for the x² term, it has a single minima, which is when . In the sample graphs above, it is easy to see the sectors - each sector is a trimmed 2^{nd} order curve. The points where *y*(*x*) is not derivable, is on the maximum of each sector.

So, by dividing the input into 3 multisets, we can find the values of *x* that minimizes *y*. This is our average.

However, the division of the input elements into 3 multisets is dependent on the value of *x*.

Here is a graphic visualization of these multisets, for different values of *x*:

If *x* = 180, all values are in multiset *B*.

If *x* is in the range [0,180), multiset *C* range is (*x* + 180,360) and multiset *D* is empty.

If *x* is in the range (180,360), multiset *D* range is [0, *x* - 180) and multiset *C* is empty.

Multisets *C* and *D* cannot coexist.

Here is the average calculation algorithm:

- Divide the circle into sectors, such that if is within a sector - the elements of multisets
*B*, *C* and *D* does not change. The maximal number of sectors equals to the number of elements in *A*.
- For each sector: Assume that (the value of x that minimizes y) is within this sector, calculate its value, and check that the value is indeed within this sector. If so, calculate
*y* as follows:
- Assuming x=180, all values are in multiset
*B*:

- Assuming
*x* is in the range [0,180), multiset *C* range is (x+180,360) and multiset *D* is empty:

- Assuming
*x* is in the range (180,360), multiset *D* range is [0,x-180) and multiset *C* is empty:

- Find the sectors(s) that has the lowest
*y*; return their *x*.

The Function `CircAverage`

calculates the average set of a given vector of circular values.

The return value is a set of all average values:

template<typename T>
set<CircVal<T>> CircAverage(vector<CircVal<T>> const& A)
{
set<CircVal<T>> MinAvrgVals ;
double fSum = 0.;
double fSumSqr = 0.;
double fMinSumSqrDiff ;
vector<double> LowerAngles ;
vector<double> UpperAngles ;
double fTestAvrg ;
auto SumSqr = [&]() -> double
{
return 32400.*A.size() - 360.*fSum + fSumSqr;
};
auto SumSqrC= [&](double x, size_t nCountC, double fSumC) -> double
{
return x*(A.size()*x - 2*fSum) + fSumSqr - 2*360.*fSumC + nCountC*( 2*360.*x + 360.*360.);
};
auto SumSqrD= [&](double x, size_t nCountD, double fSumD) -> double
{
return x*(A.size()*x - 2*fSum) + fSumSqr + 2*360.*fSumD + nCountD*(-2*360.*x + 360.*360.);
};
auto TestSum= [&](double fTestAvrg, double fTestSumDiffSqr) -> void
{
if (fTestSumDiffSqr < fMinSumSqrDiff)
{
MinAvrgVals.clear();
MinAvrgVals.insert(CircVal<UnsignedDegRange>(fTestAvrg));
fMinSumSqrDiff= fTestSumDiffSqr;
}
else if (fTestSumDiffSqr == fMinSumSqrDiff)
MinAvrgVals.insert(CircVal<UnsignedDegRange>(fTestAvrg));
};
for (const auto& a : A)
{
double v= CircVal<UnsignedDegRange>(a);
fSum += v ;
fSumSqr+= Sqr(v);
if (v < 180.) LowerAngles.push_back(v);
else if (v > 180.) UpperAngles.push_back(v);
}
sort(LowerAngles.begin(), LowerAngles.end() );
sort(UpperAngles.begin(), UpperAngles.end(), greater<double>());
MinAvrgVals.clear();
MinAvrgVals.insert(CircVal<UnsignedDegRange>(180.));
fMinSumSqrDiff= SumSqr();
double fLowerBound= 0.;
double fSumD = 0.;
auto iter= LowerAngles.begin();
for (size_t d= 0; d < LowerAngles.size(); ++d)
{
fTestAvrg= (fSum + 360.*d)/A.size();
if ((fTestAvrg > fLowerBound+180.) && (fTestAvrg <= *iter+180.))
TestSum(fTestAvrg, SumSqrD(fTestAvrg, d, fSumD));
fLowerBound= *iter ;
fSumD += fLowerBound;
++iter;
}
fTestAvrg= (fSum + 360.*LowerAngles.size())/A.size();
if ((fTestAvrg < 360.) && (fTestAvrg > fLowerBound))
TestSum(fTestAvrg, SumSqrD(fTestAvrg, LowerAngles.size(), fSumD));
double fUpperBound= 360.;
double fSumC = 0.;
iter= UpperAngles.begin();
for (size_t c= 0; c < UpperAngles.size(); ++c)
{
fTestAvrg= (fSum - 360.*c)/A.size();
if ((fTestAvrg >= *iter-180.) && (fTestAvrg < fUpperBound-180.))
TestSum(fTestAvrg, SumSqrC(fTestAvrg, c, fSumC));
fUpperBound= *iter ;
fSumC += fUpperBound;
++iter;
}
fTestAvrg= (fSum - 360.*UpperAngles.size())/A.size();
if ((fTestAvrg >= 0.) && (fTestAvrg < fUpperBound))
TestSum(fTestAvrg, SumSqrC(fTestAvrg, UpperAngles.size(), fSumC));
return MinAvrgVals;
}

Note that the complexity of this function is the complexity of the sort operation, hence *O(a log a)*. Given sorted input, the complexity can be *O(a)*.

Here is a sample use:

vector<CircVal<UnsignedDegRange>> angles;
angles.push_back( 90.);
angles.push_back(180.);
angles.push_back(270.);
auto avg= CircAverage(angles);

## 27. Weighted Average of circular values

The general definition of the weighted mean (arithmetic-average) of values with weights respectively, for any mathematical field is the value *x* that minimizes the expression .

Formally: .

*w*_{i} Is a positive real weight associated with value *a*_{i}, and *dist*(*v*1,*v*2) is a function that measures the minimal distance between *v*1 and *v*2.

For linear values, we can calculate the average as follow:

Let *y* denote the expression we want to minimize:

To find the value of *x* that minimizes *y*, we need the derivative to be 0:

We got the well-known weighted arithmetic mean formula for *linear* values.

For *circular* values, we showed that

Let *y* denote the expression we want to minimize:

Same as for non-weighted average, we will break the input into 3 distinct multi-subsets: *B*, *C* and *D*.

And the algorithm for weighted average is:

- Divide the circle into sectors, such that if is
*x* within a sector - the elements of multisets *B*, *C* and *D* does not change. The maximal number of sectors equals to the number of elements we want to average.
- For each sector: Assume that (the value of
*x* that minimizes *y*) is within this sector, calculate its value, and check that the value is indeed within this sector. If so, calculate *y* as follows:
- Assuming x=180, all values are in multiset
*B*:

- Assuming
*x* is in the range [0,180), multiset *C* range is (x+180,360) and multiset *D* is empty:

- Assuming
*x* is in the range (180,360), multiset *D* range is [0,x-180) and multiset *C* is empty:

- Find the sectors(s) that has the lowest
*y*; return their *x*.

The Function `WeightedCircAverage`

calculates the weighted average set of a given vector of *circular* values.

The return value is a set of all average values:

template<typename T>
set<CircVal<T>> WeightedCircAverage(vector<pair<CircVal<T>,double>> const& A)
{
set<CircVal<T>> MinAvrgVals ;
double fASumW = 0.;
double fASumWA = 0.;
double fASumWA2 = 0.;
double fMinSumSqrDiff ;
vector<pair<double, double>> LowerAngles ;
vector<pair<double, double>> UpperAngles ;
double fTestAvrg ;
auto SumSqr = [&]() -> double
{
return 32400.*fASumW - 360.*fASumWA + fASumWA2;
};
auto SumSqrC= [&](double x ,
double fCSumW ,
double fCSumWC ) -> double
{
return fASumWA2 + x*x*fASumW -2*x*fASumWA - 720*fCSumWC + (129600+720*x)*fCSumW;
};
auto SumSqrD= [&](double x ,
double fDSumW ,
double fDSumWD ) -> double
{
return fASumWA2 + x*x*fASumW -2*x*fASumWA + 720*fDSumWD + (129600-720*x)*fDSumW;
};
auto TestSum= [&](double fTestAvrg, double fTestSumDiffSqr) -> void
{
if (fTestSumDiffSqr < fMinSumSqrDiff)
{
MinAvrgVals.clear();
MinAvrgVals.insert(CircVal<UnsignedDegRange>(fTestAvrg));
fMinSumSqrDiff= fTestSumDiffSqr;
}
else if (fTestSumDiffSqr == fMinSumSqrDiff)
MinAvrgVals.insert(CircVal<UnsignedDegRange>(fTestAvrg));
};
for (const auto& a : A)
{
double v= CircVal<UnsignedDegRange>(a.first);
double w= a.second;
fASumW += w ;
fASumWA+= w*v ;
fASumWA2= w*v*v;
if (v < 180.) LowerAngles.push_back(pair<double,double>(v,w));
else if (v > 180.) UpperAngles.push_back(pair<double,double>(v,w));
}
sort(LowerAngles.begin(), LowerAngles.end() );
sort(UpperAngles.begin(), UpperAngles.end(), greater<pair<double,double>>());
MinAvrgVals.clear();
MinAvrgVals.insert(CircVal<UnsignedDegRange>(180.));
fMinSumSqrDiff= SumSqr();
double fLowerBound= 0.;
double fDSumW = 0.;
double fDSumWD = 0.;
auto iter= LowerAngles.begin();
for (size_t d= 0; d < LowerAngles.size(); ++d)
{
fTestAvrg= (fASumWA + 360.*fDSumW)/fASumW;
if ((fTestAvrg > fLowerBound+180.) && (fTestAvrg <= (*iter).first+180.))
TestSum(fTestAvrg, SumSqrD(fTestAvrg, fDSumW, fDSumWD));
fLowerBound= (*iter).first ;
fDSumW += (*iter).second ;
fDSumWD += (*iter).second * (*iter).first;
++iter;
}
fTestAvrg= (fASumWA + 360.*fDSumW)/fASumW;
if ((fTestAvrg < 360.) && (fTestAvrg > fLowerBound))
TestSum(fTestAvrg, SumSqrD(fTestAvrg, fDSumW, fDSumWD));
double fUpperBound= 360.;
double fCSumW = 0.;
double fCSumWC = 0.;
iter= UpperAngles.begin();
for (size_t c= 0; c < UpperAngles.size(); ++c)
{
fTestAvrg= (fASumWA - 360.*fCSumW)/fASumW;
if ((fTestAvrg >= (*iter).first-180.) && (fTestAvrg < fUpperBound-180.))
TestSum(fTestAvrg, SumSqrC(fTestAvrg, fCSumW, fCSumWC));
fUpperBound= (*iter).first ;
fCSumW += (*iter).second ;
fCSumWC += (*iter).second * (*iter).first;
++iter;
}
fTestAvrg= (fASumWA - 360.*fCSumW)/fASumW;
if ((fTestAvrg >= 0.) && (fTestAvrg < fUpperBound))
TestSum(fTestAvrg, SumSqrC(fTestAvrg, fCSumW, fCSumWC));
return MinAvrgVals;
}

Note that the complexity of this function is the complexity of the sort operation, hence

*O(a log a)*. Given sorted input, the complexity can be

*O(a)*.
Here is a sample use:

vector<pair<CircVal<UnsignedDegRange>,double>> angles;
angles.push_back(make_pair( 90., 0.3));
angles.push_back(make_pair(180., 0.5));
angles.push_back(make_pair(270., 0.7));
auto avg= WeightedCircAverage(angles);

## 28. Median of circular values

Definition 7

Given the multiset :

Let sequence be the non-decreasingly sorted permutation of

For *linear* values:

If *n* is odd:

If *n* is even:

For *circular* values:

Since there is no "first" nor "last" element, we'll evaluate all rotations of *S*:

For each k∈{1…n}:

If *n* is odd:

If *n* is even:

(The avrg set may contain more than one element)

The set of all candidate medians:

Since an optimality property of the median is the minimization of the sum of absolute distances:

(The median set may contain more than one element)

Here are some examples for *circular *values, given the range [0,360):

The function `CircMedian`

calculates the median set of a given vector of *circular* values.

The return value is a set of all median values:

template<typename T>
set<CircVal<T>> CircMedian(vector<CircVal<T>> const& A)
{
set <CircVal<T>> X;
set<CircVal<T>> B;
if (A.size() % 2 == 0)
{
auto S= A;
sort(S.begin(), S.end());
for (size_t m= 0; m < S.size(); ++m)
{
size_t n= m+1; if (n==S.size()) n= 0;
double d= CircVal<T>::Sdist(S[m],S[n]);
B.insert(((double)S[m] + d/2.));
if (d == -CircVal<T>::GetR()/2.)
B.insert(((double)S[n] + d/2.));
}
}
else
for (size_t m= 0; m < A.size(); ++m)
B.insert(A[m]);
double fMinSum= numeric_limits<double>::max();
for (const auto& b : B)
{
size_t nL = 0 ;
size_t nH = 0 ;
double fSum= 0.;
for (const auto& a : A)
{
double d= CircVal<T>::Sdist(b,a);
if (d == 0.) continue ;
if (d < 0.) ++nL; else ++nH;
fSum+= abs(d);
}
if ((nL<= A.size()/2.) && (nH<= A.size()/2.))
if (fSum==fMinSum) X.insert(b);
else if (fSum< fMinSum) { X.clear(); X.insert(b); fMinSum= fSum; }
}
return X;
}

Here is a sample use:

vector<CircVal<UnsignedDegRange>> angles;
angles.push_back( 90.);
angles.push_back(180.);
angles.push_back(270.);
auto mdn= CircMedian(angles);

## 29. Circular parameter estimation based on noisy measurements

For *linear* values:

For many observations/measurements techniques/instruments, normal (Gaussian) distribution is often used as a first approximation to model observations/measurements error. As Jonas Ferdinand Gabriel Lippmann said: “Everyone believes in the [normal] law of errors: the mathematicians, because they think it is an experimental fact; and the experimenters, because they suppose it is a theorem of mathematics.”
Hence, for our reference problem:

- Given a multiset of independent distance measurements from a stationary transmitter to a stationary receiver, using a measurement technique with a normal-distributed error – Estimate the distance.

It can be shown that given a multiset of independent measurements of a constant parameter, were the measurement technique/instrument has a normal distributed error, the maximum likelihood estimation (MLE) of the measured parameters equals to the arithmetic average of the measurements.

Hence, for our reference problem, the maximum likelihood estimation of the distance is the arithmetic average of the measurements.

For *circular* value:

Here, wrapped normal distribution, or wrapped truncated normal distribution may be used to model a measurement error (according to the physical phenomena we want to model). The maximal measurement error is bounded to half circle.

Given a multiset of measurements of a constant *circular* parameter, were the measurement technique/instrument has a wrapped/wrapped-truncated normal distributed error.

Here is our reference problem:

- Given a multiset of independent direction measurements from a stationary transmitter to a stationary receiver, using a measurement technique with a wrapped/wrapped-truncated normal-distributed error – Estimate the direction.

We will consider two methods to estimate the

*circular* parameter:

- The
*circular* average of the samples - as described above
- Consider each sample as a 2D unit vector, and calculate the weighted vector - The suggested method in most textbooks

Simulation: As a test case, we will generate 50,000 trails. Each trail will contain 1000 samples of measurements with wrapped/wrapped-truncated normal distributed error (90° truncation span). We will use the two methods to calculate the estimate for each trail, and then calculate the RMS error for each method.

We will repeat the test for different standard deviations – between 1° and 100°.

The following code will collect the data for the two methods:

std::default_random_engine rand_engine;
std::random_device rnd_device ;
rand_engine.seed(rnd_device());
ofstream f1("log1.txt");
Concurrency::parallel_for(1, 101, [&](int nStdDev)
{
uniform_real<double> ud(0., 360.);
double fSumSqrErr1= 0.;
double fSumSqrErr2= 0.;
const size_t nTrails = 50000;
const size_t nSamples= 1000;
vector<CircVal<UnsignedDegRange>> vInput(nSamples);
const double fAvrg= ud(rand_engine);
wrapped_normal_distribution <double> r_wnd1(fAvrg, nStdDev, 0., 360.);
for (size_t t= 0; t<nTrails; ++t)
{
for (auto& Sample : vInput)
Sample= r_wnd1(rand_engine);
set<CircVal<UnsignedDegRange>> sAvrg1= CircAverage(vInput);
double fSigSin= 0.;
double fSigCos= 0.;
for (const auto& Sample : vInput)
{
fSigSin+= sin(Sample);
fSigCos+= cos(Sample);
}
CircVal<UnsignedDegRange> Avrg2= atan2<UnsignedDegRange>(fSigSin, fSigCos);
const double fErr1= CircVal<UnsignedDegRange>::Sdist(*sAvrg1.begin(), fAvrg);
const double fErr2= CircVal<UnsignedDegRange>::Sdist(Avrg2 , fAvrg);
fSumSqrErr1+= Sqr(fErr1);
fSumSqrErr2+= Sqr(fErr2);
}
const double fRMS1= sqrt(fSumSqrErr1/ (nTrails-1));
const double fRMS2= sqrt(fSumSqrErr2/ (nTrails-1));
f1 << nStdDev << "\t" << fRMS1 << "\t" << fRMS2 << endl;
} );

Note that `parallel_for`

is used for parallelizing tests for different standard deviations.

So which method is better?

The graph above shows the RMS error of the average estimation for 50000 trails. Each trail averages 1000 measurements with **wrapped normal** distributed error with standard deviation according to the X axis.

Method 1 turns out to give more accurate estimations when the standard deviation is less than 65°-80° (this value grows when there are more samples per trail: 65° for 10 samples per trail, and 78° for 1000 samples per trail). If the data is noisier - method 2 becomes better. Since the standard deviation in most real scenarios is much less than 65°, method 1 is better for almost any practical purpose.

The graph above shows the RMS error of the average estimation for 50000 trails. Each trail averages 1000 measurements with **wrapped truncated normal** distributed error with standard deviation according to the X axis.

The truncation span is 90° around the parameter value, and the wrapping range is [0,360).

Method 1 proves to be better for wrapped truncated normal standard deviation error, for any value of standard deviation.

## 30. Interpolation and average estimation of sampled continuous-time circular signal

Given a sequence of samples of a continuous-time *circular* signal, sampled at times respectively, where the measurement technique/instrument is accurate - Estimate the mean value of the signal at period .

*Linear* example:

- Given airplane velocity measured each second [MPH] over a period of one hour – Estimate the mean velocity

*
Circular* example:

- Given airplane heading measured each second [
*deg*] over a period of one hour - Estimate the mean heading

For *circular* values, it is not obvious in which “direction” the signal changed between two consecutive samples. For example, if one sample was 10°, and the next sample was 300°– did the signal changed in the increasing direction, or in the decreasing direction?

Since we can’t tell, our only way is to sample fast enough to ensure that the difference between consecutive measurements will always be less than 180°, so we can always be sure about the direction in which the signal changed. If the measure-rate is not high enough, a 330° change in the measured property would be incorrectly interpreted as a 30° change in the other direction, and will insert an error into the calculation.

Basically, in order to estimate the average, we should reconstruct the signal (interpolation), calculate the area bounded between the signal and the t-axis (integration), and divide it by .

### Linear Interpolation

The simplest reconstruction technique is linear interpolation (connecting every two consecutive samples by a straight line). We will discuss the *circular* values equivalent to *linear* values linear interpolation.

For *linear* values, the average value of a signal, based on linear interpolation, is the weighted average of

where, as defined previously,

For *circular *values, it is the weighted *circular* average of

For the above example, assume that the airplane heading measurements vector is {200, 240, 320, 40, 100, 20, 300}.

Here is how to calculate the average:

- Connect the samples such that the difference between consecutive samples will always be less than 180° (see the red line intervals in the figure above).
- Calculate the average
*circular* value, and the weight of each interval *k* in the range [2,*n*]:

In our examples the intervals average values are {220, 280, 0, 70, 60, 340}
- Calculate the weighted
*circular* average of all intervals average values

(non-weighted *circular* average may be used when the samples are equally spaced).

In our examples, assuming equally spaced sampling, it is 332.857.

Note that this method is more accurate than the Mitsuta method, described in the U.S. Meteorological Monitoring Guidance for Regulatory Modeling Applications (

http://www.epa.gov/scram001/guidance/met/mmgrma.pdf), which gives equal weight to all measurements (For interpolation – the 1

^{st} and last measurements should have only half-weight) The Mitsuta method is also limited to equally-spaced samples.

The following code calculates the average estimation based on samples:

template<typename T>
class CAvrgSampledCircSignal
{
size_t m_nSamples ;
CircVal<T> m_fPrevVal ;
double m_fPrevTime;
vector<pair<CircVal<T>, double>> m_Intervals;
public:
CAvrgSampledCircSignal()
{
m_nSamples= 0;
}
void AddMeasurement(CircVal<T> fVal, double fTime)
{
if (m_nSamples)
{
assert(fTime > m_fPrevTime);
double fIntervalAvrg = CircVal<T>::Wrap((double)m_fPrevVal +
CircVal<T>::Sdist(m_fPrevVal, fVal)/2.) ;
double fIntervalWeight= fTime-m_fPrevTime ;
m_Intervals.push_back(make_pair(fIntervalAvrg, fIntervalWeight));
}
m_fPrevVal = fVal ;
m_fPrevTime= fTime;
++m_nSamples;
}
bool GetAvrg(CircVal<T>& fAvrg)
{
switch (m_nSamples)
{
case 0:
fAvrg= CircVal<T>::GetZ();
return false;
case 1:
fAvrg= m_fPrevVal;
return true;
default:
fAvrg= *WeightedCircAverage(m_Intervals).begin();
return true;
}
}
};

Here is a sample use:

CAvrgSampledCircSignal<UnsignedDegRange> A1;
A1.AddMeasurement(CircVal<UnsignedDegRange>(200.), 1);
A1.AddMeasurement(CircVal<UnsignedDegRange>(300.), 2);
A1.AddMeasurement(CircVal<UnsignedDegRange>( 20.), 6);
CircVal<UnsignedDegRange> ad1;
A1.GetAvrg(ad1);