T - バイオスーパーコンピューティング研究会

バイオスーパーコンピューティング東北2014@仙台
平成26年6月2日(月曜)
アウトライン
データ同化による
エミュレータ・デザイン学の創設
何度も講演
 内挿と外挿
 データ同化
 実例：地震音波とバイオイメージング
樋口知之（情報・システム研究機構統計数理研究所）




1
1/28
2/28
内挿と外挿問題
データ
エミュレータ
スパース回帰（気象予測）
オンライン機械学習
GPR: ガウス過程回帰
いろいろなビジネス展開が可能
データ無しの領域
移流と拡散
＋
クラウドソーシング
フォワード計算モデル
現況を捉える認識力
＝予測能力
スマホ（とGPS）
t  t  t
u
＋
内挿
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 equation
xt  f t ( xt 1 )
代入計算
xt  f t ( xt 1 , v t )
(time varying)
vt
システムモデル
Stochastic simulation model
State Vector
x
 cx 2  
t
ISM
 ξ 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 )
気象・海洋のデータ
同化の枠組み
t
yt  Ht xt  wt , wt～N (0, Robs )
 t : sampling time of observatio ns
t : simulation
time step
Missing
value (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 forecasting
p ( xt | y1:t ) or　p ( xt | y1:T )
状態（時刻依存）変数の推定
逐次（オンライン）型：集団の時間
発展を追う。つまり、Swarm Filter
非逐次（オフライン）型：ベスト初期値をもつパス
を求める
代表例： EnKF (Ensemble Kalman
Filter)
代表例：４次元変分法（Adjoint法）
■ Parameter Estimation in a simulation model
p ( | y1:T )  p ( y1:T |  )  p ( )
  x0
9/28
モデルパラメータの推定
4DVar
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 by
of 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 acoustic
waves causing oscillations in the ionosphere at the altitude of
a few hundred kilometers.
seismic fault
seismic wave country, so our governmentMagnitude
of 9.0
Japan is an earthquake-prone
has been
promoting many researches related to earthquake prediction and
ISM
earthquake disaster prevention.
When each earthquake occurred, many micro-barometer sensors
located at distances of hundreds or thousands of kilometers from the
epicenter observed “the sound of the earthquake” in the form of
atmospheric pressure perturbations.
GPS Observation
(Website of Kyoto University)
observation
seismic 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 right
dvits mechanism

1
before a cell division, but actually
  2 v is yetp
Navier-Stokes equation
unknown.
dt 

We assume that proteins “myosin” in the cell
D wall drive
physically
such 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 Image
Velocimetry) 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

t
t 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
xt1 を多数シミュレートしyi,tの
アンサンブル予測を行なう
スパース（疎性を利用した）最適化
例：多変量回帰
シミュレーション変数の計算結果（１０万～１００万次元）
同化の結果を使う場合もある
y , x 
サンプル
減らす１.
y  X a e
=
a1  22  a2  22
12
22
+
新NP問題 (N

 P)
0
2
a   min y  X a   a
a

a2
a1  a2  1
(2,2)
a1
2
(1,0)
min . y  Xa subj. to a  
ドリア
a
17/28
18/28
減らす２.
ストリーミング計算
行列、テンソル分解
顧客
ID付きPOSデータ
データは現場で産出。移動は困難。
Edge-Heavy Data (PFI, 丸山＠ISM)
日付
来店データ
ID
i軸：顧客ID
k 軸：日時
だいたい似かよるように
代数演算のみで分解
 パラメータ（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
t
時刻の
生データ
時刻 +1 の
生データ
圧縮
圧縮
影響
影響
職業的、年齢的
j 軸:商品
尤度関数
xの空間
20/28
1
エミュレータ：
オンライン計算 in 機械学習
シミュレーションパラメータを入力とする出力値を模倣
シミュレーションの結果
巨大次元
パラメータ（a）オンライン学習法
小次元
x
FOBOS: Forward Backward Splitting
RDA: Regularized Dual Averaging
ある地点の降雨量（温度）
パラメータ（特徴）変数抽出
スカラー出力
z
スパース化
ブラックボックス
気温の全球平均
二酸化炭素排出量
一般の学習器
さきほどの話（スパース回帰）
min( Xa)  (a)
w
yt , xt Tt1
時刻 t の正誤（ある，なし）
時刻 t の画像データ
エミュレー
ションの実施
要は統計モデル

1
2
at 1  min ayt , xt  a aa , a  (a) 
a  at 
t
a
2t


m(z)
パラメータ入力
z 
N
i 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
2
E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 
y
gn n1
m(z)
g
y2
回帰係数
g1
g2
y1
g が離散だったら
gn1
m(n)
gn
g0
yn
z
g  g (z)
y2
z2
共通の分散
カーネル関数：入力値
（パラメータ）間の距離
g
y1
z2
0
n
t
yn  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 2dH
m (z)  m(z)  (Y  M )K 1t(z)
Y T   y1,, yN 
事後分布
p( 2 ) 
データと統計モデルの差
規格化ファクター
tT (z)  k(z, z1),, k(z, z N )
y1

ガンマ分布
自由度 N   のt分布
1
k ( z, z ' )  k (z, z ' )  t ( z ) K t (z ' )
y2
T
エミュレータ
サンプル数
回帰モデル部分
1
2
, p(H ) 1
  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 
2
2
    (z  z' )
   2 exp  
2
2


2
2
L
S




データ同化とエミュレータ
E y g(z)  m(z)
こちらの構造ではない！
データ同化
 12 exp  
g  g(z)
a. 行列分解
b. ローカル化
g  g(z)
 L  6 S
機械学習を用いた内挿法
z
Gaussian Processes for Regression:A Quick Introduction
M. Ebden, August 2008
27/28
I. スパース回帰
II. オンライン学習
III. カーネル法
z
 ( z  z' )2 
2
2
2
   2 exp  2 sin ( ( z  z ' ))    ( z  z ' )
2 2 

 12 exp  


横幹連会誌「横幹」Vol.8,No1（2014年4月15日）
ミニ特集「数理科学の展開とその体制」
『統計数理の誕生とその広がり』樋口知之
28/28

Download Report