特集「環境リスクと統計解析―データ基盤構築と解析―」［総合報告］

統計数理（2013）
第 61 巻 第 2 号 289–305
c 2013 統計数理研究所


環境科学における方向統計学の利用
清水 邦夫

1,2

・王 敏真

3

（受付 2013 年 2 月 25 日；改訂 6 月 13 日；採択 6 月 14 日）

要旨

風向は環境科学における典型的で重要な角度観測値であり，また風向と共に風速やオゾン濃度などの別の観測値を伴う場合がある．本稿では，角度を含む観測値をモデル化・解析する学問である方向統計学における最新の結果を含んで総合的に報告を行うことを目的としている．
最初に，角度データに対する基本統計量としての標本平均方向および標本平均ベクトル長の計算法およびそれらの量の持つ意味と性質が線形な観測値の標本平均および標本分散と比較して述べられる．また，国立環境研究所が提供する霞ヶ浦における風向データを用いてデータのプロットと基本統計量の計算の例示を行う．Pewsey の分布対称性検定法の説明の後に，Jonesと Pewsey による柔軟な対称角度分布および正弦関数摂動法による対称分布の非対称化が紹介され，風向データを用いて当てはめ例が示される．さらに，対称／非対称性とは異なる性質として，分布の頂上付近で平坦／急傾斜の性質を表す Papakonstantinou 分布と Batschelet 分布が紹介される．説明変数に角度を含む場合の線形回帰がシリンダー上の分布から自然に導かれることを述べると共に，角度変数間の回帰モデルおよび構造モデルに関する最近の結果も紹介する．
キーワード： 角度回帰モデル，角度構造モデル，角度分布の当てはめ，シリンダー上の分布，トーラス上の分布．
1. はじめに

大気中のオゾンや粒子状物質が風によって運ばれて特定の地域に影響を与えることがあるな
どのように，気象要素の一つの風向が環境科学において欠くことのできない情報を与えることがあるのは周知の通りである．風向は典型的な角度の観測値であるので，角度変量をどのようにモデルに取り入れるかを研究することは環境科学研究の一つとして意義があると思われる．
本稿では，より一般に角度を含む観測値を統計的に扱う学問（方向統計学）について，最新の結果を含めて報告することを目的としている．なお，方向統計学について従来の研究をまとめた本としては Batschelet（1981）， Fisher（1993）， Mardia and Jupp（2000），Jammalamadaka and（2001），Zar（2010）などが知られている．
SenGupta
角度が線形変数とは著しく異なる性質を持つことは，単純に平均を取る操作においても現れる．例えば，方位を表すのに北を 0◦ とし，時計回りに一周（360◦ ）することにより，北（0◦ ）東1
2
3慶應義塾大学 理工学部：〒223–8522 神奈川県横浜市港北区日吉 3–14–1
統計数理研究所 客員

慶應義塾大学 理工学研究科：〒223–8522 神奈川県横浜市港北区日吉 3–14–1（現 統計数理研究所：〒190–8562 東京都立川市緑町 10–3）

290

統計数理

第 61 巻

第 2 号 2013（90◦ ）南（180◦ ）西（270◦ ）北（360◦ を法として 0◦ ）と表すことにしよう．北北東（22.5◦ ）と北北
西（337.5◦ ）の「平均」は直観的には北であり，これを南（(22.5 + 337.5)/2 = 180◦ ）とするのは不自然である．しかしながら，このような不合理性は，方位を（単位）円周上のベクトルと解釈し，「平均」を「ベクトルの平均」とすることにより解決される．こうして，角度の統計学は（単位）円
周上のベクトルの統計学と言うことができる．環境科学では，1 変量の角度観測値だけでなく，
風向・風速や風向・汚染物質濃度のように角度・線形変量の形式の 2 変量の観測値や 2 地点において同時に風向を観測するときのように 2 変量の角度観測値がよく現れる．したがって，方向統計学は環境科学を研究する上で極めて有用な方法の一つと考えられる．
本稿の第 2 節では，角度観測値が得られたときの基本統計量の計算の仕方と統計量の意味・性質について述べ，また，いくつかの角度分布の紹介と共にデータへの角度分布の当てはめについて紹介する．第 3 節では角度変数を含む回帰モデルについて述べる．第 4 節では，回帰とは異なるモデル化として，角度変数間の構造モデルについて現在知られている結果について紹介する．第 5 節はまとめと課題である．
2. 基本統計量の計算と分布の当てはめ得られた角度データを θ1 , θ2 , . . . , θn とする．本節では，角度データ θj (j = 1, 2, . . . , n) からどのような統計量を算出するとデータに内在するどのような情報を抽出することができるのか，
データの分布状況をモデルで表すにはどのような分布を当てはめればよいのか，について記述する．
2.1 基本統計量
数学においては，平面上の横軸の正方向の角度を 0（ラジアン）とし反時計回りに 1 周して 2πまで，区間で表せば [0, 2π) で角度（の主値）を表現するのが普通である．しかしながら，諸科学技術においては角度を区間 [−π, π) で表したり，縦軸の正方向から時計回りに計測することもある．しかし表し方がたとえ異なっていたとしても方向としての実体は変わらないはずである
ので，それぞれの表現の仕方に即して適切な解釈が可能となるように角度データの統計理論を組み立てることが必要である．そのために，角度 θ を単位円周上の点 P(cos θ, sin θ) と同一視
（図 1（a））し，点 P を大きさ（長さ）が 1 のベクトル (cos θ, sin θ) と捉える．なお， は転置を表す記号である．角度データ θj (j = 1, 2, . . . , n) を点 (cos θj , sin θj ) もしくはベクトル (cos θj , sin θj )
と同一視し，nnC=
cos θj , S =sin θjj=1

j=1

とおくとき，ベクトル (C, S) は n 個のベクトル (cos θj , sin θj ) の合成ベクトルを表し，また，(C/n, S/n) は平均合成ベクトルを表す．図 1（b）は n = 2 のとき，2 つの点もしくはベクトル
P1 と P2 から平均合成ベクトル P を得る様子を示している．合成ベクトル (C, S) もしくは平均合成ベクトル (C/n, S/n) の方向 θ は標本平均方向（sample mean direction）と呼ばれ，合成√ベクトル (C, S) の大きさ R = C 2 + S 2 の値が 0 でないとき
cos θ =C,R

