こんにちは，株式会社Nospareの小林です．最近個人的に方向データに関する研究に関心があるので，ここで方向統計（directional statistics），特に円周（circular）データに関する入門的な内容を簡潔にまとめておきたいと思います．
方向データとは
方向データはある基準方向からの角度で表されます．データの次元によって円周データと球面（spherical）データに分けられます．円周データは単位円上に値が観測され，それはあるひとつの基準方向に対する角度となっています．球面データでは単位球上のに値が観測され，直交する2つの基準方向に対する角度となっています．ある点から円周上をたどっていくとまた元の点にまた戻るといったような性質をもっているため，方向データに対して通常の（線形の）分析手法を適用するのは適切ではありません．また角度はあるが方向がないデータを軸データと呼びます．軸データは$[0,\pi)$に値をとります．
方向データの例
方向統計の標準的なテキストであるMarida and Jupp (2000)には各分野における方向データの例として以下のようなものが挙げられています．
- 地学
- 地震の震源（球面データ）
- 岩石の残留磁化の方向（球面データ）
- 気象学
- 風向（円周データ）
- 1日で雷雨が発生する時間（円周データ）
- 1年で豪雨が発生する時期（円周データ）
- 生物学
- 動物の移動方向（円周データ）
- 物理学
- 結晶の光学軸の方向（軸データ）
- 医学
- ベクトル心電図（球面データ）
円周データの可視化
ここではRパッケージcircularを使用します．データはpegionsという，鳩をホームから離れたところで離したときに飛んでいく方向が記録されたものを使用します．0度がホームの方向となっているので，観測値が180度の場合はその鳩はホームとは正反対の方向に飛んでいったことになります．
library(circular)
data(pigeons)
head(pigeons$bearing)
線形ヒストグラム
通常のデータの場合と同じように円周データに対してもヒストグラムを作成することは可能です．
hist(pigeons$bearing, breaks=30, xlim=c(0,360), main="", xlab="", xaxt='n')
axis(side=1, at=seq(0,360, by=45))
ただし注意しないといけないのは，円周データは360度（$2\pi$）で循環するということです．下のヒストグラムでは0度から50度あたりと300度から360度あたりまで値がたくさん観測されていますが，二峰データというわけではありません．円周をどこで切るか（基準をどこに設定するか）に依存してヒストグラムの見え方が変わってしまいます．特にモード付近で切ってしまうと適切な見方ができなくなってしまいます．
特にはっきりとモードがないデータに対しては，データをもう一周繰り返して0度から720度までの範囲で表示すると比較的見やすくなります．下のヒストグラムのように，多くの鳩が360度（0度）方向，つまりホームに向かって飛んでいっていることがわかります．
hist(c(pigeons$bearing,pigeons$bearing+360), breaks=30, xlim=c(0,720), main="", xlab="", xaxt='n')
axis(side=1, at=seq(0,720, by=45))
円周プロット（circular plot）
円周上に観測値を点としてプロットしたもので，データの構造をそのまま可視化することができます．同じ値を持つデータがある場合にはスタックして表示させることもできます．circularパッケージではデータをcircularオブジェクトに変換したものに対してplot関数を適用すると描画できます．
x <- circular::circular(pigeons$bearing, units='degrees')
plot(x, shrink=1, stack=TRUE, col=rgb(0,0,0,alpha=0.5))
ローズダイアグラム
ローズダイアグラムはヒストグラムの円周バージョンで，扇形を円状に並べたもので，各扇形の半径は対応する階級の頻度の平方根に比例したものになっています．以下は各階級の幅が15度のローズダイアグラムに円周プロットを重ねたものになります．
circular::rose.diag(x, shrink=1, prop=2, bins=24)
par(new=TRUE)
plot(x, shrink=1, stack=TRUE, col=rgb(0,0,0,alpha=0.5))
記述統計
平面上の方向は単位ベクトル，または単位円上の点として見ることができます．角度$\theta$に対する単位円上の点$x$は以下の関係で表せます（下図）：
x = (\cos\theta, \sin\theta)'
平均方向（mean direction）
単位ベクトルが$x_1,\dots,x_n$で与えられ，またそれらに対応する角度が$\theta_1,\dots,\theta_n$で与えられる．$\theta_1,\dots,\theta_n$の平均方向$\bar{\theta}$は$x_1+\cdots+x_n$の方向（重心の角度）として定義されます．$x_i$の座標は$(\cos\theta_i, \sin\theta_i)$なので，
\bar{C}=\frac{1}{n}\sum_{i=1}^n\cos\theta_i,\quad \bar{S}=\frac{1}{n}\sum_{i=1}^n\sin\theta_i
とすると，重心の座標は$(\bar{C},\bar{S})$で与えられます．よって$\bar{\theta}$は以下の方程式の解になります：
\bar{C}=\bar{R}\cos\bar{\theta}, \quad \bar{S}=\bar{R}\sin\bar{\theta},
ここで，$\bar{R}=\sqrt{\bar{C}^2+\bar{S}^2}$は平均合成ベクトル長（mean resultant length）と呼ばれます．$\bar{R}=0$の場合は平均方向が定義されず，$\bar{R}>0$の場合は
\bar{\theta}=\left\{
\begin{array}{ll}
\tan^{-1}(\bar{S}/\bar{C})&\quad \bar{C}>0,\quad \bar{S}\geq0\\
\pi/2&\quad \bar{C}=0,\quad \bar{S}>0\\
\tan^{-1}(\bar{S}/\bar{C})+\pi&\quad \bar{C}<0\\
3\pi/2&\quad \bar{C}=0,\quad \bar{S}<0\\
\tan^{-1}(\bar{S}/\bar{C})+2\pi&\quad \bar{C}>0,\quad \bar{S}<0\\
\end{array}
\right.
となります．円周データについては$\bar{\theta}$は$(\theta_1+\cdots+\theta_n)/n$とならないことに注意しないといけません．mean.circular関数を使うとpigeonsデータの平均方向は約10度となります．
また上の関係から，
\frac{1}{n}\sum_{i=1}^n\cos(\theta_j-\bar{\theta})=\bar{R},\quad \frac{1}{n}\sum_{i=1}^n\sin(\theta_j-\bar{\theta})=0
と変形できます．これは通常の算術平均における$\sum_{i=1}^n(y_i-\bar{y})=0$と同様になっています．
データを回転させる（基準方向を変える）と，平均方向も同じだけ回転します．データを$\theta_i'=\theta_i-\alpha$というように変換してやると，$\alpha$が基準方向となります．この新しい座標における平均方向$\bar{\theta}'$は$\bar{\theta}-\alpha$と等しくなります．また平均合成ベクトル長は変わりません．
円周分散（circular variance）
平均合成ベクトル長$\bar{R}$は単位ベクトルの重心ベクトル$\bar{x}$の長さなので，$0\leq \bar{R}\leq 1$を満たします．もし方向$\theta_1,\dots,\theta_n$のばらつきが小さい場合には$\bar{R}$は$1$に近くなり，ばらつきが大きい場合には$\bar{R}$は$0$に近くなります．この考え方をもとに，円周分散は
V=1-\bar{R},\quad 0\leq V\leq 1
で定義されます（$2(1-\bar{R})$という定義もあります）．pigeonsデータの円周分散は，var.circular関数を使うと約0.43となります．
２つの角度$\theta$と$\xi$があったとき，これらの距離として次ののようなものが考えられます：
1-\cos(\theta-\xi)
同様に角度データ$\theta_1,\dots,\theta_n$とある角度$\alpha$との距離として
D(\alpha)=\frac{1}{n}\sum_{i=1}^n(1-\cos(\theta_i-\alpha))
というものが考えられます．ここで$D(\bar{\theta})=V$となります．また
D(\alpha)=V+2\bar{R}\left\{\sin\left(\frac{\bar{\theta}-\alpha}{2}\right)\right\}^2
というような分解が可能になります．また$D(\alpha)$は$\alpha=\bar{\theta}$で最小となります．これも通常のデータにおける
\frac{1}{n}\sum_{i=1}^n(x_i-u)^2=\frac{1}{n}\sum_{i=1}^n(x_i-\bar{x})^2 + (\bar{x}-u)^2
という関係と同様なものになっており，これは$u=\bar{x}$で最小となります．
円周標準偏差（circular standard deviation）
円周標準偏差は以下のように定義されます．
v = \sqrt{-2\log(1-V})=\sqrt{-2\log\bar{R}},\quad v\in[0,\infty]
sd.circular関数を使うとpigeonsデータの円周標準偏差は約1.06となります．
おわりに
方向データに遭遇する場面は意外と少なくないかもしれません．方向データを扱う際にはその構造に適した可視化方法および記述統計を使用する必要があります．これは確率分布をモデルとして使用する際にも同様で，普段慣れ親しんでいる標準的な確率分布ではなく，方向データ用の確率分布を使用する必要があります．次の記事では円周データに対する確率分布に関して取り扱いたいと思います．
【参考文献】
- Mardia, K.V. and Jupp, P.E. (2000). Directional Statistics. John Wiley & Sons Ltd.
一緒にお仕事しませんか！
株式会社Nospareではベイズ統計学に限らず統計学の様々な分野を専門とする研究者が所属しており，新たな知見を日々追求しています．統計アドバイザリーやビジネスデータの分析につきましては弊社までお問い合わせください．インターンや正社員も随時募集しています！
