バイオスーパーコンピューティング東北2014@仙台
平成26年6月2日(月曜)

アウトライン

データ同化によるエミュレータ・デザイン学の創設

何度も講演 内挿と外挿 データ同化 実例： 地震音波とバイオイメージング

樋口知之 （情報・システム研究機構 統計数理研究所）

 エミュレータ スパース回帰（気象予測）
 オンライン機械学習
 GPR: ガウス過程回帰
1

1/28

2/28

内挿と外挿問題データ

いろいろなビジネス展開が可能

データ無しの領域

移流と拡散 ＋ クラウドソーシング ＝ 予測能力
フォワード計算モデル

現況を捉える認識力 スマホ（とGPS）t  t  tu＋

内挿D
情報の不確実性（多人数からのレポート：多様なノイズの混在）

3/28

4/28

データ同化の目的：気象・海洋学の観点からつなぐ：データ同化[1] 予報を行うための最適な初期条件を求める。これは既に、現業の天気予報で実用化されていることである。
[2] シミュレーションモデルを構成する際の最適な境界条件を求める。連成現
象を取り扱う際の適応的な境界条件設定もこの作業に含まれる。

データ同化[3] スケールが異なるシミュレーションモデル間の橋渡しを行うスキーム内に含まれる諸パラメータの最適な値を求める。経験的に与えられるモデル内のパラメータ値の検証も一つの具体例である。
[4] シミュレーション（物理）モデルにもとづいた、観測されていない時間・空間点における観測値の補間を行う。この作業は再解析データセットの生成とも呼ばれる。このデータセットから新しい科学的発見をもくろむ。 ダウンスケーリング(2011年9月刊行）

経験[5] 時間・経費を節約できる効率的な観測システムを構築するための仮想観測ネットワークシミュレーション実験や感度解析を行う。
（参考文献： 蒲地 他、「統計数理」、54(2), 223-245, 2006.)
6

5/28

6/28

データ同化と一般状態空間モデル

システムモデルとしてのシミュレーションモデル
（simplified meteorological model around Japan）

State Vector（Simulation variables）

PDE : Partial differential equationxt  f t ( xt 1 )代入計算

xt  f t ( xt 1 , v t )(time varying)vt

システムモデル

Stochastic simulation modelState Vectorx cx 2  tISM ξ 1,t    ξ m ,t x t   ξ m  1,t 


  
 ξ M ,t  θx t  f t ( x t  1 , v t ),v t ～ p ( v | θ sys )y t  h t ( x t , w t ),w t ～ p ( w | θ obs )気象・海洋のデータ
同化の枠組みtyt  Ht xt  wt , wt～N (0, Robs ) t : sampling time of observatio ns
t : simulation
time step
Missingvalue (Outlier) t  1  t

Observation model
Measurement model

観測モデル

t
time integration

7/28

8/28

逐次（アンサンブル）vs. 非逐次（最適パス）

Our interest in a statistical inference by DA

データ同化のイメージ
データ空間にシミュレーション解を射影した時

ht ( x t )■ State Estimation and forecastingp ( xt | y1:t )orp ( xt | y1:T )状態（時刻依存）変数の推定

逐次（オンライン）型： 集団の時間発展を追う。つまり、Swarm Filter

非逐次（オフライン）型：ベスト初期値をもつパスを求める

代表例： EnKF (Ensemble Kalman
Filter)

代表例： ４次元変分法（Adjoint法）■ Parameter Estimation in a simulation modelp ( | y1:T )  p ( y1:T |  )  p ( )モデルパラメータの推定  x 0 4DVar
9/28

10/28

Japan, An Earthquake-Prone Country (Nagao et al., 2011)

Ionospheric Disturbances Excited by the 2011 Great East Japan EQ
Nagao

acoustic wave

Observation (GPS-TEC)

Nagao

From website of Kyoto Univ.
Magnitude byof 7.2
Application to analyze an acoustic wave which is generated
earthquake followed
417km by Tsunami.
seismic fault

In the case of the Great East Japan Earthquake, the dense
array of GPS receivers in Japan clearly observed the acousticwaves causing oscillations in the ionosphere at the altitude ofa few hundred kilometers.

seismic fault

seismic wave country, so our governmentMagnitudeof 9.0Japan is an earthquake-prone
has been
promoting many researches related to earthquake prediction andISMearthquake disaster prevention.

When each earthquake occurred, many micro-barometer sensors
located at distances of hundreds or thousands of kilometers from theepicenter observed “the sound of the earthquake” in the form ofatmospheric pressure perturbations.

GPS Observation(Website of Kyoto University)

observationseismic wave

acoustic wave

ionosphere