sin θ =

S,
R

すなわち tan θ = S/C を満たす θ によって定義される．なお，θ は次のように陽に表すことがで
きる：

環境科学における方向統計学の利用⎧⎪
tan−1 (S/C),
⎪
⎪
⎪
⎪
⎪
⎨ π/2,
θ = tan−1 (S/C) + π,
⎪
⎪
⎪
3π/2,
⎪
⎪
⎪
⎩ tan−1 (S/C) + 2π,

291

C > 0, S ≥ 0,
C = 0, S > 0,
C < 0,
C = 0, S < 0,
C > 0, S < 0.

ここで，tan−1 (•) ∈ (−π/2, π/2) である．R の値が 0 のとき平均方向は定義されない．平均合√成ベクトル (C/n, S/n) の大きさ R = C 2 + S 2 /n は標本平均合成ベクトル長 (sample mean
resultant length) と呼ばれ，0 ≤ R ≤ 1 を満たす．R の値が 1 となるのは，すべての θj が同じ角度を表すときでそのときに限る．このことから ν = 1 − R を分布の分散（拡がりを表す量）として採用するのは自然である．実軸上のデータのときと異なり，角度データの分散 ν は 0 ≤ ν ≤ 1と上に有界である．2 つの角度はどんなに離れても原点に関して対称な点までしか離れようがないことから，角度データの分散が有界な量として定義されるのは妥当と言える．なお，R は
R = n−1 n
j=1 cos(θj − θ) と表すことができる．
平均方向は実軸上のデータの場合の算術平均に類似の性質を持つ．実際，実軸上のデータ
nx1 , x2 , . . . , xn の算術平均 x = n
j=1 xj /n に対しては偏差 xj − x の和をとるとj=1 (xj − x) = 0
nが成立し，一方，平均方向 θ に対しては j=1 sin(θj − θ) = 0 が成立する．また，角度分散は
実軸上のデータの分散 S 2 = n−1 nj=1 (xj − x)2 と類似の性質を持つことを説明しよう．実数 cnに対して U (c) = n−1 j=1 (xj − c)2 を最小にするには U (c) = n−1 nj=1 (xj − x)2 + (x − c)2 であるから c = x とすればよいことが分る．そのときの値 U (x) は分散 S 2 となる．一方，単位円
周上の 2 点の角度を θ, α とすると，2 点間の弦の長さの 2 乗が 2{1 − cos(θ − α)} と表される

ことから，V (α) = 1 − n−1 nj=1 cos(θj − α) とおくと，V (α) = ν + 2R sin2 {(θ − α)/2} を得るので，V (α) は α = θ のときに最小となる．そのときの最小値 V (θ) は角度分散 ν となる．なお，
tj = (cos θj , sin θj ) , j = 1, . . . , n, とし，t = (C, S) とすると，恒等式
2
n 

1 
tj − t 
2(1 − R) =n j=1  t 
が成り立つ．ただし， •  は a = (a1 , a2 ) のときノルム  a = a21 + a22 を表す．この恒等式の右辺は，大きさ 1 のベクトル tj (j = 1, 2, . . . , n) と大きさ 1 のベクトル t/  t （標本平均方向を表す合成ベクトル t の基準化ベクトル）の差の 2 乗ノルムの平均であることを示している．
2.2 プロット
環境科学において典型的な角度データである風向データを用いて例示する．データとして国立環境研究所が提供する霞ヶ浦臨湖実験施設気象データサイト http://www.cger.nies.go.jp/db/kasumi/
index12.html に掲載されている 2012 年 4 月 1 日 ∼ 5 月 31 日の毎日朝 6 時と正午 12 時に本館

図 1. （a）角度と単位円周上の点の対応，（b）n = 2 のときの平均合成ベクトル P．

統計数理

292

図 2.

第 61 巻

第 2 号 2013

霞ヶ浦朝 6 時（左）と正午 12 時（右）の風向データ（2012 年 4 月 1 日 ∼ 5 月 31 日）の円周プロット．

（地上）において観測された風向風速データ（時間平均データ）中の最多風向を用いる．風向は 16方位を用いて表されている．図 2 は，北の風を 0◦ とし時計回りに方位を配置（北 東 南 西 北）したときの円周プロットを表している．標本平均方向 θ の値は，朝 6 時の風向データに対して
0.541（ラジアン）= 31.03◦ で，正午 12 時に対して 1.702（ラジアン）= 97.51◦ である．また，標本
平均合成ベクトル長 R の値は，朝 6 時に対して 0.237，正午 12 時に対して 0.506 となる．したがって，標本角度分散はそれぞれ 0.763 と 0.494 であり，朝 6 時の方が分布の拡がりが大きい．
言い換えれば，正午 12 時の分布の方が一方向に集中する傾向がある．
2.3 対称分布の当てはめ
2.3.1 対称性の検定√
Θ を角度確率変数とするとき，Θ の p 次三角モーメントは E(eipΘ ) (i = −1) で定義され，1次の三角モーメント E(eiΘ ) = ρ eiμ (ρ > 0) における μ は Θ の平均方向（mean direction），ρ は
Θ の平均合成ベクトル長（mean resultant length）と呼ばれる．ρ = 0 のとき平均方向は定義されない．μ の周りの 2 次正弦モーメント（second sine moment）は β 2 = E{sin 2(Θ − μ)} で定義され，分布の歪みを量る歪度（skewness）は β 2 と ρ を使って β 2 /(1 − ρ)3/2 で定義される．Pewsey（2002）は角度データ θ1 , θ2 , . . . , θn から「対称性」の仮説 H0 : β 2 = 0 をノンパラメトリック検定する方法を以下のように与えている．

標本平均方向を θ とし，θ 周りの標本 2 次正弦モーメントを b2 = n−1 nj=1 sin 2(θj − θ) とする．また，b2 の分散の推定値を代入法によりvar(bˆ 2) =

1 1 − a4
2a2
− 2a2 +n2
Rとする．ここで，1
ap =cos p(θj − θ)n j=1n(p = 2, 3, 4),a3 +a2 (1 − a2 )R2  n
2
 n


1

R=
cos θj+sin θjnj=1
j=1を表す．検定統計量を
b2
Z= 
var(b
ˆ 2)

