新領域融合ＰＪから：融合研究の苦労と成果３「データ同化」（PD 統計数理研究所 中野純司）
２０１４年１０月１７日
情報・システム研究機構シンポジウム２０１４

アウトライン
∗はじめに∗データ同化とは∗具体的な成果∗研究基盤整備∗融合研究の難しさとおもしろさアウトライン
∗はじめに∗データ同化とは∗具体的な成果∗研究基盤整備∗融合研究の難しさとおもしろさプロジェクトの変遷

新領域融合プロジェクト 統計数理基盤：
データ同化による複雑システムの定量的理解と計測デザイン

データ中心科学 リサーチコモンズ モデリング・解析基盤整備：
データ同化・シミュレーション支援技術

本プロジェクトの目的∗ 広範囲な分野の複雑なシステムを対象とする研究の具体的な問題解決に向けた共同研究をとおして、シミュレーションと
データ解析の作業を一体化する手法の高度化と一般化を実現する。あわせて、未適用分野の発掘を戦略的にすすめ、個別科学を横断的につなぐ新しい学問領域を創る。

∗ データに基づいて学術研究や社会における課題解決を目指
す研究コミュニティに対して、データ基盤整備事業によって提供するデータを高度活用し、課題の背後にある構造を可視化・モデリングし、シミュレーションなどを通じて、適切な課題解決を促進する研究環境を提供する

本プロジェクトの研究体制∗ ３チームで構成∗ 数理・計算チーム：
理論・ツールの開発∗ モデリングチーム：
既存の多量データのモデリング∗ データデザインチーム： データ取得からの解析∗ データ同化研究開発センター（統数研）が運営のとりまとめ∗ データ同化を用いる具体的明確な研究テーマ∗ 地震による微気圧変動∗ 感染症伝搬∗ グローバルＭＨＤシミュレーション∗ 大気レーダーの３次元高分解能観測∗ 動物細胞の分裂の動態モデル∗ 「ほ乳類体節形成と精子形成」のモデル化に向けたデータ取得

アウトライン
∗はじめに∗データ同化とは∗具体的な成果∗研究基盤整備∗融合研究の難しさとおもしろさ

・科学研究推進に 演繹的
は４輪が必要！

シミュレー
・前輪の軸が今ま ション

統計数理の実績

データ同化

帰納的

大量データ処理でなかった
・データ同化がそ
の軸となる

実験

理論
科学の駆動力本プロジェクトで解明する現象
・オーロラ発生時の磁気圏の状況
・大気レーダーの３次元高分解能観測
・動物細胞の分裂の動態モデル
・ほ乳類体節形成と精子形成 など

データ同化とは？（１）データ同化とは？（２）• 気象学・海洋学の分野で発達 （1990年代中頃から）• 物理数値シミュレーションモデルと実際の観測を統合する手法– シミュレーションのみでは適切に現実の物理現象を再現できない• シミュレーションモデルには，モデルの不完全性や境界条件が正確にはわからないなどの不確かさが存在• （たとえば）正確な気象予測には適切な初期条件の構成が必要– 観測データは物理的・社会的制約により得られる情報に限界がある

観測データを用い数値シミュレーション内の変数を修正＝データ同化シミュレーションモデルと誤差∂x= cx 2 + ∂txt = f t ( xt −1 )xt = f t ( xt −1 , vt )

データ同化の定式化• データ同化は非線形非ガウス状態空間表現での状態 xn の推定：(システムモデル)

xt = f t ( xt −1 , vt | x0 )(観測モデル)yt = H t xt + wt

または

yt = ht ( xt ) + wt

∗ xt : シミュレーションモデル内の全ての変数∗ yt : 観測される全ての変数∗ vt : モデルの不確かさ (境界条件, …)∗ wt : 観測誤差
∗ vt , wt : （通常）ガウス分布，
超高次元

x0: 初期条件

104~106 yt : 102~105 dim( xt ) >> dim( yt )xt : （1万～100万）（100～10万）分布の表現：モンテカルロ近似p( xt | y1:t −1 ), p( xt | y1:t ), p( xt | y1:T )True distribution

モンテカルロ近似
実現値の集合で分布を表現する

