角度の平均や分散を複素数を用いて求める

## はじめに

### 角度の平均を単純に計算してはいけない

突然ですが、以下の角度値の平均を求める方法をご存知でしょうか？

$$
\lbrace10^\circ, 30^\circ, 350^\circ\rbrace
$$

これが普通の数値データであれば、$(10 + 30 + 350) / 3 = 130$のように平均が求まります。
しかし、角度の場合は**平均10°**となります。図示してみると、確かに10°を中心に分布しています。
なぜ単純計算では求まらないのでしょうか？

<img alt="circular-sample.png" src="https://qiita-image-store.s3.amazonaws.com/0/150557/fe3d1e9f-945c-a928-240a-fb69828df69b.png" width=300px>

これには**角度の周期性**が関係しています。350°は数値上大きく見えますが、実際には0°に近い角度です。本来連続して続いているものを1周期で切断しているため、**角度に対して素朴な統計処理をすることはできない**のです。

角度データを扱う統計の分野を「**角度統計**（Circular Statistics, Directional Statistics）」といいます。日本ではあまりメジャーでないようで、書籍も英語のものがほとんどです。とはいえ、以下の記事やPDF資料を参照することで概要はつかめることでしょう。

- [角度の平均　角度の分散について - 中田　亨 （産業技術総合研究所）](https://staff.aist.go.jp/toru-nakata/averageangle.html)
- [角度データの統計処理基礎 - 2012/1/7 第４回 定量生物学䛾会　チュートリアル](http://q-bio.jp/images/5/53/%E8%A7%92%E5%BA%A6%E7%B5%B1%E8%A8%88%E9%85%8D%E5%B8%83_qbio4th.pdf)
- [NCSS Statistical Software - Circular Data Analysis](https://ncss-wpengine.netdna-ssl.com/wp-content/themes/ncss/pdf/Procedures/NCSS/Circular_Data_Analysis.pdf)

### 角度の統計計算に複素数を使う
「角度データの統計処理基礎」のp.5にも出ているとおり、角度の平均は**それぞれの角度を取る単位ベクトルを足し合わせたベクトルの角度**です。
角度データの集合を$\lbrace\theta_j\mid j=0,1,2,...,N\rbrace$として、$\cos\theta_j, \sin\theta_j$の総和をそれぞれとれば、足し合わせたベクトルが求まります。

$$
\left( \sum_j\cos\theta_j, \sum_j\sin\theta_j \right)
$$

複素数平面の考えを用いると、この計算はもっと楽に書けます。
オイラーの公式$e^{i\theta}=\cos\theta + i\sin\theta$より、合計したベクトルは$\sum_j e^{ i\theta_j}$で表せます。この偏角が求めたい平均値なので、以下のようにシンプルな式で角度の平均を求めることができます。

$$\rm{mean}(\lbrace\theta_j\mid j=0,1,2,...,N\rbrace) = \arg\left(\sum_j e^{i\theta_j}\right)$$

### 本記事の内容と実装したコード

本記事では、複素数を使う方法で、角度における基本的な統計量を算出します。
Rには既に[circular](https://www.rdocumentation.org/packages/circular/versions/0.4-93)や[CircStats](https://www.rdocumentation.org/packages/CircStats/versions/0.2-6)のような高機能な角度統計ライブラリが存在するため、本記事ではPythonを使って実装することにします。

実装したコードの全体は以下のGistからご覧いただけます。モジュールとしてimportして使うことも可能です。

- [**kn1cht/circular.py**](https://gist.github.com/kn1cht/89dc4f877a90ab3de4ddef84ad91124e)


## 準備
以下のパッケージ・モジュールを使用します。`cmath`は複素数の計算をするために必要です。
基本的な数学関数は`numpy`に実装されているものを使用します。

```python:circular.py
import cmath
import numpy as np
import scipy.stats as stats
```

## 平均
冒頭で述べた通り、偏角が$\theta_j$である複素数の総和をとり、偏角を求めます。データの単位が違うと大変なので、度であるかラジアンであるかを引数`deg`で受け取るようにしています。単位が度の場合は最初に`np.deg2rad`でラジアンに統一し、返す際に`np.rad2deg`で度に戻します。

```python:circular.py
def mean(angles, deg=True):
    '''Circular mean of angle data(default to degree)
    '''
    a = np.deg2rad(angles) if deg else np.array(angles)
    angles_complex = np.frompyfunc(cmath.exp, 1, 1)(a * 1j)
    mean = cmath.phase(angles_complex.sum()) % (2 * np.pi)
    return round(np.rad2deg(mean) if deg else mean, 7)
```

Pythonでは`j`をつけて虚数を表す点にご注意ください。また、[`cmath.phase`は結果を[-π, π]の範囲で返す](https://docs.python.org/a/3/library/cmath.html#cmath.phase)ため、2πで割った余りを計算することで0から2πに直しました。

## 分散
以下のように$R$という数値を定めます。絶対値1である複素数の総和をデータ数$N$で割っているので、$0 \leq R \leq 1$です。

$$R = \rm{abs}\left(\frac{1}{N}\sum _ { j } e ^ { i \theta _ { j } }\right)$$

分散$V=1-R$です。$V$が大きい（$=R$が小さい）ほど、データが様々な方向にばらついているといえます。

```python:circular.py
def var(angles, deg=True):
    '''Circular variance of angle data(default to degree)

    0 <= var <= 1
    '''
    a = np.deg2rad(angles) if deg else np.array(angles)
    angles_complex = np.frompyfunc(cmath.exp, 1, 1)(a * 1j)
    r =abs(angles_complex.sum()) / len(angles)
    return round(1 - r, 7)
```

## 標準偏差
標準偏差$S = \sqrt{-2 \ln R}$です。突然logが出てきましたが、定義がこのようになっています[^1]。
[Google検索などで](https://www.google.com/search?q=sqrt(-2ln(x)))プロットを見てみると分かる通り、$R$が小さいほど大きくなり、$R=0$で無限大に発散します。

[^1]: Kanti V Mardia and Peter E Jupp. _Directional statistics_, Vol. 494. John Wiley & Sons, Inc., 2000.

明確な文献は見つけられていませんが、単位はラジアンと考えて良いようです[^2]。

[^2]: [circular statistics - How do I interpret the standard deviation of a directional dataset?](https://stats.stackexchange.com/questions/185361/how-do-i-interpret-the-standard-deviation-of-a-directional-dataset#comment352753_185511)

```python:circular.py
def std(angles, deg=True):
    '''Circular standard deviation of angle data(default to degree)

    0 <= std
    '''
    a = np.deg2rad(angles) if deg else np.array(angles)
    angles_complex = np.frompyfunc(cmath.exp, 1, 1)(a * 1j)
    r = abs(angles_complex.sum()) / len(angles)
    std = np.sqrt(-2 * np.log(r))
    return round(np.rad2deg(std) if deg else std, 7)
```

## 相関係数
2つの角度データの相関係数を求めます。[NCSSのスライド](https://ncss-wpengine.netdna-ssl.com/wp-content/themes/ncss/pdf/Procedures/NCSS/Circular_Data_Analysis.pdf)のp.12に式が載っているのでそのまま実装します。

相関が有意であるかどうかの検定も可能です。実装には[Rのcircularパッケージの該当部分](https://github.com/cran/circular/blob/master/R/cor.circular.R#L119-L125)を参考にしました。`l20`,`l02`,`l22`と相関係数$r$から検定統計量が求まり、検定統計量が標準正規分布に従うと仮定してp値を求めることができます。

```python:circular.py
def corrcoef(x, y, deg=True, test=False):
    '''Circular correlation coefficient of two angle data(default to degree)

    Set `test=True` to perform a significance test.
    '''
    convert = np.pi / 180.0 if deg else 1
    sx = np.frompyfunc(np.sin, 1, 1)((x - mean(x, deg)) * convert)
    sy = np.frompyfunc(np.sin, 1, 1)((y - mean(y, deg)) * convert)
    r = (sx * sy).sum() / np.sqrt((sx ** 2).sum() * (sy ** 2).sum())

    if test:
        l20, l02, l22 = (sx ** 2).sum(),(sy ** 2).sum(), ((sx ** 2) * (sy ** 2)).sum()
        test_stat = r * np.sqrt(l20 * l02 / l22)
        p_value = 2 * (1 - stats.norm.cdf(abs(test_stat)))
        return tuple(round(v, 7) for v in (r, test_stat, p_value))
    return round(r, 7)
```

## 様々な統計検定
角度統計では、検定も独自のものを使用します。冒頭で挙げた資料に検定も載っているものの、さすがに実装が重いため本記事では再実装しません。

もし必要になった場合は、Rに用意されている角度統計パッケージを使いましょう。PypeRというパッケージを使うと、Pythonからデータを渡してRに処理してもらうこともできます。

- [PypeRを使ってみた - Qiita](https://qiita.com/m__k/items/df83a03bc7aa1ad0a59e)

以下は、平均角度に関するパラメトリック検定であるWatson-Williams Testを使用している例です。

```python:circular_data_test.py
import pyper

r = pyper.R()
r('library(circular)')

r.assign('data1', [35, 45, 50, 55, 60, 70, 85, 95, 105, 120])
r.assign('data2', [75, 80, 90, 100, 110, 130, 135, 140, 150, 160, 165])
r('data1 <- circular(data1, units="degrees")')
r('data2 <- circular(data2, units="degrees")')

# Watson-Williams Test
print(r('watson.williams.test(list(data1, data2))'))
```

```:result
try({watson.williams.test(list(data1, data2))})

        Watson-Williams test for homogeneity of means

data:  1 and 2
F = 14.558, df1 = 1, df2 = 19, p-value = 0.001168
```