とすると，Z は帰無仮説の下で近似的に標準正規分布に従うので，データから計算される Z の
値 z に対し片側 P 値 P (U > |z|) を計算する．U は標準正規分布に従う確率変数を表す．対称性の仮説の採択／棄却を P 値によって判断する．仮説 H0 の有意水準 α の検定では，uα を標準正規分布の上側 100α％点として |z| > uα のとき H0 を棄却する．

環境科学における方向統計学の利用

293

2.3.2 Jones と Pewsey の分布
データへの分布の当てはめが可能であるように「データのサイズは十分に大きい」ものとしておく．柔軟性を持つ対称分布の一つに Jones and Pewsey（2005）の分布があり，確率密度関数は次式で与えられる：(2.1)

fψ (θ) ={cosh(κψ) + sinh(κψ) cos(θ − μ)}1/ψ,
2πP1/ψ (cosh(κψ))

0 ≤ θ < 2π.

ここで，κ ≥ 0，0 ≤ μ < 2π ，−∞ < ψ < ∞（極限の意味で ψ = 0 も含む）であり，P1/ψ は Legendreの陪関数
 π(z + z 2 − 1 cos x)1/ψ dx = πP1/ψ (z), ψ > 0, π01√dx = πP−1/ψ−1 (z) = πP1/ψ (z), ψ < 0
