移動関数k-means法の提案とその応用

要　旨

本論文では，関数データにおける部分的な特徴を考慮した「移動関数k-means法」を提案し，実際の適用例を通してその有効性を示す．具体的には，従来の関数データクラスタリングのように定義域を固定せずに，部分的な区間を逐次的に移動させながら関数クラスタリングを行う．また本手法を，福島県における空間線量率を測定した大規模センシングデータに適用し，従来手法では検出しにくいようなクラスターの同定に対して有効であることを示す．

1. はじめに

近年のデバイス技術の著しい進歩によって，以前までは得られなかったような非常に高い頻度で，大量のデータを得られる場面が増えている．特に時間に依存したデータが大量に収集され，その解析が社会的に望まれている．

このようなデータは，時間帯により特異な挙動を示すことがあり，従来の統計解析手法のように各データ点に基づいて解析を行うのではなく，データの集まりを一つの要素と考えて解析する手法が提案されている．その中で，連続的なデータを「関数」と見て，関数として表されたデータの集まりを一つの要素として扱う関数データ解析（Functional Data Analysis, FDA）が注目されている．FDA は，1990 年代にRamsay らによって提唱され，これまでの成果をまとめた成書も出版されている(Ramsay & Silverman 2002, 2005; Ramsay, Hooker & Graves 2009)．関数データ解析では，従来の多変量解析法を関数に対応した手法に拡張して用いることが多く，これまでに関数回帰分析や関数主成分分析，関数クラスター分析など多くの有用な手法が開発・応用されている(Ramsay et al., 2005)．

ところで，これらの関数データ解析法においては，定義域をあらかじめ固定している．しかし，実際のデータ解析では部分的な定義域に限定された特徴を導出するなど，定義域を予め仮定しない関数データを分析する場合も考えられる．例えば，Mizuta (2004) は定義域を固定しない階層的なクラスタリングについて提案した．

本論文では，定義域を考慮した関数データに対する非階層的なクラスタリング法を提案する．また実際の適用例として，福島県における空間線量率を測定したセンシングデータを用いて，本手法の有効性を示す．

2. 関数クラスタリングについて

関数クラスタリングは，似ている関数同士をグループに分類するための一連の手法であり，これまでに多くの手法が開発されている(Ferraty & Vieu, 2006; Jacques & Preda, 2013)．

以下ではn 個のデータを{x_i(t); i = 1, 2, . . . , n}, t ∈ [a, b]; a, b ∈ R とおき，n 個の個体を，互いに素な集合C₁, C₂, . . . , C_K に分割することを考える．

2.1. 階層的関数クラスタリング

i とj の間の非類似度を定義することで，それに基づく手法として階層的関数クラスタリングが提案されている．例えば，

により非類似度を定義することで，通常の階層的クラスター解析法のアルゴリズムを用いることができる．

非階層的関数クラスタリング

代表的な非階層的手法として関数k-means 法がある．これはクラスターC_k の重心となる関数をとおき，目的関数

を最大化するようなクラスターを求める手法である．この他に，関数ファジィクラスタリングなども提案されている(Tokushige, Yadohisa & Inada, 2007)．

3. 移動関数k-means 法の提案

経時的なデータを扱う分野では，適切な期間に対する移動平均を用いて季節変動などを調整する手法がしばしば用いられ，それにより大きなトレンドを見つけることができる．移動平均と同様に，期間に対応して関数データの定義域を柔軟に扱うことは重要である．

一般に，2 節で述べた関数クラスタリングでは，関数の定義域は固定して扱われる．そこで，定義域を固定せずに，部分的な定義域を逐次移動させながら関数データをクラスタリングする移動関数k-means 法（Moving Functional k-means Clustering）を提案する．

r に依存して変化する定義域をI(r) = [r, r + w] ⊆ [a, b]; r,w ∈ R とおき，その上で定義される関数データ{x₁(t), x₂(t), . . . , x_n(t)}; t ∈ I (r)を考える．I(r) ごとに，関数k-means 法を適用する．すなわち，n 個の個体をr に依存したK 個の互いに素な集合に分割する．

で重心となる関数をとおいて，目的関数

表1 各クラスターに属する個体の推移

を最大にするクラスターを各r ごとに計算し，を得る．

各r ごとに独立してクラスタリングを行うため，各r でつけられたラベルC^(r) には一意性が担保されないが，できるだけ整合するよう再ラベリングする．具体的には，以下の指標