N: 粒子数[]≡ [x , x , , , x ]p( xt | y1:t −1 ) ≅ X t|t −1 ≡ x , x , , , xp( xt | y1:t ) ≅ X t|t(1)
t |t −1(2)
t |t −1(N )t |t −1

(1)
t |t

(2)
t |t

(N )
t |t

13

逐次データ同化の手法∗ 逐次データ同化においては，アンサンブルカルマンフィ
ルタ(EnKF)がよく用いられる．
∗ 条件付分布を実現値の集合（アンサンブル）として保持∗ フィルタリングにはカルマンフィルタの式を利用∗ われわれは粒子フィルタ（Particle Filter (PF))も応用する．
∗ 実現値集合を用いる点でEnKFと同様
1 N
p ( xt | y1:t −1 ) ≅ ∑ δ ( xt − xt(|it )−1 )N i =1
1 N
p ( xt | y1:t ) ≅ ∑ δ ( xt − xt(|it ) )N i =1i 番目の実現値

{x }
{x }

N

(i )
t |t −1 i =1
(i ) Nt |t i =1

(i )
t |t

x
時刻

使用した観測の最後の時刻

アンサンブルカルマンフィルタ
の手続き
1.
2.

t←１ とし，初期状態の実現値集合 {x }以下を行う(i ) N0|0 i =1

A) (一期先予測)
i.
ii.

{ }

(i ) Nを作る．
ti =1

システムノイズの集合 v
(i )
(i )
(i )
各 i に対して， xt |t −1 = f t ( xt −1|t −1 , vt ) を計算する．

B) (フィルタリング)
i.
ii.
iii.

を作る．

{ }

(i ) Nti =1を作る．

観測ノイズの実現値集合 w
分散共分散行列 Vˆt|t −1 , Rˆ t ，カルマンゲイン K̂ t を計算する．
各 i に対して， xt(|it ) = xt(|it )−1 + Kˆ t ( yt − H t xt(|it )−1 + wt( i ) ) を計算する．

3. t=T ならば停止．それ以外は，t←t+1 として２．に戻る．

アンサンブルカルマンフィルタの手続き（続）ただし

1 N (i )xˆt|t −1 = ∑ xt|t −1
N i =1

1 N (i )xˆt|t = ∑ xt|t
N i =1

N
1
(i )
(i )
ˆ
ˆ
( xt|t −1 − xt|t −1 )( xt|t −1 − xˆt|t −1 )'Vt|t −1 =∑N − 1 i =1N1(i ) (i )ˆRt =wt wt ' （所与の Rt で代替可能）∑
N − 1 i =1

−1
ˆ
ˆ
ˆ
ˆ
K t = Vt|t −1 H t ' ( H tVt|t −1 H t '+ Rt )

EnKF vs PF∗ どちらも状態変数の分布を多数の粒子の集合（アンサンブル）で表すことは同じ∗ 異なるのは、時間の経過ごとの粒子の決め方∗ アンサンブルカルマンフィルタ∗ それぞれの粒子を観測に合うような方向に寄せる∗ 粒子フィルタ
∗ 観測への当てはまりの悪い粒子を破棄し、減った分だけ当て
はまりの良い粒子を複製して増やす∗ なお、当てはまりの良さは尤度ではかるアウトライン
∗はじめに∗データ同化とは∗具体的な成果∗研究基盤整備∗融合研究の難しさとおもしろさ

地震音波検出を目的とした微気圧観測計画（京都大学・極地研究所・電気通信大学との共同研究）

微気圧観測点設置を計画中東日本大震災に伴う地震波・地震音波伝搬シミュレーション

東日本大震災による地震音波強度や伝搬特性を評価した

地震音波を利用した津波早期警戒システム
概念図 地震音波が津波よりもずっと速い伝播速度で沿岸に到達することを利用した津波早期警戒システムを提案。日本経済新聞2012年1月26日夕
刊第一面で取り上げられた。
 菅平宇宙電波観測所と潮岬風力実験所に微気圧観測装置を設置し、それぞれ2011年12月および2013年3月から測定を開始。
 津波早期警戒システムの基盤となる地震音波データ同化システムを開発。震源に関するモデルパラメータと各時刻の地震音波を高速に推定できることを検証した。その研究結果は、日本地球惑星科学関連学会
