R Advent Calendar 2016の3日目です。
アドベントカレンダー初参加、Qiita投稿も初めてで、勝手が分からずすみません。おかしな点があればご指摘願います……。
Rというより統計分析の小ネタです。方向統計っぽいタイトルですが違います。方向統計は方向を目的変数にする場合で、この記事で扱うのは説明変数にする場合です。
角度変数の分類
まず角度変数について。
角度は大まかに分けて二つ、回転量と方向を表します。
回転量とは、45度の回転とか720度の回転のようにどれだけ回ったかを表します。
方向とは、基準方向を0度として、90度の方向とか240度の方向のように相対的な向きを表します。
この二つを判別するには、
- 720度が360度の2倍と言えるか（言えるなら回転量）
- 0度、360度、720度を同一視するか（するなら方向）
のように検討すればいいですね。
方向を扱う分析は用途が限られるかもしれませんが、たとえば生態学や環境学で風向の関係するケースが考えられます。あとは実験心理学とかでもあるかも……？
線形回帰の適用は？
さて、線形回帰を適用するのが妥当かどうか考えてみます。
回転量を表す角度は、周期性がなく単なる連続量です。なので、目的変数としても説明変数としても線形回帰にそのまま組み込んで問題ないでしょう。
問題は方向を表す角度です。こっちは周期性、つまりθとθ+360°×nが同等になるのを考慮する必要があります。
目的変数の場合、周期性をもつ分布（フォン・ミーゼス分布など）を用いたモデルが提案されています。方向統計とか方位統計と呼ばれているやつですね。
説明変数の場合は……？ θとθ+360°×nに対する目的変数の応答が同じということで、周期関数をあてはめるとよさそうです。
直線をあてはめるのは、θの範囲を一周分に切り取っても、両端で回帰予測値が違ってしまいます。
【図1 直線のあてはめは妥当？】
一旦ここまでを表にまとめておきます。
| 目的変数／説明変数 | 分析方法 | |
|---|---|---|
| 回転量 | 目的変数 | 連続量として扱える |
| 回転量 | 説明変数 | 連続量として扱える |
| 方向 | 目的変数 | 方向統計モデル |
| 方向 | 説明変数 | 周期関数のあてはめ （記事の続きへ！） |
正弦曲線のあてはめ
周期関数といっても無数にあるので、一番簡単な正弦曲線の場合を扱います。
もし目的変数の応答が0度を中心に対称なら、偶関数のコサイン曲線のあてはめが考えられます。
【図2 0度を中心に対称ならcos θが候補】
基準方向の0度に特別な意味がないケースでは、横方向の位置も可変がよさそうですね。
北を0度にしたり、測定者が任意に0度を決めたときは、0度が中心の応答になりそうにありません。
【図3 基準方向が任意のケースは一般の正弦曲線が候補】
高校で習った三角関数の合成公式
a\sin\theta + b\cos\theta = A\sin(\theta + \alpha) \\
ただし A = \sqrt{a^2 + b^2},\ \tan\alpha = \frac{b}{a}
を思い出してください。この公式により、一般の正弦曲線は、
y = \beta_0 + A\sin(\theta + \alpha) = \beta_0 + a\sin\theta + b\cos\theta
と表すことができます。
つまりどういうことかというと……？
一般の正弦曲線のあてはめは、sin θとcos θを線形回帰に投入するだけでOKです。
また、cos θだけにすると0度を中心に対称なあてはめになります。
線形回帰に組み込む
そういえば、この記事はRのアドベントカレンダー用でした！
Rを使って実践してみましょうか。
まず模擬データから。データサイズは50、正弦曲線に等分散の正規誤差を乗せたもの。
なお、θは度数法（度）ではなく弧度法（ラジアン）で表します。
intercept = 1
A = 1
alpha = -0.5
theta = 0:49 * (2 * pi / 50)
sigma = 0.3
y = intercept + A * sin(theta + alpha) + sigma * rnorm(50)
次に、作ったデータをプロットします。背後にある正弦曲線も点線で描きます。
plot(theta, y)
par(lty=2)
curve(intercept+A*sin(x+alpha), add=T, type="l")
【図4 わざとらしい模擬データ……^^】
このデータにモデルをフィッティングします。たった1行。ついでにsummaryも。
model = lm(y~sin(theta)+cos(theta))
summary(model)
パラメータの推定結果が表示されます。
(Intercept)の行が切片の、sin(theta)とcos(theta)の行が回帰係数の推定値にその他諸々です。
最後の仕上げです。あてはめた正弦曲線を実線で描きます。
par(lty=1)
curve(coef(model)[1]+coef(model)[2]*sin(x)+coef(model)[3]*cos(x), add=T, type="l")
【図5 あてはめた正弦曲線】
どうでしょうか？ この例では真の分布とモデルが同一なのであてはまりがいいのは当然ですが……。
まとめ
- 方向を表す角度を説明変数にするときは正弦曲線のあてはめがお手軽
- sin θとcos θを線形回帰にそれぞれ投入するだけ
- ただし、「周期的だから正弦曲線」は「増加・減少だから直線」と同じくらいに乱暴
こんなところです。最後までお読みいただきありがとうございました。
宣伝です！ 機械学習に必要な高校数学やり直しアドベントカレンダー Advent Calendar 2016にも参加しています。
9日目に離散型分布（二項分布、多項分布、超幾何分布など）について、高校数学の確率の問題と対応させてざっと紹介します。ご興味のある方はぜひ！
参考文献
- 清水邦夫, 王敏真. 環境科学における方向統計学の利用. 統計数理. 2013, 61(2). http://www.ism.ac.jp/editsec/toukei/pdf/61-2-289.pdf, (参照2016-12-02).