が最大となるようにC^(s) のラベルを換えて，新たなクラスターラベルを得る．

実際の関数データにおける計算では，等間隔に離散化した関数データを扱う．すなわちt をt_j ; j = 1, 2, . . .,M として離散化した関数データx_i(t_j) を考える．I(rl) = [rl, rl + w]; l =1, 2, . . . , L; r₁ = t₁; r_L + w ≤ t_M とおく．w = t_M とした場合には従来の関数k-means 法の結果と一致する．w < t_M とした場合，定義域の移動幅r_l+1 −r_l はクラスラベルの連続的な変化に対しての解像度であり，w を小さくするほど解像度は高くなりクラスラベルの変化を詳しく分析できる．本手法では，t_j と同じ解像度（l = j,L = M）とする．

(4) はこのとき，

となる．

(5) の具体例として1757 個体を3 つのクラスターに分ける場合を考える．表1 はI(r_m) からI(r_m+1) にかけての，各クラスターに属する個体の推移を示している．例えば，1 行目はクラスター1 に属する555 個体のうち415 個体が同じクラスター残り，140 個体はクラスター3 に移ったことを表している．指標(6) の最大化は，表の対角成分の和（合計を除く）を最大化するように列を入れ替える操作に相当し，この例では，最大化の基準を満たしていることがわかる．あとは，対応するクラスターのラベルに付け替えればよい．

4. 数値実験

提案手法を人工データに適用して，その有効性を確かめる．まず，150 個の関数データx_i(t); i = 1, 2, . . . , 150 を生成させる．本実験ではB-スプライン基底関数を用いた．全体の定義域をI = [1,5000] とし，節点（ノット）を等間隔に50 個おき，また次数は5 とした．基底関数Φ₁(t), Φ₂(t), . . . , Φ₂₄(t) に対する係数行列{c_ib}; i = 1, 2, . . . , 150; b = 1, 2, . . . , 24 を正規分布c_ib～N(μ_b, 1) により生成させた（表2）．t をt_j = j; j = 1, 2, . . . , 5000 として離散化し，関数データ

を得る．図1 は生成した関数データを示しており，係数を発生させた分布ごとに色分けされている．

図2 に従来の関数k-means 法を適用した結果を示す．また，図3a, b に提案手法を，それぞれw = 400 およびw = 200 として適用した結果を示す．t_j = 1000 からt_j = 3000 付近の間で，クラスターの構成が大きく変化するが，提案手法ではその変化を捉えていることがわかる．一方，従来の関数k-means 法では，部分的な区間におけるクラスター構成の変化は捉えられていない．定義域の幅であるw については，w = 200 はw = 400 に比べ，局所的な変化を反映したクラスター構成の変化を表せている（例えばt_j = 3000 付近）．

表2 実験で用いた分布のパラメータ

図1 実験に用いた関数データ

図2 従来の関数k-means 法

図3a 提案手法（w = 400）

図3b 提案手法（w = 200）

5. 実際例

提案した移動関数クラスタリング法を，福島県の空間線量率を測定したセンシングデータに対して適用する．本データは原子力規制委員会によって公開されているモニタリングデータ（http://radioactivity.nsr.go.jp/map/ja/）である．

このデータは，2011 年3 月11 日に起きた東北地方太平洋沖地震の影響で，東京電力福島第一原子力発電所で発生した事故による環境への影響を調べるためであり，全国3921 箇所に順次設置されたセンサー（線量計）を用いて，空間線量率（μSv/h）を10 分間隔で測定している．福島県には最も多くのセンサーが設置され（2013 年8 月15 日時点で2380 箇所），各地域ごとにラベル（会津，南会津，いわき，相双，県北，県南，県中）がつけられている．

本論文では，福島県の会津地域およびいわき市を合わせた計1757 箇所のセンサー（会津地域：468 箇所，県中地域：800 箇所，いわき市：458 箇所）において，2012 年10 月1 日から2013 年5 月1 日までに測定されたデータを用いる．測定は10 分ごとになされているため，欠損データがなければ各センサーで30529 個のデータが得られている．ここではそれらの原データをとおく．なお，この期間中2 週間以上にわたり欠損しているセンサーは114 か所あり，1757 箇所には含めずに予め除外した．