2012年大会のハイライト論文に選出された。

微気圧観測拠点京都大学 潮岬風力実験所

地震音波データ同化計算

電気通信大学 宇宙電波観測所

統計数理研究所設置PCクラスタ

システムの詳細と動作検証実験
MCMCによる推定

津波早期警戒システム

パラメータの収斂過程推定されたパラメータ分布

burn-inモデルパラメータ
断層破
震源位置
壊伝搬
速度

θ = (C V

地震音波の状態推定
風によるマグニ
チュード 時間遅延効
果

m1:N

観測ノイズ

t0 σ 1:2K )事後分布
p (θ | δ p1:obs
K ,1:T ) ∝TKp (θ ) ∏∏k 1=t 1 =N 
exp  −  δ pkobs− ∑ mnδ p sim
 2σ k2 tnL,n , k ,t − t 0 −

 

n 1
=
2πσ k2( N −1)V 

1

赤: 推定値
黒: 疑似観測データ(真値+ノイズ)

感染症数理モデルによる感染症介入政策の評価
 都市における個人の行動を考慮に入れた大規
模感染症シミュレーション(人口:122万人)を多数の計算機を接続して実行する、並列計算によって実現した。
 ワクチンを会社員などの活動的な集団に優先的に接種することで集団免疫を高め、高リスク者を含む住人全体の感染割合を低下させ得ることを示唆する結果を得た。
シミュレータの模式図。鉄道で結ばれた複数の地区からなる都市を扱う。それぞれの住民は各自のスケジュールによって行動し、その結果として形成される集団内でのインフルエンザの確率的伝播をシミュレーションする。

# of patients

infected
exposed 本研究の成果は、PLOS ONEに掲載され、毎日新聞平成25年9月19日夕刊で取り上げられた。

シミュレータの可視化。
エージェント・シミュレー
ションの内部で起っている
状況を2次元グラフで捉えるのは困難である。また、感染の進行を非専門家にもわかりやすく伝えることが感染症予防に役立てるために重要である。このような背景からシミュレータの開発と並行して、可視化システムの開発にも取り組んできた。

シミュレーション例

利用事例: ワクチン優先接種対象群の選択∗∗

通常ウィルス上陸前に当該地域の全住民に接種することは不可能である。そのため、どのような集団に優先接種すれば、ワクチンを受けられなかった集団内での感染を最小化できるかが焦点になる。活動性の高い会社員に集中投与すると、非接種者
感染割合を最小化できるというシミュレーション結果を得た。
最適な優先接種のもとでは、感染者数のするどいピークが消滅し、シーズン全体に分散していることが確認できる。これは、感染路が消滅しない程度に伝達確率が低く抑えられたことによる。

優先接種対象群と感染割合

最適な選択の下での感染者数の推移

地球磁気圏グローバルMHDシミュレーションモデルによるデータ同化の研究
才田聡子（ROIS）、門倉昭（極地研）、藤田茂（気象大）、田中高史（九州大学）、他

地球磁気圏グローバルMHDシミュレーションモデル

地球・電
離圏

磁気圏内部
境界

MHD (MagnetoHydroDynamics)：電磁流体力学
磁気圏モデル内部境界～接続（任意性）～電離圏モデル（３Re（地球半径））（１Re)
データ同化対象パラメータ
同化対象データ（候補）１．電離圏電位分布： SuperDARNレーダー、衛星観測
２．電離圏電流： 地上磁場、衛星観測
３．電離圏電気伝導度： EISCATレーダー、オーロラ観測＜全体計画＞
平成22年度非同化モデルの検証
データ獲得戦略立案
平成23年度
データ同化初期実験
データ獲得・収集
平成24年度
初期実験の解析と評価
データ獲得・適用
平成２5年度
最適値推定・評価（続き）
データ獲得・適用・結果検討
平成２6年度同化手法の実装・改善
データ獲得・適用・結果検討
平成２7年度同化手法の実装・改善・運用
成果公開オーロラ発生時の共役点の移動データの取得について