In the case of the Iwate-Miyagi Nairiku Earthquake, for example,
micro-barometer sensors 417 km away from the epicenter observed
acoustic waves excited by the earthquake.
Because acoustic velocity is much faster than the propagation
velocity of tsunami, such acoustic signals may be useful in predicting
the magnitude of approaching.

11/28

Assimilation
12/28

Video (10min) uploaded to YouTube
Introduction to R&D Center for Data Assimilation

Data Assimilation on Intracellular Fluid in C.elegans embryo(Collaboration with Dr. Kimuta at NIG)

Nagao

Observation

Assimilation

Research Projects

Methods:■ Applications of Sequential Data Assimilation Methods－ Meteorology－ Space Science－ Life Science■ Advanced Monte Carlo algorithm and Applications■ Statistical Algorithms Dedicated to Massively Parallel Computers

Data assimilation works powerfully also in the cell biology.

MPS method

A large dynamic intra-cellular flow occurs in a fertile egg rightdvits mechanism
1
before a cell division, but actually  2 v is yetpNavier-Stokes equation
unknown.
dt We assume that proteins “myosin” in the cellD wall drive
physicallysuch intra-cellularand v we perform data
mass conservation
law flow,Dt that solves
assimilation by combining a numerical simulation
simultaneously these physical equations and the observation
obtained by the image analyses.

PIV (Particle ImageVelocimetry) imaging

Hardware:■ Hardware Random Number Generators

Computational Services:

Niwayama et al. [PNAS 2011]■ Statistical Analysis System for Massively Parallel Computing Environment■ Cloud Computing Services

13

13/28

機能のモデル化：エミュレータ
高度情報社会におけるユニバーサルな研究課題の表現形■ Visualization Software

14/28

機能のモデル化：“見よう見まね”を科学する

達成度：予測可能性，再現性などのパーフォーマンスからの評価
基礎方程式群

実体的に積み上げ式でモデル化・“見よう見まね”のプロセスを加速する。

？

未知の出力

ー体系化されていない研究分野において有効未知の入力・“見よう見まね”による完成の域がお手本を超える。

複雑なシステム入力

出力

過去の入力データ

未知の入力

日本人は伝統的に“贋作”を心から嫌う傾向が強いが、同じく漢字文化圏にある中国人で心ある人は、過去の文物に目を転じるとき、それがたとえ”贋作”であろうと判断できても、”本物”を超える出来栄えであるならば、自分の心と目を満足させるため、”本物”以上にその”贋作”を尊重して手に入れると伝えられている。

過去の出力データ機能自体を模倣する（入出力関係を近似する）

数理モデル
エミュレータ15/28

ー贋作が“本物”を超える

未知の出力“贋作”と”本物”の、どちらが本物か実は分からない
16/28

気象予測：気象エミュレータ
観測データベクトル（数千～１万次元）ttt 1~ yi,t  Hi xt 1  et
スパース学習xt1  yi,t予測

シミュレーション変数ベクトルの新しい入力

yi  a0  a1xi,1  a2 xi,2  aP xi,P  ei (i  1,, N) y1,t    yt   yi,t     yL,t  T

ある地点の降雨量（温度）（スパース）線形回帰

例：多変量回帰

シミュレーション変数の計算結果（１０万～１００万次元）同化の結果を使う場合もあるy , x 

サンプル減らす１. スパース（疎性を利用した）最適化

xt1 を多数シミュレートしyi,tのアンサンブル予測を行なう

y  X a e=a1  22  a2  22  1
12
22+新NP問題 (N  P)

02a   min y  X a   aa

a2
a1  a2  1(2,2)

a12(1,0)min . y  Xa subj. to a  ドリア

a

17/28

18/28

減らす２.

ストリーミング計算

行列、テンソル分解
顧客

ID付きPOSデータ

データは現場で産出。移動は困難。
Edge-Heavy Data (PFI, 丸山＠ISM)

日付
来店データ

ID

i軸：顧客IDk 軸：日時

だいたい似かよるように代数演算のみで分解 パラメータ（w）オンライン学習法
M2Mの増大↓ スパース学習
FOBOS: Forward Backward Splitting
トラフィックの大渋滞
RDA: Regularized Dual Averaging
ベイズの定理と情報循環

p( y | x )  p( x )p( x | y)  p( y | x )  p( x )p( y)

週末的、月末的yijk =0,1,2…

確率分布『バラバラを癖で束ねる』を実現
19/28

確率分布

t

時刻 +1 の生データ

圧縮

圧縮
影響
影響

職業的、年齢的j 軸:商品尤度関数t

時刻 の生データ

xの空間
20/28

エミュレータ：

オンライン計算 in 機械学習

シミュレーションパラメータを入力とする出力値を模倣
シミュレーションの結果
巨大次元パラメータ（a）オンライン学習法

小次元x