関数データ解析の前処理として，B-スプラインによる平滑化を施し，データが欠損している時刻（センサーごとに約50～500 個程度）では補間を行った．平滑化基準には一般化クロスバリデーションを用い(Hastie et al., 2009)，また各関数データに対して，によって標準化を施した．

に対し，提案手法を適用する．定義域の大きさw は2 週間（6×24×14 = 2016）とし，10 分ごとに定義域を移動させてI(r_j) = [r_j, r_j +w](j =1, 2, . . . , 30529) とした．また10 分ごとに再ラベリングを行った．

図4a～7a に従来の関数k-means 法（クラスター数，K = 2）を適用した結果を示す．図4b～7b に，提案手法である移動関数k-means 法（K = 2, 3, 4, 5）を適用した結果を示す．図4a では，12 月から4 月にかけて大きく変化している関数とそうではない関数とのクラスターの特徴を捉えているが，10 月から12 月および4 月から5 月におけるクラスターの特徴は解釈できなかった．一方，移動関数k-means 法（図4b）では，12 月から4 月では従来の関数k-means 法と同様なクラスターを検出しているが，10 月から12 月および4 月から5 月の期間においても，解析に用いた全期間における平均値（x_i(t_j) = 0）よりも高いクラスターと低いクラスターに大きく分かれていることが分かる．図4a と図4b から，12 月から4 月におけるクラスターの構成と10 月から12月および4 月から5 月におけるクラスターの構成は異なるものであると考えられる．さらに12月から4 月の期間について詳細に調べると，この期間周辺では，会津地域は積雪量が多くなる反面，いわき市ではほとんどない．また，県中地域は会津地域といわき市の中間的な状況となっている（気象庁データベース，2012 年12 月から2013 年5 月にかけての「最深積雪」データを参照：http://www.data.jma.go.jp/obd/stats/etrn/）．クラスターの構成がこの期間で大きく変化した要因に，積雪の影響が考えられる．

図4a 従来の関数k-means 法を適用した結果（K = 2）

図4b 移動関数k-means 法を適用した結果（K = 2）

図5a 従来の関数k-means 法を適用した結果（K = 3）

図5b 移動関数k-means 法を適用した結果（K = 3）

図6a 従来の関数k-means 法を適用した結果（K = 4）

図6b 移動関数k-means 法を適用した結果（K = 4）

図7a 従来の関数k-means 法を適用した結果（K = 5）

図7b 移動関数k-means 法を適用した結果（K = 5）

同様に，K が3 以上の場合でも12 月から5 月におけるクラスターの構成が変化している様子が表れている．例えば，先に示した表1 は，図5b 中の2013-3～2013-4 付近におけるラベルの急激な変化を説明しており，具体的には2014 年3 月23 日前後で，赤に対応したクラスターから140個体が緑に対応したクラスターに移り，さらに青に対応したクラスターから緑に対応したクラスターに300 個体が移っており，結果として緑に対応したクラスターのサイズが186 個体から415個体となっていることがわかる．他のクラスターラベルの大きな変化についても，表1 と同様な表を作ることで，クラスター構成の大きな変化を確かめることができる．

図8, 9 にクラスターと地理的な対応関係を示す．横軸は経度（東経），縦軸は緯度（北緯）であり，図中の三角形（紫）は東京電力福島第一原子力発電所を指している．図8a, 図8b, 図8c はそれぞれ，K = 2 における降雪前の2012 年11 月と降雪後の2013 年1 月，雪解け期の2013 年4 月におけるクラスターを示している．積雪の多い会津地域（青色）と，積雪が少ないか，あるいは積雪がないような地域とでクラスターが分かれている（表3a）．同様に，図9 にK = 3の場合を示す．図9a では主に二つのクラスターに分かれているが，降雪後の2013 年1 月から3 つのクラスターに分かれはじめ，2013 年4 月では3 つのクラスターに大きく分かれている（表3b）．前半では主に2 つのクラスターから構成されているのに対して，後半では主に3 つのクラスターに分かれており，前半と後半とでは大きく異なっていることが見られる．これは，雪解けタイミングが地域によって異なるためであると解釈できる．ただし，クラスターが急激に変わるところや異常に低い外れ値の部分など，解釈が困難な部分も残されている．

表3a 図8におけるクラスターごとの地域別センサー数

表3b 図9 におけるクラスターごとの地域別センサー数

6. おわりに