2 − 1 cos x)−1/ψ(z+z0を表す．分布族（2.1）は特別な場合として以下の分布を含むことが知られている．
：（1）一様分布（κ = 0; ψ → ±∞, κ 有限）f (θ) =

1.
2π（2）von Mises 分布 VM(μ, κ) (ψ → 0)：
f (θ) =

1eκ cos(θ−μ) ,2πI0 (κ)

I0 (κ) =

1
2π

 2π

eκ cos θ dθ.

0

I0 は 0 位の第 1 種変形 Bessel 関数を表す．
（3）巻き込み Cauchy 分布 WC(μ, ρ) (ψ = −1, ρ = tanh(κ/2))：f (θ) =

1
1 − ρ2=.
2
2π cosh(κ){1 − tanh(κ) cos(θ − μ)} 2π{1 + ρ − 2ρ cos(θ − μ)}（4）ハート型分布 C(μ, ρ) (ψ = 1, ρ = tanh(κ)/2)：f (θ) =

1
{1 + 2ρ cos(θ − μ)} .
2π（5）Cartwright（1963）の分布（ψ = 1/m, κ → ∞, m は負でない整数）：
f (θ) =

Γ(m + 1)√
{1 + cos(θ − μ)}m .
2m+1 π Γ(m + 1/2)

（6）自由度 d (> 0) の円周 t 分布 (−1 < ψ < 0)（Shimizu and Iida, 2002）：
f (θ) =

{1 + tanh(κψ) cos(θ − μ)}1/ψ,
2π 2 F1 (−1/(2ψ), −1/(2ψ) + 1/2; 1; tanh2 (κψ))

d11=− − .
ψ
2
2

2 F1 は Gauss の超幾何関数を表す．

図 3 に，μ = π, κ = 1.5 のとき，ψ = −1（巻き込み Cauchy 分布），0（von Mises 分布），1（ハート型分布）の値に対して Jones と Pewsey の分布の確率密度関数（2.1）のプロットを与える．
データ θ1 , θ2 , . . . , θn から分布族（2.1）のパラメータ μ, κ, ψ を最尤推定するには，例えば最適化の手法により，尤度n
L(μ, κ, ψ) =fψ (θj )
j=1

統計数理

294

図 3.

第 61 巻

第 2 号 2013

Jones と Pewsey の分布の確率密度関数（2.1）のプロット：μ = π, κ = 1.5；ψ = −1
（巻き込み Cauchy：破線），0（von Mises：実線），1（ハート型：点線）．

を μ, κ, ψ に関して最大化する．モデル選択のためには，分布族（2.1）をフルモデル，（1）∼（6）を（2.1）のサブモデルとして AIC（= −2 ×（最大対数尤度）+ 2 ×（自由パラメータ数））や BIC（= −2 ×（最大対数尤度）+ log n ×（自由パラメータ数））の値を最小とするモデルを採用すればよい．n は標本の大きさを表す．また，選択されたモデルの適合の良さを調べるには， χ2 適合
度統計量の P 値を計算する．
Abe et al.（2010）は Jones and Pewsey（2005）の分布の拡張形を与え，例として Jones andPewsey（2005）で取り上げられた Mardia and Jupp（2000, Table 1）の渡り鳥の飛翔方向データ(n = 714) を用いて解析を行っている．しかし，渡り鳥の飛翔方向データへの分布の適合に関しては僅かに改善が見られる程度であった．
2.3.3 Papakonstantinou 分布と Batschelet 分布
分布の頂上付近が平坦であったり，反対に急傾斜であったりすることが観察されるときには，
Papakonstantinou 分布（Papakonstantinou, 1979）や Batschelet 分布（Batschelet, 1981）の利用が考えられる．
Papakonstantinou 分布の確率密度関数は(2.2)

f (θ) =

1
[1 + κ cos{(θ − μ) + ν sin(θ − μ)}],2π{1 − κJ1 (ν)}

で与えられる．Jr は r 位の第 1 種 Bessel 関数

1 π
Jr (z) =cos(rθ − z sin θ)dθ,
π 0−π ≤ θ < π

r = 0, 1, 2, . . .

を表す．μ (−π ≤ μ < π) は平均方向，κ (0 ≤ κ ≤ 1) は分布の集中度を表すパラメータで，ν (−∞ <（一様分布でない）のとき，
ν < ∞) は分布の形と単峰・多峰を制御するパラメータである．κ = 0
分布が単峰であるための必要十分条件は |ν| ≤ 1 である．ν = 0 であれば，（2.2）がハート型分布に帰着することは明らかであろう． −1 ≤ ν < 0 ならばハート型分布より平坦であり， 0 < ν ≤ 1
ならば，より急傾斜となる． Papakonstantinou 分布の性質や適用例は Abe et al.（2009）に見ることができる．図 4 に，μ = π ，κ = 0.75，ν = −1，0（ハート型分布），1 の組合せに対して
Papakonstantinou 分布の確率密度関数のプロットを与える．
Batschelet 分布の確率密度関数は(2.3)

f (θ) = C −1 exp[κ cos{(θ − μ) + ν sin(θ − μ)}],−π ≤ θ < π環境科学における方向統計学の利用

図 4.

295

Papakonstantinou 分布の確率密度関数（2.2）のプロット：μ = π, κ = 0.75；ν = −1（点線），0（ハート型：実線），1（破線）．

である．ここで，−π ≤ μ < π, κ ≥ 0, −∞ < ν < ∞ であり，C は積分 πC=
exp{κ cos(θ + ν sin θ)}dθ−πで与えられる正規化定数を表す．ν = 0 であれば，（2.3）は von Mises 分布に帰着する．Batschelet
分布の単峰性の必要十分条件は Papakonstantinou 分布の場合と同じ |ν| ≤ 1 であり，−1 ≤ ν < 0
ならば von Mises 分布より平坦であり，0 < ν ≤ 1 ならば，より急傾斜となる．Batschelet 分布のパラメータ推定と適用例は Pewsey et al.（2011）に与えられている．Batschelet 分布の確率密度関数の形状は Papakonstantinou 分布と似たものになる．
2.4 非対称分布の当てはめ
第 2.3.1 節で対称性の検定を述べ，第 2.3.2 節で柔軟性を持つ対称分布の当てはめについて紹介し，第 2.3.3 節では別の形の対称分布の利用可能性について簡単に触れた．対称性の検定によって対称性が棄却された場合は非対称分布を当てはめることになるであろうし，また対称性に拘ることなくデータに対して良い適合を与える分布に関心が持たれる場合は，初めから非対称分布の適合を考えるかも知れない．本節では，第 2.3.2 節および第 2.3.3 節の分布の非対称化について述べる．得られる分布はそれぞれの対称分布を特別な場合として含んでいる．
2.4.1 Jones と Pewsey の分布の正弦関数摂動
Abe and Pewsey（2011）において Jones and Pewsey（2005）の分布の非対称化が与えられている．非対称分布構成のアイディアは単純明快であり，（2.1）の fψ (θ) を正弦 (sin) 関数によって(2.4)

f (θ) = {1 + λ sin(θ − μ)}fψ (θ),−π ≤ θ < π, −π ≤ μ < π, −1 ≤ λ ≤ 1と摂動するものである．確率密度関数（2.4）は一般に非対称（λ = 0 ならば対称）であり，κ > 0 お√
よび −1 ≤ λ ≤ 1 に対して 4 − 3 3 ≤ ψ ≤ −0.5 ならば単峰である． Jones と Pewsey の分布が特別な場合として von Mises，巻き込み Cauchy，ハート型の各分布を含んでいることから，Abeと Pewsey の分布は正弦関数摂動 von Mises，巻き込み Cauchy，ハート型の各分布を特別な場合として持つことは明らかである．図 5 に，μ = π ，κ = 1，λ = 0（von Mises 分布），0.5，1 の組合せに対して von Mises 分布の正弦関数摂動確率密度関数のプロットを与える．
（2.4）の分布をフルモデルとし，正弦関数摂動 von Mises，巻き込み Cauchy，ハート型の各分
布および Jones と Pewsey の対称分布をサブモデルとして AIC，BIC 最小化法により最適モデルを選択する方法が考えられる．Abe et al.（2012b）では，北西フィンランド Pallas-Ylläs National
Park と北西ロシア Dvina-Pinega 域における Picea abies 倒木方向データに対して上記の対称お

統計数理

296

第 61 巻

第 2 号 2013

図 5.

von Mises 分布の正弦関数摂動確率密度関数のプロット：μ = π, κ = 1；λ = 0（von Mises：実線），0.5（破線），1（点線）．

表 1.

霞ヶ浦風向データへの正弦関数摂動分布の当てはめ（上段：朝 6 時，下段：正午 12 時）．

よび非対称 von Mises と巻き込み Cauchy の分布を当てはめて対称 von Mises 分布を最適な分布として選択している．また，北西ロシアの Kazkim 川流域のデータでは正弦関数摂動巻き込
み Cauchy 分布を選択している．
2.4.2 風向データへの当てはめ例
第 2.2 節の霞ヶ浦風向データに対する Pewsey（2002）の対称性の検定（第 2.3.1 節）の P 値は，
朝 6 時の風向データに対して 0.013 で，正午 12 時に対して 0.008 であった．この結果から，データに非対称分布を当てはめることとし，第 2.4.1 節の方法を適用してみよう．表 1 に，正弦関数
摂動 Jones–Pewsey 分布をフルモデルとし，正弦関数摂動 von Mises，巻き込み Cauchy，ハート型分布をサブモデルとして霞ヶ浦朝 6 時と正午 12 時の風向データへの当てはめを行ったときの最尤推定値，最大対数尤度の値，AIC と BIC の値を与える．6 時と 12 時の双方とも AIC
最小の意味ではフルモデルが選択される．また，図 6 は，円周プロット上で分布の当てはめを示している．
2.4.3 非対称 Papakonstantinou 分布と非対称 Batschelet 分布
第 2.3.3 節の Papakonstantinou 分布と Batschelet 分布の確率密度関数における cos-sin の組
合せを sin-sin もしくは本質的に同等であるが cos-cos とするだけで非対称性を表すようにできる．具体的な確率密度関数は次の通りである：
(2.5)
(2.6)

1
[1 + κ sin{(θ − μ) + ν sin(θ − μ)}],2πf (θ) = C −1 exp[κ sin{(θ − μ) + ν sin(θ − μ)}].

f (θ) =

非対称 Papakonstantinou 分布では正規化定数に特殊関数を含まないが，非対称 Batschelet 分

環境科学における方向統計学の利用

図 6.

朝6時
（上段）
と正午 12 時（下段）の風向データに対して，左から正弦関数摂動 Jones-Pewsey
分布，von Mises 分布，巻き込み Cauchy 分布，ハート型分布の当てはめ．

図 7.

非対称 Papakonstantinou 分布の確率密度関数のプロット：μ = π/2, κ = 0.75；ν = 0（ハート型：実線），0.5（破線），1（点線）．

布の正規化定数では一般に積分

297

 π

C=

exp{κ sin(θ + ν sin θ)}dθ−πで表現されることになる．ν = 0 であれば，（2.5）と（2.6）はそれぞれ明らかにハート型と von（ハート型分布），0.5，1 の組合せにMises 分布に帰着する．図 7 に，μ = π/2，κ = 0.75，ν = 0
対して非対称 Papakonstantinou 分布の確率密度関数のプロットを与える．非対称 Batschelet 分
布の確率密度関数の形状は非対称 Papakonstantinou 分布の場合に似たものになる．
3. 角度変数を含む回帰モデル
3.1 角度／線形説明変数・線形目的変数の回帰モデル
説明変数に三角関数を含み，目的変数が線形であるような回帰モデルは Johnson and Wehrly（1978）によって与えられている．モデルを導くための基礎となるのは次の事実である．

統計数理

298

第 61 巻

第 2 号 2013

角度ベクトルを Θ = (Θ1 , Θ2 , . . . , Θp ) とし，⎞⎛cos Θ1 · · · cos nΘ1 sin Θ1 · · · sin nΘ1
⎟
⎜ .
..
..
..
⎟
H(Θ) = ⎜
.
.
.
⎠
⎝ ..
cos Θp · · · cos nΘp sin Θp · · · sin nΘpとおき，線形ベクトルを X = (X1 , X2 , . . . , Xq ) とする．Θ と X は結合確率密度関数(3.1)

1
f (θ, x) = c−1 exp − x Σ−1 x + λ Σ−1 x + a(θ) Σ−1 x
2を持つとする．ここで，c は正規化定数で，θ ∈ [0, 2π)p ，x ∈ Rq ，λ = (λ1 , λ2 , . . . , λq ) ∈ Rq ，a(θ) = (a1 (θ), a2 (θ), . . . , aq (θ)) ∈ Rq ，ai (θ) =pnaijk cos{k(θj − μijk )}j=1 k=1

=

p
n



{αijk cos(kθj ) + βijk sin(kθj )},i = 1, . . . , q

j=1 k=1である．αijk と βijk は αijk = aijk cos(kμijk )，βijk = aijk sin(kμijk ) を表す．また，Σ−1 は q × q
正定値行列である．
（3.1）の結合確率密度関数 f (θ, x) は E(XX  )，E(X)，E{X ⊗ H(Θ)} が指
定された値を持つという制約条件の下にエントロピー − f (θ, x) log f (θ, x)dθdx を最大化する
確率密度関数として特徴付けられる．なお，⊗ はクロネッカー積を表す．
いま，X = (X 1 ; X 2 ) と分割 (X 1 : r × 1, X 2 : (q − r) × 1) し，対応して λ = (λ1 ; λ2 ) ，a(θ) =Σ11 Σ12(a1 (θ) ; a2 (θ) ) ，Σ =と分割すると，X 2 = x2 と Θ = θ が所与のときの X 1 の条
Σ21 Σ22
件付き分布は平均ベクトル λ1 + a1 (θ) + Σ12 Σ−1
22 {x2 − (λ2 + a2 (θ))} で分散共分散行列 Σ11 −Σ12 Σ−1
22 Σ21 の r 次元正規分布に従うことが分る．特に，各成分 Xi (i = 1, 2, . . . , r) の平均は
ν0 +

q


νi x i +

i=r+1

q
p
n




{γijk cos(kθj ) + δijk sin(kθj )}i=r+1 j=1 k=1の形となる．上の事実より，角度変数 θ1 , θ2 , . . . , θp に対して正弦・余弦変換を施して，角度変数ベクトル θ と線形変数ベクトル x2 を説明変数とする線形回帰モデルを構成できたことになる．したがって，通常の回帰モデル分析の手続きを適用することにより，モデルのパラメータを推定することができる．
特別な場合として，角度変数 Θ と線形変数 X の次の結合確率密度関数(3.2)

f (θ, x) = c−1 exp −

x2
λx
κx
+ 2 + 2 cos(θ − μ) ,
2σ 2
σ
σ

0 ≤ θ < 2π, −∞ < x < ∞を持つシリンダー上の分布 (Johnson and Wehrly, 1978) を考えてみる．ここで，c は正規化定数であり，μ, σ 2 , λ, κ はパラメータで 0 ≤ μ < 2π, σ > 0, −∞ < λ < ∞, κ > 0 である．周辺分
布はよく知られた形ではないが， X = x が所与のときの Θ の条件付き分布は von Mises 分布
VM(μ, κx/σ 2 ) となり，また Θ = θ が所与のときの X の条件付き分布は平均 λ + κ cos(θ − μ)，
分散 σ 2 の正規分布となる．したがって，説明変数 x（線形ベクトル）と θ（角度），目的変数 y（線形）に対してy = x β + κ cos(θ − μ)の形の回帰モデルの利用が自然となることが分る．ここで，β はパラメータベクトル，κ と μはパラメータである．なお，Johnson and Wehrly（1978）では，解析例として大気質指標 x1 を

環境科学における方向統計学の利用

299

温度 x2 と風向 θ で説明する回帰式 x1 = a + bx2 + c1 cos θ + c2 sin θ が与えられている．ここにおいて a, b, c1 , c2 はパラメータを表す．
3.2 角度／線形説明変数・線形目的変数の変形回帰モデル
対称 Papakonstantinou 分布もしくは対称 Batschelet 分布の確率密度関数に現れる形を参考にして，前節の（3.2）式をx2
λx
κx+ 2 + 2 cos{θ − μ + ν sin(θ − μ)} ,2σ 2σσ
0 ≤ θ < 2π, −∞ < x < ∞f (θ, x) = c−1 exp −と修正してみよう．そうすると，Θ = θ が所与のときの X の条件付き分布は平均 λ + κ cos{θ −μ + ν sin(θ − μ)}，分散 σ 2 の正規分布となることが分る．同様に非対称 Papakonstantinou 分布
もしくは非対称 Batschelet 分布の確率密度関数に現れる形を参考に cos-sin の組合せを sin-sin
もしくは cos-cos の組合せに修正することができる．こうして，SenGupta and Ugwuowo（2006）
は，説明変数 x（線形ベクトル）と θ（角度），目的変数 y（線形）に対して次の形の回帰式の使用を提案している：
y = x β 1 + a1 cos{θ − μ1 + ν1 sin(θ − μ1 )},y = x β 2 + a2 cos{θ − μ2 + ν2 cos(θ − μ2 )}.

ここで，(β1 , a1 , μ1 , ν1 ) と (β 2 , a2 , μ2 , ν2 ) は，それぞれのモデルのパラメータである．誤差が正規分布 N (0, σ 2 ) に従うと仮定するときには，パラメータ β 1 , a1 , μ1 , ν1 , σ 2 と β 2 , a2 , μ2 , ν2 , σ 2 の最尤推定値を最適化の手法で求めることができる．
3.3 線形説明変数・角度目的変数の回帰モデル
Fisher and Lee（1992）と Gill and Hangartner（2010）は， von Mises 分布の平均方向パラメー
タが共変数の一次結合として表されると仮定するモデル化について考察を与えている．角度観測値 θ1 , θ2 , . . . , θn は von Mises 分布 VM(μj , κ) からのものとし，各平均方向 μj は
μj = μ + g(β xj ),g(z) = 2 tan−1 zと表されると仮定する．集中度を表すパラメータ κ は j に依らず一定と仮定している．ここで，xj は p 次元共変数ベクトルであり，β は p 次元パラメータベクトルを表す．データ θ1 , θ2 , . . . , θnに基づいてパラメータ μ, κ, β を最尤推定するためには，尤度関数
L(μ, κ, β|θ1 , θ2 , . . . , θn ; x1 , x2 , . . . , xn ) = −n log{2πI0 (κ)} + κn

cos{θj − μ − g(β xj )}j=1

を μ, κ, β に関して最大化する．
いま，uj = sin{θj − μ − g(β xj )}, u = (u1 , u2 , . . . , un ) , X = (x1 , x2 , . . . , xn ) ,








dg(t) dg(t) dg(t) 
,
,
.
.
.
,
,
G = diagdt  
dt  
dt  
S=

n
1

n j=1

t=β x1t=β x2t=β xnsin{θj − g(β xj )},C=cos{θj − g(β xj )},

n
1

n j=1

R = (S 2 + C 2 )1/2

とおくとき，最尤推定方程式はX  Gu = 0,

R sin μ̂ = S,

R cos μ̂ = C,A(κ̂) = R

統計数理

300

第 61 巻

第 2 号 2013

 2πで与えられる．ここで，A(κ) は，Ir (κ) を r 位の第 1 種変形 Bessel 関数 Ir (κ) = (2π)−1 0
cos(rθ)eκ cos θ dθ として，A(κ) = I1 (κ)/I0 (κ) を表す．Fisher and Lee（1992）は，最尤推定方程式の数値解を求めるために反復再重み付け最小二乗法を用いている．
3.4 角度説明変数・角度目的変数の回帰モデル
説明変数と目的変数が共に角度であるときの回帰モデルは Downs and Mardia（2002）， Katoet al.
（2008），Kato and Jones（2010）で取り扱われている．回帰モデルを表現するために Möbius
変換を用いるので，複素数を持ち出すと数学的扱いが簡潔になる部分がある．以下では，通常の角度表現による解釈も与える．
複素平面上の単位円周を Ω = {z ∈ C; |z| = 1} とし，j = 1, 2, . . . , n に対し Ω 上の共変数を√
xj = eiξj ，独立同一分布する Ω 上の誤差を εj = eiΦj と複素数表示する．i は虚数単位 i = −1を表す．また，Ω 上にない複素定数 α = ρα eiμα ∈ C (ρα ≥ 0, ρα = 1, −π ≤ μα < π) と Ω 上の複
素定数 β = eiμβ ∈ Ω (−π ≤ μβ < π) を取る．このとき，Ω 上の目的変数を Yj とする回帰モデル(3.3)

Yj = βxj + αεj1 + αxj

を考えることができる．ここで，α は α の共役複素数を表す．変換 M (z) = (z + α)/(1 + αz) は，単位円周上の点を単位円周上に写す Möbius 変換として知られている．
（3.3）の表現を用いれば複素関数論の知識を使って単位円周上の回帰モデル理論を展開することが可能となるが，ここでは（3.3）の角度表現Θj ≡ arg(Yj ) = ν(μα , μβ , ρα , ξj ) + Φj ,
11 − ραtan (ξj − μα ) + (μα + μβ )(3.4)
ν(μα , μβ , ρα , ξj ) = 2 arctan
1 + ρα2を用いることにする．
（3.4）は，単位円（点 (0, 1) を除く）と実軸の 1 対 1 対応を表す立体射影（図 8）を用いて次のように解釈することができる．角度 ξj − μα を立体射影により実軸に写した点を (1 − ρα )/(1 + ρα ) 倍して尺度を変え，次に逆立体射影により実軸から単位円に写した後に角度 μα + μβ 分を回転することにより角度が ν となる．
Downs and Mardia（2002）は Φj に von Mises 分布 VM(0, κ) を仮定し，データ (θj , ξj ) から対数尤度関数
L(μα , μβ , ρα , κ) = −n log I0 (κ) + κn

cos{θj − ν(μα , μβ , ρα , ξj )} + const.

j=1を最大化し，パラメータ μα , μβ , ρα , κ を最尤推定する方法について研究した．また，Kato et al.
（2008）は Φj に巻き込み Cauchy 分布 WC(0, ρ) を仮定することを提案している．数値計算の手

図 8.

立体射影 x = tan(θ/2), −π < θ < π ．

環境科学における方向統計学の利用

301

間としては誤差分布として VM を仮定しても WC を仮定しても大差はないが，WC では WC
変数の Möbius 変換は再び WC 分布に従い，WC 分布は再生性を持ち，また，ξj が所与のときの Θj の条件付き分布は WC 分布となることが知られている．これらのように，誤差分布として WC を仮定すると Möbius 変換と相性が良いので理論展開がスムーズに進むというメリットがある．
さらに，Kato and Jones（2010）は VM 変数を Möbius 変換した分布を誤差分布として採用することを提案した．その誤差分布は一般に非対称，2 峰性であり，VM と WC 分布を特別な場合として含むという大きなメリットを持つ．少し複雑な分布ではあるが，角度変数間回帰モデ
ルによるデータ解析においてはこの分布を誤差分布として採用することが推奨される．
4. 角度変数間の構造モデル

線形変数に関する構造モデルに類似して，角度変数間の構造モデル（Wang and Shimizu, 2012）θ = ξ + δ,

φ=η + εを考えよう．ξ と η は角度確率変数であり，関係 η = M (ξ) を持つとする．また，δ と ε は観測誤差を表す独立な角度確率変数で，ξ とも独立とする．簡単のために，ξ はハート型分布 C(μ, ρ)
に従い，δ と ε はそれぞれ平均方向 0 のハート型分布 C(0, ρ1 ) と C(0, ρ2 ) に従うと仮定する．関係 η = M (ξ) は角度と角度の間の関係であるので，（3.4）の右辺で現れている Möbius 変換を用いることにする．
確率ベクトル (θ, φ) の分布は 2 変量角度分布（トーラス上の分布）を表す．この分布を MC 2 (μ, ρ,
μα , ρα , μβ , ρ1 , ρ2 ) と表記することにしよう． θ の周辺分布はハート型 C(μ, ρρ1 ) であり，φ の周辺分布はハート型 C(μφ , ρφ ) である．ここで，μφ = μα + μβ + tan−1ρ(1 − ρ2α ) sin(μ − μα )
ρα + ρ|1 − ρ2α | cos(μ − μα ),

ρφ = ρ2 ρ2 (1 − ρ2α )2 + ρ2α + 2ρρα (1 − ρ2α ) cos(μ − μα )!1/2を表す．MC2 (μ, ρ, μα , ρα , μβ , ρ1 , ρ2 ) に従う (θ, φ) の結合確率密度関数は幾分面倒な形であるが，陽に表すことができる．最初に仮定した ξ, δ, ε の分布はハート型と限定されるわけではなく vonMises や巻き込み Cauchy 分布でも構わないが，ハート型以外のときにモデルの統計的性質を導くことは現在のところ成功していない．
Wang and Shimizu（2012）の方法に倣い，本稿では例として第 2.2 節の霞ヶ浦風向データ中，
朝 6 時と正午 12 時の風向データの解析例を以下に示す．表 2 に，（f）MC2 (μ, ρ, μα , ρα , μβ , ρ1 , ρ2 )をフルモデルとし，（a）恒等変換 M (ξ) = ξ ，(b) 角度変換 M (ξ) = ξ + μβ ，2 つの独立モデル（c）
MC2 (μ, ρ, μα , ρα , μβ , 0, ρ2 ) と（d）MC2 (μ, ρ, −, −, −, ρ1 , 0)，（e）モデル MC2 (μ, ρ, μα , ρα , μβ , ρ1 , ρ2 )，ρ1 = ρ2 をサブモデルとするときの最尤推定値，最大対数尤度の値，AIC と BIC の値を与える．
AIC と BIC 最小の意味では（e）MC2 (μ, ρ, μα , ρα , μβ , ρ1 , ρ2 )，ρ1 = ρ2 ，のモデルが選択される．
図 9 は選択されたモデルの等高線と構造曲線の図と 3 次元プロットを示している． μβ の値は0.933 と推定されている．この値は標本平均方向の差 1.161 (= 1.702 − 0.541) に近い．しかしな
がら，選択されたモデルは朝 6 時と正午 12 時の間に平均方向のズレだけでなく，点 α = ρα eiμαを含む変換を通した構造の存在を示唆している．なお，選択されたモデルでは ρ̂1 = ρ̂2 = 0.5 であり，ハート型分布 C(μ, ρ) の ρ の範囲 0 ≤ ρ ≤ 1/2 の上限に達している．このことはハート型
よりも急傾斜な性質を持つ分布が採用されるべきであることを示唆しているので，ξ, δ, ε の分
布がハート型以外のときの理論構築が望まれる．

統計数理

302

第 61 巻

第 2 号 2013

表 2. （a）恒等変換 M (ξ) = ξ ，（b）角度変換 M (ξ) = ξ + μβ , 独立モデル（c）MC2 (μ, ρ,
μα , ρα , μβ , 0, ρ2 ) と（d）MC2 (μ, ρ, −, −, −, ρ1, 0)，（e）モデル MC2 (μ, ρ, μα , ρα ,μβ , ρ1 , ρ2 ), ρ1 = ρ2 ，（f）フルモデル MC2 (μ, ρ, μα , ρα , μβ , ρ1 , ρ2 ) におけるモデル
選択．

図 9.

霞ヶ浦朝 6 時と正午 12 時の風向データに対する構造モデルの推定（左：推定された等高線と構造曲線，右：3 次元プロット）．

5. おわりに

本稿では，環境科学においてしばしば現れる角度データに対して，データに含まれる情報を抽出するためにはどのような統計量を算出すればよいのか，また計算された統計量はどのような意味・意義を持つのかから始めて，角度分布とその当てはめ，シリンダー上の分布から導かれる角度回帰モデル，角度構造モデルとトーラス上の分布について述べた。これらは，角度を含むデータを扱う方向統計学の一部をなしている．とは言え，環境統計学において現れる角度を含むデータは本稿で述べたものに限られるわけでないことは明らかである．他にも興味のあるデータの型はさまざまな場面で現れる．例えば，Abe et al.（2012a）は倒木データの向きを軸データとして捉えたときの解析方法とデータ解析例を示している．関心のある読者のために，
軸分布に関する文献として Arnold and SenGupta（2006, 2011）を挙げておこう．また Abe et al.
（2012b）は，第 2.4.1 節にあるように，通常の角度分布を用いて倒木データを解析している．島谷編（2011）は統計数理研究所で行われた共同利用研究集会「生物群集の多様性と統計数理」の講演のテープを起して作成された共同研究リポートである．その中に阿部と加藤のそれぞれの講演が第 4 章に方向統計学として文章化されており，倒木データの角度統計解析と第 3.4 節のKato and Jones（2010）の角度回帰モデルが分り易く解説されているので，環境統計学・生態学
への方向統計学の利用に関心を持つ読者に一読を薦めたい．また，島谷（2012）の第 8 章に鳥の飛翔方向データへ角度分布を適用した例が紹介されている．

環境科学における方向統計学の利用

303

環境科学では典型的な（風向，風速）のモデル化にはシリンダー上の分布を欠くことができな
いが，未だ研究途上のものが多いため本稿では多くを述べることをしなかった．ディスク上の分布の理論と実際についてや，外れ値もしくはモデルに強く影響を与える観測値の検出に有効な影響診断法の開発等についても機会があれば稿を改めて紹介できればと考えている．
謝辞

原稿は，2 名の査読者による有益なコメントに基づいて改訂されました．ここに両査読者に感謝します．

参 考 文 献Abe, T. and Pewsey, A.（2011）. Sine-skewed circular distributions, Statistical Papers, 52, 683–707.
Abe, T., Pewsey, A. and Shimizu, K.（2009）. On Papakonstantinou’s extension of the cardioid distribution, Statistics & Probability Letters, 79, 2138–2147.
Abe, T., Shimizu, K. and Pewsey, A.（2010）. Symmetric unimodal models for directional data motivated by inverse stereographic projection, Journal of the Japan Statistical Society, 40, 45–61.
Abe, T., Shimizu, K., Kuuluvainen, T. and Aakala, T.（2012a）. Sine-skewed axial distributions with
an application for fallen tree data, Environmental and Ecological Statistics, 19, 295–307.
Abe, T., Kubota, Y., Shimatani, K., Aakara, T. and Kuuluvainen, T.（2012b）. Circular distributionsof fallen logs as an indicator of forest disturbance regimes, Ecological Indicators, 18, 559–566.
Arnold, B. C. and SenGupta, A.（2006）. Probability distributions and statistical inference for axial
data, Environmental and Ecological Statistics, 13, 271–285.
Arnold, B. C. and SenGupta, A.（2011）. Models for axial data, Advances in Directional and Linear
Statistics (eds. M. T. Wells and A. SenGupta), 1–9, Springer-Verlag, Berlin.
Batschelet, E.（1981）. Circular Statistics in Biology, Academic Press, London.
Cartwright, D. E.（1963）. The use of directional spectra in studying the output of a wave recorder
on a moving ship, Ocean Wave Spectra, 203–218, Prentice Hall, Englewood Cliﬀs.
Downs, T. D. and Mardia, K. V.（2002）. Circular regression, Biometrika, 89, 683–697.
Fisher, N. I.（1993）. Statistical Analysis of Circular Data, Cambridge University Press, Cambridge.
Fisher, N. I. and Lee, A. J.（1992）. Regression models for an angular response, Biometrics, 48,
665–677.
Gill, J. and Hangartner, D.（2010）. Circular data in political science and how to handle it, Political
Analysis, 18, 316–336.
Jammalamadaka, S. R. and SenGupta, A.（2001）. Topics in Circular Statistics, World Scientiﬁc,
Singapore.
Johnson, R. A. and Wehrly, T. E.（1978）. Some angular-linear distributions and related regression
models, Journal of the American Statistical Association, 73, 602–606.
Jones, M. C. and Pewsey, A.（2005）. A family of symmetric distributions on the circle, Journal of
the American Statistical Association, 100, 1422–1428.
Kato, S. and Jones, M. C.（2010）. A family of distributions on the circle with links to, and applications
arising from, Möbius transformation, Journal of the American Statistical Association, 105,
249–262.
Kato, S., Shimizu, K. and Shieh, G. S.（2008）. A circular-circular regression model, Statistica Sinica,
18, 633–645.
Mardia, K. V. and Jupp, P. E.（2000）. Directional Statistics, Wiley, Chichester.

304

統計数理

第 61 巻

第 2 号 2013

Papakonstantinou, V.（1979）. Beiträge zur zirkulären Statistik, Ph.D. Dissertation, University of
Zurich, Switzerland.
Pewsey, A.（2002）. Testing circular symmetry, Canadian Journal of Statistics, 30, 591–600.
Pewsey, A., Shimizu, K. and de la Cruz, R.（2011）. On an extension of the von Mises distribution
due to Batschelet, Journal of Applied Statistics, 38, 1073–1085.
SenGupta, A. and Ugwuowo, F. I.（2006）. Asymmetric circular-linear multivariate regression models
with applications to environmental data, Environmental and Ecological Statistics, 13, 299–309.
島谷健一郎 編（2011）. 生物群集の多様性と統計数理，統計数理研究所共同研究リポート，No. 263,
177–211.
島谷健一郎（2012）. 『フィールドデータによる統計モデリングと AIC，ISM シリーズ：進化する統計数理 2』，近代科学社，東京．
Shimizu, K. and Iida, K.（2002）. Pearson type VII distributions on spheres, Communications in
Statistics–Theory and Methods, 31, 513–526.
Wang, M.-Z. and Shimizu, K.（2012）. On applying Möbius transformation to cardioid random variables, Statistical Methodology, 9, 604–614.
Zar, J. H.（2010）. Biostatistical Analysis, ﬁfth ed., Prentice Hall, Upper Saddle River.

Proceedings of the Institute of Statistical Mathematics Vol. 61, No. 2, 289–305 (2013)305

Use of Directional Statistics in Environmental Science
Kunio Shimizu1,2 and Minzhen Wang3
1 Department of Mathematics, Keio University
2 Visiting Professor, The Institute of Statistical Mathematics
3 School of Fundamental Science and Technology, Keio University;

Now at The Institute of Statistical Mathematics

Wind direction is a typical and important angular variable in environmental science.
Sometimes observations include wind speed and ozone concentration together with winddirection. This article is aimed at obtaining fundamental results in directional statistics
that deal with modeling and analysis of data, including angular observations and a review
of recent studies in this area.
The diﬀerence between the sample mean direction for angular data and the sample
mean for linear data, and that between the mean resultant length and the variance, areemphasized. We use wind direction data available at the site of the National Institutefor Environmental Studies to illustrate circular plots and calculations of basic statistics. Pewsey’s test for symmetry, a ﬂexible symmetric angular distribution by Jones andPewsey, and the method of sine-skewing are presented, and an illustration of ﬁtting angular distributions using the wind direction data is shown. Papakonstantinou and Batschelet
distributions, which have the ﬂat-topped and sharply-peaked properties, are introduced.
The article also discusses circular-linear regression models derived from distributions on
the cylinder, circular-circular regression models, and circular-circular structural models.

Key words: Angular regression, angular structural model, distributions on the cylinder, distributions
on the torus, ﬁtting angular distributions.