FOBOS: Forward Backward Splitting
RDA: Regularized Dual Averaging

ある地点の降雨量（温度）パラメータ（特徴）変数抽出

スカラー出力z

スパース化

ブラックボックス

気温の全球平均二酸化炭素排出量

一般の学習器

さきほどの話（スパース回帰）min( Xa)  (a)w

yt , xt Tt1

時刻 t の正誤（ある，なし）
時刻 t の画像データ
エミュレー
ションの実施要は統計モデル12
at 1  min ayt , xt  a aa , a  (a) a  at t
a2t

m(z)

パラメータ入力z Ni i 1
シミュレー
ションの実施降下法

モンテカルロ実験
実験計画

スパース性

Fused LASSO (Tibshirani et al., 2005)のオンライン版
21/28

結果（予測値） m(zi )iN1

補正と補間m (z)モデルの改良結果（実験データと思えば良い） yi iN1

22/28

エミュレータ：統計モデル（要は、線形回帰モデル＋GP）

エミュレータ：時系列のトレンドモデルとほぼ同じ

システムモデルGaussian Process：g は無限次元（連続変数）

p(g)  p(g(z))  N(m(z), k(z, z)  2 )g  g(z)12Eg(z)  m(z)  H    Cov g ( z ), g ( z ' )   k ( z , z ' )  Covariance function
統計モデルは所与z
カーネル関数k (z, z ' )  exp  (z  z ' )T R(z  z ' )観測モデル

yi  gi  e, e～N(0, 2 ) (i 1,, N)
事後分布

g 連続変数p(g | Y )  p(Y | g) p(g)Y   y1,, yN ygn n1m(z)

gy2

回帰係数

g1

g2y1

g が離散だったら

gn1m(n)gn
g0ynzg  g (z)

y2z2

共通の分散

カーネル関数：入力値（パラメータ）間の距離

gy1z20ntyn  gn  e, e～N(0, 2 )

g1

g2

z1

gn  m(n)  gn1  vn , vn～N(0,~2 )

m(z)

T

23/28

g(z)

システム

カーネル関数は定数値zk (n, n  1)  exp  (1)T R(1)

z1
24/28

エミュレータ：補正（キャリブレーション）と補間
簡単のため

エミュレータ：フルベイズへ
簡単のため

事後分布

2
 0 の場合
2

p(g(z) | Y )～N(m (z), k (z, z' ) )


M T  m(z1),, m(zN )

2

2
 0 の場合
2

入力値とのずれの程度

p(g(z) | Y , 2 , H )～N(m (z), k  (z, z' ) 2 )p(g | Y )   p(g | Y , 2 , H ) p( 2 ) p(H )d 2dHm (z)  m(z)  (Y  M )K 1t(z)Y T   y1,, yN 

事後分布1p( 2 )  2 , p(H ) 1

データと統計モデルの差規格化ファクターt (z)  k(z, z1),, k(z, z N )

ガンマ分布

Ty1



自由度 N   のt分布

1k ( z, z ' )  k (z, z ' )  t ( z ) K t (z ' )y2Tエミュレータサンプル数

回帰モデル部分  dim(z)  1

データがもたらす情報により不確実性は必ず減少する k(z1, z1)  k(z1, zN ) K    
k(zN , z1)  k(zN , zN )

参考文献
J. Sacks et al., “Design and analysis of computer experiments,” Statistical Science, 1989.
M. C. Kennedy and A. O’hagan, “Bayesian calibration of computer models,” J. Roy. Statist. Soc. Ser. B, 2001.
中野、樋口、“地球科学におけるシミュレーションとビッグデータ―データ同化とエミュレーション―、” 信学会会報、to appear
Rasmussen, C. and C. Williams, Gaussian Processes for Machine Learning, MIT Press, 2006.

グラム行列

25/28

26/28

エミュレータ：カーネル関数の設計（問題毎に与える）

Covy g(z), g(z' )   2k (z, z' )   2 (z  z' ) ( z  z' )2  ( z  z' )2 22    (z  z' )   2 exp  2222LS


データ同化とエミュレータ

E y g(z)  m(z)こちらの構造ではない！

データ同化

 12 exp  g  g(z)a. 行列分解b. ローカル化g  g(z) L  6 S

機械学習を用いた内挿法

zGaussian Processes for Regression:A Quick Introduction
M. Ebden, August 2008

27/28

I. スパース回帰
II. オンライン学習
III. カーネル法z ( z  z' )2 222
   2 exp  2 sin ( ( z  z ' ))    ( z  z ' )2 2  12 exp  

横幹連会誌「横幹」Vol.8,No1（2014年4月15日）
ミニ特集「数理科学の展開とその体制」『統計数理の誕生とその広がり』 樋口知之

28/28