図１（下）：共役点移動の様子
図２（右）：緯度・経度方向の共役点の移動

共役点の移動方向

共役点の移動方向
北半球（上）と南半球（下）で反対方向に共役点が移動する

地球磁気圏グローバルMHDシミュレーションモデルによるデータ同化の研究
才田聡子（ROIS）、門倉昭（極地研）、藤田茂（気象大）、田中高史（九州大学）、他

平成２５年度の成果
SuperDARNレーダーネットワーク観測データとシミュレーション結果の比較
電離圏電位分布とプラズマ対流分布

SuperDARNレーダー観測データ

シミュレーション結果

大気レーダーの３次元高分解能観測と同化の研究
西村耕司,中村卓司(極地研),佐藤亨(京大),佐藤薫(東大),堤雅基,冨川喜弘,山内恭(極地研)

１．新アルゴリズムによる鉛直風速推定
精度向上を実証– 同化データのための物理量として精度の高い鉛直風速を得るため、ノＭＵレーダー
ルム拘束付ＤＣＭＰ法を用いる手法を開発。これにより推定精度が大幅に向上することを京都大学MUレーダーによる実験で実証した。
南極昭和基地大型大気（ＰＡＮＳＹ）レーダー

２．PANSYレーダーとの組み合わせによる超高精度風速推定• PANSYでは積雪対策のためにアンテナを分散化しており、細いメ
インビームと高いサイドローブを有する。数値最適化により静的にサイドローブを抑圧した上で、上記手法と組み合わせることにより、mm/sオーダーの理論精度が得られた。これは地球大循環
における極域での鉛直速度と同程度であり全球モデルの同化への利用が可能性が見えてきた。

高精度推定アルゴリズムと分散アレイによる超高精度風速推定
NC-DCMPアルゴリズムにより適応最適ビーム制御

ＰＡＮＳＹのアンテ
ナ配列
他の大型
レーダーと
異なる分
散型アレ
イ配列
高サイドローブは精度上NG

アニーリング法
により最適化さらにNishimura et al., JAOT, 2012.
推定精度
NC=1.0(not adapted)
1.2

1.5

２つの技術の相乗効果によりmm/sオーダーの超高精度推定が可能にデータデザインチーム （遺伝研中心）
目標：データ同化手法を発生細胞生物学分野の問題に適用し、同分野における有力な手法として確立する
これまで、バイオインフォマティクスの分野で、データ同化の手法が適用されてきた。
例えば、遺伝子発現量の網羅的計測の時間的変化の背後にある遺伝子発現制御ネットワークの推定など。
Yoshida R, Nagasaki M, Yamaguchi R, Imoto S, Miyano S, Higuchi T. Bayesian learning of biological pathways on genomic data
assimilation. Bioinformatics. 24 (22),2592-2601 (2008)

データデザインチームは、発生生物学、細胞生物学の分野においてデータ同化法、あるいは関連手法を適用し、生物学と統計学の融合による、同分野でのブレークスルーを目指す。

動物細胞の分裂、細胞質流動

ほ乳類の体節形成、精子形成

細胞分裂（線虫）：力学的変化

主な参加研究者：木村（NIG）、樋口(ISM)、
小山(元TRIC研究員)【目的】 細胞分裂時の細胞形状の時系列データから細胞表層の硬さを推定する。
【結果】

細胞分裂時の細胞形状(座標、曲率)

データの可視化ソフトウェア開発：
H22融合研究シーズ探索（長尾代表）

表層
硬い

表層
硬い

① 細胞分裂時の細胞表層の硬さを推定した→分裂面が柔らかくなっていることを支持する（従来、分裂面の軟化は提案されてきたが、定量的な証拠はなかった）② 細胞分裂に異常をきたす遺伝子変異体では、表層の軟化が阻害されていた→推定の信頼度を支持する、表層の軟化の機能的意義を支持する
Koyama H., Umeda T., Nakamura K., Higuchi T., and Kimura A. A high-resolution shape fitting demonstrated
equatorial cell surface softening during cytokinesis. Submitted.

細胞質流動（線虫）

