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.));
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,c2).
Subtraction: Starting from c1, walking with length -pdist(z,c2).
- 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
c1,
c2, for each comparison operator
, we can use:
Though all definitions satisfy all the requirements, they mean different things. The 1st compares the pdist from l, the 2nd compares the pdisht from z, and the 3rd 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 x1..xN – calculate their mean
-or-
Given a random sample x1..xN 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 x1..xN of a continuous-time [circular] signal, sampled at times t1..tN respectively, where the measurement technique/instrument is accurate - Estimate the mean value of the signal at period [t1, tN]
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 2nd 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 2nd 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: .
wi Is a positive real weight associated with value ai, and dist(v1,v2) is a function that measures the minimal distance between v1 and v2.
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);