本論文では，関数クラスタリングにおいて定義域を逐次的に移動させて解析する移動関数k-means法を提案した．従来の関数クラスタリングにおいて，定義域についてはあまり考慮されていなかったが，移動平均における期間のように，部分的な区間ごとの関数データに着目して分析することで，従来のクラスタリングでは検出がしにくいようなクラスターの特徴を検出することができる．本手法を適用した，福島県における空間線量率を測定したセンシングデータの例では，積雪のある冬期間にクラスターの構成が変化する様子を捉えることができた．

適用例では，定義域の幅w として2 週間をとったが，他にも1 ヶ月，3 か月なども考えられる．w を大きく（小さく）とる程，結果として大域的（局所的）な変動を捉えたクラスターが得られる．本論文では，積雪の前後における，線量率の局所的な変化が，非常に大きいことを考慮して，2 週間とした．また，r_j についても10 分間隔以外の場合も考えられるが，詳細なクラスター構成の変化を分析するために，最小の時間間隔である10 分とした．

関連した研究として，Mizuta (2004), Mizuta and Kato (2007), Tokushige et al. (2007) では，個体間の非類似度を関数s_ij(t) として定義し，各t ごとにクラスタリングを行う関数クラスタリング法や多次元尺度構成法を行う関数多次元尺度構成法を提案している．それらに対して本手法では，関数データの部分的な変化を捉えることができる点に特徴がある．

また，本手法において，グラフィカルな表現も結果を解釈において非常に重要となる．特にアニメーションを活用することで，定義域を移動させていく過程におけるクラスターの連続的な変化の様子をより効果的に表現することができる．

図8a 提案手法によるクラスターと地図との対応関係（K = 2, 2012 年11 月）

図8b 提案手法によるクラスターと地図との対応関係（K = 2, 2013 年1 月）

図8c 提案手法によるクラスターと地図との対応関係（K = 2, 2013 年4 月）

図9a 提案手法によるクラスターと地図との対応関係（K = 3, 2012 年11 月）

図9b 提案手法によるクラスターと地図との対応関係（K = 3, 2013 年1 月）

図9c 提案手法によるクラスターと地図との対応関係（K = 3, 2013 年4 月）

最近の関数データ解析法の研究では，多次元の関数データに対する新たな手法も提案されている(Berrendero, Justel & Svarc, 2011; Jacques et al., 2014)．適用例で用いたデータの分析では，空間線量率の他に気温や降雨量，積雪量などの複数事象に基づく関数データも合わせて，より詳細な分析を進めていくことが必要で，その手段として多次元関数データの解析手法は重要と考えられる．今後の課題として，本手法の拡張を含め，多次元関数データ解析における具体的手法の開発をしていきたい．

References

• Berrendero, J.R., Justel, A. & Svarc, M. (2011). Principal components for multivariate functional data. Computational Statistics & Data Analysis, 55(9), 2619-2634, Elsevier.
• Ferraty, F. & Vieu, P. (2006). Nonparametric Functional Data Analysis. Springer.
• Hastie, T., Tibshirani, R. & Friedman, J. (2009). The Elements of Statistical Learning (2nd ed.). Springer.
• Jacques, J. & Preda, C. (2013). Functional data clustering: a survey. Advances in Data Analysis and Classification, 1-25, Springer.
• Jacques, J. & Preda, C. (2014). Model-based clustering for multivariate functional data. Computational Statistics & Data Analysis, 71, 92-106, Elsevier.
• Mizuta, M. (2004). Clustering methods for functional data. Proceedings in Computational Statistics, 1503-1510, Springer.
• Mizuta, M. & Kato, J. (2007). Functional Data Analysis and Its Application. Lecture Notes in Artificial Intelligence 4481, 228-235, Springer.
• Ramsay, J. O., Hooker, G. & Graves, S. (2009). Functional Data Analysis with R and MATLAB. Springer.
• Ramsay, J. O. & Silverman, B. W. (2002). Applied Functional Data Analysis. Springer.
• Ramsay, J. O. & Silverman, B. W. (2005). Functional Data Analysis (2nd ed.). Springer.
• Tokushige, S., Yadohisa, H. & Inada, K. (2007). Crisp and fuzzy k-means clustering algorithms for multivariate functional data. Computational Statistics, 22(1), 1-16, Springer.
• 松井佑介・南　弘征・水田正弘 (2013). 放射性物質の分布状況等調査データベースを用いた空間線量率分布に関する考察, 2013年度統計関連学会連合大会講演報告集, 145.