主な参加研究者：木村、庭山（NIG）、長尾(ISM)【目的】 細胞質流動の流速分布から、流動の原動力の分布を推定する
細胞質流動とその定量化
anterior

posterior

従来：手作業（試行錯誤）でのパラメータ推定
Niwayama, Shinohara, Kimura A, PNAS, 2011
流体シミュレーション(MPS法)

cortical myosin
GFP::VIT-2 (Yolk granule)

PIV methodによる流速の定量化

データ同化を用いたパラメータ推定→力の分布が高い精度で推定できる→力を生み出す機構についての考察が可能に-0.1 μm/s
map>

0.1 μm/s

Preliminary trial

Kimura, A. The 49th Annual Meeting of the Biophysical Society of Japan. [S2F-3], 2011.9.16

線虫初期胚における細胞内流動データ同化（遺伝学研究所との共同研究）

priorp (θ )

likelihoodp( y |θ )posteriorp (θ | y )Model parameters• velocity of the cell wall
how to design uwall ( L, t )• viscosity coefficient• density

Particle filter• search model parameters thatmaximize the likelihood• observed velocity field (2D)MCMC
• estimate posterior
distribution

Niwayama et al. [PNAS 2011]データ同化によって推定された細胞内流動の可視化

粒子フィルタ法によるデータ同化アルゴリズムを実装し、理研RICCや京速コンピュータ「京」上で、ジョブ投入実験を実施した。

体節形成（マウス）：体節時計の同調機構

主な参加研究者：相賀、木村、菅原（NIG）、
大久保(元TRIC研究員)• 周期的な体節の形成機構に関する研究では、理論的解析もすすんでおり、数値シミュレーションも多く行われている。→データ同化の適用に適したモデル系？
• 体節は、各細胞にある「体節時計」の時間的振動が、空間的な縞模様に翻訳されることによって形成される。
• きちんと「縞模様」ができるためには、多数の細胞の「体節時計」が同調していなければならない。
• 同調機構の理解を進めるために、正常細胞と遺伝子機能欠損細胞をランダムに混合した「キメラ」胚を作製した。複雑なキメラ胚実験を解釈するには、データ同化を利用した数理モデルでの実験観察結果の再現が効果的である可能性を考えた。
• 現時点では、データ同化（パラメータ推定）は行っていない。データ同化に向けてシミュレーションモデルを作成したところ、キメラ胚実験の解釈はパラメータ値に左右されない定性的な性質で可能であることがわかってきた。

体節時計の振動と縞模様の形成 (Oginuma…Saga, 2010)

実験データ：キメラ率vs同調度

遺伝子発現量Deltaキメラ

前
領域

体節
後

複雑なキメラ個体の結果を解釈したい

同調度の違いを説明するモデル

Lfngキメラアウトライン
∗はじめに∗データ同化とは∗具体的な成果∗研究基盤整備∗融合研究の難しさとおもしろさデータ同化のための4次元可視化ソフトウェア
バーチャルオーロラ

細胞内流動

地震音波伝搬（岩手・宮城内陸地震）

地震音波伝搬（東日本大震災）

超高並列計算機のためのデータ同化アルゴリズムの高度化~トーラス型ネットワークトポロジーに適した粒子交換アルゴリズム~
京速コンピュータ「京」は3Dトーラス型のネットワークトポロジーを持つため、粒子フィルタを適用するためには、効率的に粒子交換を実施するためのアルゴリズムが必要。

Lorenz 96モデルに対し、 統計数理研究所が所有するCray XT6m
（2Dトーラス）で実験
32x32

32x32

 The total number of particles is 4194304 in the left panel and 262144 in the right panel.
 If the number of particles for each node is sufficient, local degeneracy problem can begenerally avoided.

数値実験結果
本年度上半期は、各ノードあたりの粒子数を一定としてノード数を増加させること
で，計算速度を落とすことなく推定性能を向上できることを実験的に確認した．
PEs

Particles

RMSE

Elapsedtimes

Effective
sample size

Effective
number ofnodes

Effective
sample size
for a node

4x4

65536

1.61

00m54s22467.1

15.6

1440.18x82621441.18

01m00s90762.2

59.6

1523.316x1610485761.11

01m04s315334.6

203.8

