角度の平均
問題設定
角度を出力するセンサーがあるとします。
ただし、そのセンサーは完全に正確な値を出す理想的なセンサーではなくて、一定の誤差を含んだようなデータを出力する、現実世界のセンサーであるとします。
たとえば、本当の角度は60°なんだけど、そのセンサーで測定すると、データが58°, 59°, 63°, 61°, 57°, 64°, …と各測定で少しずつずれたデータが取得できる、という感じです。
こういった現実世界の角度センサーの出力値の平均値を求めてみましょう。
角度データは0°から360°の範囲で出力されるとします。
簡単な問題のようですが、意外な落とし穴があるものです。
うまくいかない方法 : 足してデータ数で割る
平均というのは「足してデータ数で割る」。定義に従った実装はこんな感じでしょうか。
def mean_angle(angles: list[float]) -> float:
return sum(angles) / len(angles)
60°付近でブレのあるデータ
angles = [58, 59, 63, 61, 57, 64]
に対して、この方法で平均を計算すると次のようになります。
m = mean_angle(angles)
print(m) # 60.333333333333336
うまく計算できているようです。
問題になるようなデータ
しかし、この素朴な方法には問題があるのです。
それは次のようなデータです。
angles = [358, 359, 3, 1, 357, 4]
これは0°付近でブレのあるデータです。
このデータに対して同じように計算すると次のようになります。
m = mean_angle(angles)
print(m) # 180.333333333333336
本来は、0°付近で動いているので平均値としては0°あたりに出るのが正しいはずですが、真反対の180°に出てしまっています。
問題の原因
原因は非常に簡単で、角度データが「0°から360°の範囲で出力されている」ということにあります。
0°と359°は実際には1°しか違わないのに、この表し方にすると359という数値的に大きな差が出てしまいます。
先程の方法はそれを考慮していないので、こんな困ったことが起こるのですね。
解決法?
では、先程の角度データを「-180°から180°の範囲」に直してみてはどうでしょうか。
実装としては次のようになります。
def mean_angle_2(angles: list[float]) -> float:
# anglesを-180°から180°の範囲に調整する
angles = [angle - 360 if angle > 180 else angle for angle in angles]
return sum(angles) / len(angles)
m = mean_angle_2(angles)
print(m) # 0.3333333333333333
うまくいきました。
別の問題
勘の良い方はお気づきと思いますが、この方法だと別の問題が発生します。
つまり、同じ問題が180°付近にブレを持つデータで発生するわけです。例えば
angles = [178, 179, 183, 181, 177, 184]
を考えると、
m = mean_angle_2(angles)
print(m) # 0.3333333333333333
うまくいきません。
この素朴な方法を取る場合、角度をどんな範囲で表現しても、その範囲の端っこでブレるデータについては問題が発生してしまいます。
うまくいく方法1 : 両方やって良さそうな方を取る
0°から360°の範囲で表したパターンの平均と、-180°から180°の範囲で表したパターンの平均の両方をやってみて、良さそうな方を取るというのはどうでしょうか。
ブレが常識的な範囲に収まっているなら、これでうまくいきそうです。
ただ、「良さそう」というのは、どういうことでしょうか?
ここを明確にしないと、計算はできませんね。
実は、「足してデータ数で割る」という平均(相加平均と呼ばれます)には、次のような性質があります。
相加平均の性質
n個のデータ列$x_1, x_2, … x_n$の相加平均を $$ \bar{x} = \frac{x_1 + x_2 + … + x_n}{n} $$ とする。 また、このデータ列に対し、平均二乗誤差(MSE)を次のような関数として定める。 $$ f(x) = \frac{(x - x_1)^2 + (x - x_2)^2 + … + (x - x_n)^2}{n} $$ $f(x)$は$\bar{x}$において最小値を取る。
実際、$f(x)$はただの2次関数ですから、少し変形して平方完成してあげれば高校生の知識で$\bar{x}$ で最小となることがわかります。ちなみに、最小値$f(\bar{x})$は分散と呼ばれます。
角度の平均を考えるうえで「良さそう」を数値化して考えるには、この性質を応用すればよいです。
平均二乗誤差$f(x)$の分母の各項に注目すると、これは$x$と各データ$x_i$との差(距離)の2乗になっています。
角度データで考える場合は、この「差」を単純な引き算ではなくて、図形的に考えた「差」に直せばよいのですね。
たとえば、0°と359°の差は359ではなく、1です。
このように定めた平均二乗誤差$f(x)$がより小さい方を「良さそう」と考えれば、計算することができるようになります。
具体的なコードは次のようになります。
def dist_angle(x: float, y: float) -> float:
"""2つの角度の差を求める"""
a = abs(x - y) % 360
if a > 180:
return 360 - a
else:
return a
def mse_angle(x: float, angles: list[float]) -> float:
"""角度データの平均二乗誤差(MSE)"""
return sum([dist_angle(x, a) ** 2 for a in angles]) / len(angles)
def mean_angle_3(angles: list[float]) -> float:
"""0°から360°の範囲で取った相加平均と、
-180°から180°の範囲で取った相加平均の良さそうな方を選ぶ"""
m1 = mean_angle(angles)
m2 = mean_angle_2(angles)
s1 = mse_angle(m1, angles)
s2 = mse_angle(m2, angles)
if s1 < s2:
return m1
else:
return m2
これで、うまくいくようになりました。
うまくいく方法2 : 2次元空間で考える
うまくいくようになりましたが、実装が結構複雑です。
こんな複雑な計算をしないといけないのでしょうか？
実は別の、全然違ったアプローチもあります。
そもそも、「角度をどんな範囲で表現しても、その範囲の端っこでブレるデータが問題となる」のでした。
こういった端っこの部分を「特異点」なんて言ったりしますが、そもそも角度自体には特異点はないはずなのです。
円周を思い浮かべてみてください。どこにも特別な点はありません。ただ円周があるだけです。
特異点がどうして発生するかと言うと、この円周のどこかの点を我々が勝手に決めて、その点で円周を切って、直線上に切り開いて見ているからにほかなりません。
問題を発生させている特異点は人工的なものなのです。
なので、円周を切り開かずそのままにして、角度を2次元空間で考えることはできないでしょうか？
できます。
具体的には、角度データが与えられたとき、それぞれの角度に対応する方向(単位)ベクトルを考えて、その平均を求めます。そしてその平均のベクトルに対応する角度を、平均の角度と考えるのです。
コードとしては次のようになります。
import math
def mean_angle_4(angles: list[float]) -> float:
"""2次元空間で考える"""
c = [math.cos(math.radians(a)) for a in angles]
s = [math.sin(math.radians(a)) for a in angles]
mc = sum(c) / len(c)
ms = sum(s) / len(s)
return math.degrees(math.atan2(ms, mc))
ここで計算される「平均」は、もちろんこれまで計算してきた相加平均とは異なるものです。 実際ちょっと違った値が出ます。
angles_around_0 = [358.0, 359.0, 3.0, 1.0, 357.0, 4.0]
angles_around_60 = [58.0, 59.0, 63.0, 61.0, 57.0, 64.0]
angles_around_180 = [178.0, 179.0, 183.0, 181.0, 177.0, 184.0]
print(mean_angle_3(angles_around_0)) # 0.3333333333333333
print(mean_angle_4(angles_around_0)) # 0.3331940884210799
print(mean_angle_3(angles_around_180)) # 180.33333333333334
print(mean_angle_4(angles_around_180)) # -179.66680591157893
print(mean_angle_3(angles_around_60)) # 60.333333333333336
print(mean_angle_4(angles_around_60)) # 60.33319408842108
しかし、実用的にこの差が問題になることはほぼないでしょう。
結局どうする？
「うまくいく方法2」のほうが応用はききます。
たとえば、3次元的な方向を表すのには2つの角度を使用しますが、こういった角度の組のデータを出すようなセンサーについて、同じ問題 = 平均値を計算する問題を考えるときには、「うまくいく方法2」で考えるのが良いでしょう (ぜひ考えてみてください)。
あと実装も楽です。numpyを使用するとこんな感じになります。
import numpy as np
def mean_angle_5(angles: list[float]) -> float:
"""numpyを使う"""
a = np.asarray(angles)
return np.degrees(np.angle(np.sum(np.exp(np.radians(a) * 1j))))
複素数を使っていますが、同じことです。
ただ、三角関数の計算が入っているので、リソースの厳しい組み込み環境とかパフォーマンスをカリカリにチューニングしたような場面では、少し厳しいかもしれません。
そういった場合は、「うまくいく方法1」を使うとよいでしょう。 大量のデータがあるのであれば、例えば最初のnサンプルを取って、「うまくいく方法1」でおおよその平均を計算したあとに、それと反対方向で円周を「切る」。そして全サンプルで単純に相加平均を使うなどするのが最も速くなりそうです。
MSEの代わりにMAEを使うという選択肢もあります。MAE(平均絶対誤差)は次のようなものです。
$$ g(x) = \frac{|x - x_1| + |x - x_2| + … + |x - x_n|}{n} $$
まとめとしては、まずは「うまくいく方法2」でサクッと計算する。高速化したくなってきたら、「うまくいく方法1」ベースで考える、というのが良いでしょう。