1547.1

32x32

4194304

1.09

01m07s993516.8

638.6

1548.5

 If the number of particles for each node is fixed, the increase of the number ofnodes would result in the increase of the effective number of particles.
 The increase of the number of nodes would make sense if the aim is toimprove the accuracy or to address high-dimensional problems.
Λ

Λ

λ

=
S −∑ Ω(kµ ) log Ω(kµ ) + ∑ Ω(kµ ) (−∑ ωk( µ ,ν ) log ωk( µ ,ν ) )
=
µ 1

=
µ 1=
ν 1

並列計算ソフトウェア∗ CloCK-TiME (Cloud Computing Kernel for Time-series Modeling Engine)∗ http://sheep.ism.ac.jp/CloCK-TiME/∗ Web 上で多変量時系列データを粒子フィルタ法によって分析するシステム∗ （現在，統計解析システムRのパッケージに移植中）∗ Rhpc (An R package for High-Performance Computing)∗ 統計解析システムRでMPIを使って並列計算を効率的に行うためのパッケージ∗ スーパーコンピュータ上でのRの利用，データ同化の実行を念頭に置いて開発∗ “The Comprehensive R Archive Network”の世界中のミラーサイトより入手できる∗ 日本では例えば http://cran.ism.ac.jp

データ同化の普及活動∗

樋口知之著, 「予測にいかす統計モデリングの基本―ベイズ統計入門から応用まで」, 講談社, 2011.

∗

樋口知之編著, 上野玄太，中野慎也，中村和幸，吉田亮，シリーズ＜予測と発見の科学＞6「データ同化入門－次世代のシミュレーション技術―」, 朝倉書店, 2011.

新分野への適用など∗ データ同化の多くの新分野への適応の試み∗ 地球・宇宙科学∗ 多数（特にこれまであまりデータ同化が使われていない領域）∗ 生命科学∗ 細胞∗ 個体形成∗ ゲノム科学∗ 医学∗ 感染症の予防政策∗ 若手人材の育成∗ 融合研究における∗ 統計数理研と極地研，遺伝研との密な人材交流∗ 育った若手人材の外部への供給

アウトライン
∗はじめに∗データ同化とは∗具体的な成果∗研究基盤整備∗融合研究の難しさとおもしろさ

融合研究の難しさ∗ 目的意識の違い∗ （データ同化を含む）統計数理は「データ」を扱うための学問∗ 不確定性を含むデータから有用な情報を取り出すことがその目的∗ ただし、特定のデータだけではなく汎用性があることが必要∗ 各分野の研究者は自分の問題が解決すればよい∗ すでに結論を（自分の中で）決めている場合がある∗ 統計的な推論法がまわりくどいと感じることが多い∗ 統計的な考え方への理解不足∗ 日本では統計学がきちんと教えられていない∗ データを軽んじる∗ 白黒をはっきりとつけたがる∗ 日本には統計学科がない（先進国では特異な状況）∗ 分野が異なると（専門）用語が異なる∗ 意思疎通に苦労する＝時間がかかる

融合研究のおもしろさ∗ 学問の幅が広がる∗ 統計数理側： データの重要性・価値の再認識∗ 個別分野側： 統計的な考え方・ツールを知る∗ コミュニケーション能力の向上∗ 「井の中の蛙」「専門馬鹿」からの脱却∗ 社会へのアピール能力も向上∗ 人間的な成長ができる（？）∗ お互いを尊重することが必要∗ 単なる「下請け」「丸投げ」「いいとこ取り」ではうまくいかない今後の課題と期待∗ 融合研究は学問的には重要、でも参入障壁がある∗ それを突破するためには内的・外的な助けが必要∗ 本機構の融合プロジェクトは非常に有効な外的補助だった∗ このような補助がなんらかの形で続くことを期待する∗ データ同化は応用の広い統計手法、でも参入障壁がある∗ われわれでその参入障壁を低くしたい∗ 多くの分野への応用∗ スパコンやクラウドを有効利用できる使い易いソフトウェアの開発と配布∗ データ中心科学の有用なツールとしてのデータ同化がより広く使えることでビッグデータの科学的な有効活用が可能となると期待している