科研費シンポジウム「統計学と機械学習における数理とモデリング」報告書
日時：２０１６年２月２１日（日）～２月２２日（月）場所：東京工業大学大岡山キャンパス西八号館 W８３３号室
開催責任者：鈴木大慈・野村俊一（東京工業大学）
本シンポジウムは下記の助成を受け開催された。
科学研究費 基盤研究（A）15H01678 「大規模複雑データの理論と方法論の総合的研究」
研究代表者：青嶋誠（筑波大学）【プログラム】

１日目 ２０１６年２月２１日（日）13:25～13:30：オープニング
13:30～15:30：機械学習と計算量セッション
1. Nonlinear Feature Selection for Large and Ultra High-dimensional data
2. Bayesian Network の構造学習における分枝限定法の適用に関する考察
3. 部分モニタリング問題における漸近最適アルゴリズム

山田 誠（京都大学）

岩本 俊亮（大阪大学）

小宮山 純平（東京大学）

4. マルコフ連鎖モンテカルロ法の高次元漸近論 鎌谷 研吾（大阪大学）※病欠により講演はキャンセル
16:00～18:00：機械学習理論セッション
5. 変分鞍点近似によるベイズ推論とモデル選択

林 浩平（国立情報学研究所）

6. グラフ表現を用いた異常変数同定における一致推定量

原 聡（IBM）

7. Non-convex Compressed Sensing with the Sum-of-Squares Method
8. 教師付き学習における Barron and Cover 理論

相馬 輔（東京大学）

川喜田 雅則（九州大学）

２日目 ２０１６年２月２２日（月）10:00～12:00：統計理論セッション
9. 状態空間モデルを用いた時系列の周期成分分解と位相推定
10. 関数予測における漸近ミニマックス予測分布
11. 低次元多様体上の Data Sharpening

松田 孟留（東京大学）

矢野 恵佑（東京大学）

工藤 雅紀（島根大学）

12. 確率的交互方向乗数法とマルチクラスグラフ型正則化学習への応用 鈴木 大慈（東京工業大学）

13:30～15:30：ベイズセッション
13. 余震活動の確率予測：ベイズ統計を用いたアプローチ

近江 崇宏（東京大学）

14. 点過程モデルによる地震活動異常の解析と群発地震への応用

熊澤 貴雄（統計数理研究所）

15. 統計的言語モデルに基づく化学構造の特徴抽出と分子設計への応用 池端 久貴（総合研究大学院大学）
16. Bay-Bridge : local linear approximated bridge and its model building procedure via the Bayesianmodel

保科 架風（中央大学）

16:00～17:30：生体ログデータ分析セッション
17. 潜在的歩数増加パターンと体組成指標との統計的関連性についての分析

大橋 洸太郎（立教大学）

18. 生体ログデータを使用した、歩行パターンと体組成指標との関連を検討する分析方法の提案
真弓（慶應義塾大学）
19. 体型を示す新たな体組成指標の提案とそのダイナミクス
17:30～：クロージング

野村 俊一（東京工業大学）

渡辺Nonlinear Feature Selection for Large and Ultra High-dimensional
Data
Makoto Yamada
Kyoto University
myamada@kuicr.kyoto-u.ac.jp
January 18, 2016
Abstract
We propose nonlinear non-negative least-angle regression (N3 LARS), that can select non-redundant features from a large and ultra high-dimensional data in nonlinear way. A key advantage of N3 LARS is that
it is a convex method and can find a globally optimal solution. Moreover, N3 LARS can easily incorporate
with map-reduce frameworks such as Hadoop and Spark. Thus, with the help of distributed computing, a setof features can be efficiently selected from a large and ultra high-dimensional data. The effectiveness of the
proposed method is demonstrated through an ultra high-dimensional biology dataset.

1

Problem Formulation

Let X = [x1 , . . . , xn ] = [u1 , . . . , ud ]> ∈ Rd×n denotes the input data, y = [y1 , . . . , yn ]> ∈ Rn denotes the output
data or labels and suppose that we are given n independent and identically distributed (i.i.d.) paired samples{(xi , yi ) | xi ∈ X , yi ∈ Y, i = 1, . . . , n} drawn from a joint distribution with density p(x, y). Note, Y could beeither continuous (i.e., regression) or categorical (i.e., classification). For the purposes of this paper consider x
to be a dense vector.
The goal of supervised feature selection is to find m features (m < d) of input vector x that are responsible
for predicting output y.

2

Nonlinear Non-Negative LARS (N3 LARS)The optimization problem of N3 LARS is given asemin kL−α∈RddX

f(k) k2 , kαk1 ≤ η, s.t. α1 , . . . , αd ≥ 0,
αk K
Frob

(1)

k=1

f = K̄/kK̄kFro , L
e = L̄/kL̄kFro such that 1> K1
f n = 0, 1> L1
e n = 0, and kKk
f 2 = 1., K̄ (k) = ΓK (k) Γ
where Knn
Fro(k)and L̄ = ΓLΓ are centered Gram matrices, Ki,j = K(xk,i , xk,j ) and Li,j = L(yi , yj ) are Gram matrices,K(x, x0 ) and L(y, y 0 ) are kernel functions, Γ = In − n1 1n 1>
n is the centering matrix, In is the n-dimensional
identity matrix, 1n is the n-dimensional vector with all ones, η ≥ 0 is the regularization parameter, and k · kFrob
is the Frobenius norm.
It is worth pointing out the objective function of N3 LARS can be reframed as
1−2

d
X
k=1

αk NHSIC(uk , y) +

d
X
k,l=1

αk αl NHSIC(uk , ul ),fL)
e is the normalized version of Hilbert-Schmidt Independence Criterion (HSIC)
where NHSIC(u, y) = tr(K(Gretton et al., 2005) and NHSIC(y, y) = 1. Note that NHSIC always takes a non-negative value, and is zero if
and only if two random variables are statistically independent when a universal reproducing kernel (Steinwart,2001) such as the Gaussian kernel is used.
If the k-th feature uk has high dependence on output y, NHSIC(uk , y) becomes a large value and thus αkshould also be large. On the other hand, if uk and y are independent, NHSIC(uk , y) is close to zero; uk tends tobe not selected by the `1 -regularizer. Furthermore, if uk and ul are strongly dependent (i.e., redundant features),NHSIC(uk , ul ) is large and thus either of αk and αl tends to be zero. That is, non-redundant features that havea strong dependence on output y tend to be selected by the N3 LARS. Thus, N3 LARS can be regarded as awidely used minimum redundancy maximum relevance (mRMR) feature selection method (Peng et al., 2005).

2.1

Nyström Approximation for NHSIC

We start introducing kernels used in this paper.
For input x,we first normalize
the input x to have unit standard deviation, and then use the Gaussian kernel,

0 2
)
,
where
σx is the Gaussian kernel width. In regression cases (i.e., y ∈ R), we normalize
K(x, x0 ) = exp − (x−x
2σx2


0 2
)
an output y to have unit standard deviation, and then use the Gaussian kernel, L(y, y 0 ) = exp − (y−y, where2
2σyσy is the Gaussian
kernel width. In classification cases (i.e., y is categorical), we use the delta kernel for y,
1/ny if y = y 0 ,0L(y, y ) =where ny is the number of samples in class y.
0otherwise,
Nyström approximation: To reduce the computational cost of Gram matrices, we use the Nyström approximation for NHSIC asNHSIC(u, y) = tr((F > G)2 ),−1/2

F = ΓKnb Kbb−1/2

G = ΓLnb Lbb−1/2/(tr((Kbb−1/2/(tr((Lbb−1/2 2>Knb
Knb Kbb) ))1/4 ,−1/2 2L>
nb Lnb Lbb) ))1/4 ,where Knb ∈ Rn×b , [Knb ]ij = K(ui , ub,j ), Kbb ∈ Rb×b , [Kbb ]ij = K(ub,i , ub,j ), Lnb ∈ Rn×b , Lbb ∈ Rb×b , and b isthe number of basis function.

2.2

Distributed computation with Map-Reduce

The Nyström approximation is useful for large n. However, for ultra high-dimensional cases (e.g., d > 106 ),the computation cost of N3 LARS is still large. To deal with this issue, we separately compute ck in N3 LARS
in parallel with distributed computing such as Spark1 . Since we can independently compute ck , we can easily
speed up N3 LARS by simply increasing the number of servers.

References
Gretton, A., Bousquet, O., Smola, A., and Schölkopf, B. (2005). Measuring statistical dependence with Hilbert-Schmidt
norms. In ALT .
Peng, H., Long, F., and Ding, C. (2005). Feature selection based on mutual information: Criteria of max-dependency,
max-relevance, and min-redundancy. IEEE Transactions on Pattern Analysis and Machine Intelligence, 27, 1226–1237.
Song, L., Smola, A., Gretton, A., Bedo, J., and Borgwardt, K. (2012). Feature selection via dependence maximization.
JMLR, 13, 1393–1434.
Steinwart, I. (2001). On the influence of the kernel on the consistency of support vector machines. JMLR, 2, 67–93.
1 https://spark.apache.org/

BD スコアによるベイジアンネットワーク
構造学習の効率化
岩本 俊亮

鈴木 譲

大阪大学大学院理学研究科

本研究ではベイジアンネットワーク (BN) の構造学習における BD スコアの新たな上界を導き，構造学習にかかる計算量をこの上界を用いることで削減できることを示した．
離散型確率変数 X = (X1 , . . . , Xn ) について，各変数を頂点とする有向非巡回グラフ (Directed

Acyclic Graph, DAG)Bs を考える．このとき BN の構造 Bs は同時分布 P (X1 , . . . , Xn ) の条件付き分布による因数分解を与える．つまり，Πi を Xi の親 (依存対象の変数，DAG における Xi を指す矢∏n
印の元にある変数) の集合とすると P (X1 , . . . , Xn ) = i=1 P (Xi |Πi ) となる．Bp = {P (Xi |Πi )}ni=1とすると，BN は (Bs , Bp ) で定義することができる．ただし，条件付き分布による因数分解が同じものを表現する DAG は区別しない．

X = (X1 , . . . , Xn ) について，データ D = {C1 , . . . , Cm } を与える．このとき新たなデータ∑
C が生じる予測分布は P (C|D) =allBs P (C|Bs ) · P (Bs |D) となる．ベイズ的アプローチによる BN の構造学習とは事後確率 P (Bs |D) が最大となる Bs を探索すること．ある (Bs , Bp ) に∏ついて，ri = |Xi |，qi = |Πi | = Xj ∈Πi rj とし，Xi の k 番目の値を yk ，Πi の j 番目の組をwj で表す．また，θijk = P (Xi = yk |Πi = wj )，Θij = (θij1 , . . . , θijri )，Θi = (Θi1 , . . . , Θiqi )，ΘBs = (Θ1 , . . . , Θn ) とする．このとき，Θij がディリクレ分布であること，つまり確率密度関数が
′
∑ri
∏ri
′
Γ(Nij
)
Nijk
−1
′
′
′
ρ(Θij |Bs ) = ∏ri Γ(N
Nijk
，Nijk
はディリクレ分布のパラメー
= k=1
(Nij
′ ) ·
k=1 θijk
k=1

ijkタ) となることを仮定すると∫P (D, Bs )=P (Bs )= P (Bs ) ·P (D|ΘBs , Bs )ρ(ΘBs |Bs )dΘBs
qi
n ∏
∏

ri
′
′
∏
Γ(Nijk+ Nijk )Γ(Nij)·′ +N )
′
Γ(N
Γ(N
ij
ij
ijk )
i=1 j=1

(1)

k=1

∑ ri

′
k=1 Nijk ，Nijk は Xi = yk かつ Πi = wj となるデータの数である．この同時確率 P (D, Bs ) を BD(Bayesian Dirichlet) スコアという．また，式 (1) において Nijk = qiϵri とを得る．ここで Nij =するとき，同時確率 P (D, Bs ) を BDeu スコアという．ϵ は定数で，ESS(Equivalent Sample Size)という．事後確率 P (Bs |D) が最大となる Bs については同時確率 P (D, Bs ) も最大となるので，実際の構造学習では同時確率が最大となるような Bs を探索する．これを BD スコア (BDeu スコア)による構造学習という．
式 (1) について，ri′′∏
Γ(Nijk+ Nijk )Γ(Nij)·
Q(Xi , Πi ) =
′
′
Γ(N
Γ(N
+
N
)
ij
ij
ijk )
j=1
qi
∏

k=1

を局所 BD スコアという．局所 BD スコアについては次の定理が成り立つ．

定理 1 離散型確率変数 X について，ΠX , Π′X を X の親集合で ΠX ⊂ Π′X なるものとする．

Q(X, ΠX ) > Q(X, Π′X )(2)となるとき，Π′X を X の親集合とした DAG による BD スコアは最大とならない．

BD スコアによる構造学習のアルゴリズムとしては Dynamic Programming，Branch and Bound，
AStar などがあるが，これらは
1. 局所 BD スコアの計算
2. BD スコアを最大にする DAG の探索というステップからなる．1 にかかる計算量は n · 2n−1 であるが，定理 1 式 (2) 右辺の局所 BD ス
コアの代わりに上界を比べることで，特定の条件下では局所 BD スコアの計算量を減らすことができる．また，この上界は小さいほどよい．de Campos と Ji は局所 BDeu スコアに対して上界を導き [1]，鈴木はこの上界を局所 BD スコアに一般化した [2]．定理 2 はその上界である．
定理 2 X を与えられたデータの i 番目の確率変数，ΠXi を Xi の親集合とする．このとき，Q(Xi , ΠXi ) ≤qi′∏max Nijk
j=1

′
Nij

が成り立つ．
本研究ではこれよりも小さい上界を導いた．
主定理 X を与えられたデータの i 番目の確率変数，ΠXi を Xi の親集合とする．このとき，Q(Xi , ΠXi ) ≤

qi N∏
ij −1
′
∏
α + max Nijk
j=1 α=0′
α + Nijが成り立つ．
この上界により BDeu スコアによる構造学習において計算実験を行い，定理 2 の上界を用いた場合と比較したところ，効率化できることが確認できた．

参考文献[1] de Campos, C. P., & Ji, Q, ”Properties of Bayesian Dirichlet scores to learn Bayesian network
structures”, Fox, M., & Poole, D. (Eds.), AAAI, pages 431-436, (2010).
[2] Suzuki, J., ”Eﬃciently Learning Bayesian Network Structures Based on the B&B Strategy:
A Theoretical Analysis”, Advanced Methodologies for Bayesian Networks 2015, pages 1-14,(2015).

部分モニタリング問題における漸近最適アルゴリズム
小宮山純平（東京大学）

1

部分モニタリング問題
以下の定式化は，基本的に部分モニタリング問題の既存研究 [1] に準拠している。部分モニタリング問題は，繰り返しゲームとして定義される。各ラウンド t = 1, . . . , T に学習者 ( アルゴリズム）はアクション i(t) ∈ [N ] := {1, . . . , N }を選択し，また，結果 j(t) ∈ [M ] := {1, . . . , M } が未知の確率分布 p∗ ∈ PM ( PM は M 次元確率単体）から引かれる。最初のラウンドの開始前に，損失行列 L = (li,j ) ∈ RN ×M とフィードバック行列 H = (hi,j ) ∈ [A]N ×M （ここで [A] は有限個のシグナル集合）が学習者に与えられる。各ラウンドごとに，学習者は損失 li(t),j(t) ∈ R を被り，
シグナル hi(t),j(t) ∈ [A] を受け取る。ここで，結果や損失について学習者は知ることができず，学習者が知りえる
情報はシグナルが [A] のいずれであったかという情報のみである。つまり，結果に関する部分的な情報が（アクションに依存して）与えられる状況で，如何にして結果の分布 p∗ を学習するかという問題である。
学習者の目的は，T ラウンドを通した損失の最小化である。結果の分布に対して，損失の期待値を最小化するようなアクションを多く選択したい。このアクションは，∗i∗ = arg mini∈[N ] L⊤i p と書くことができる。ここで，Li は L の i 行目である. ここで，i∗ は unique だと仮定し，一般性を失わずに i∗ = 1 とする. ここで，∆i =⊤ ∗(Li − L1 ) p ∈ R と定義し，Ni (t) を t ラウンド目までにアクション i を選択した数とすると，損失の最小化問題は，以下の Regret（最適なアクションを選択し続け+||p*-C1c||M

C2

p* C4
C1
C5
C3

た場合とアルゴリズムの期待損失の差分）の最小化問題と等価である。

Regret(T ) =T∑∆i(t) =

t=1∑

図 1: N = 5, M = 3 での∆i Ni (T + 1),

セル分解での一例.

i∈[N ]

この Regret のうち，本研究ではラウンド数 T を十分大きくとったときの挙動に興味がある。
部分モニタリング問題での結果分布に関し，セル分解が定義できる。アクション i が最適である結果分布の集合

Ci = {q ∈ PM : ∀j̸=i (Li − Lj )⊤ q ≤ 0}.
を i の最適セルと呼ぶ。すべてのアクションの最適セルは閉凸多面体であり，結果の分布の仮説空間である PM の各セルへの分解（図 1）が定義できる。
最適なアクションが発見できることの必要十分条件は，２つのアクションのどちらが期待損失が大きいか，すべてのアクションから得られるフィードバックを総合して判断可能であることである。つまり，任意のアクションのペア i, j に対し，Li − Lj ∈ ⊕k∈[N ] ImSk⊤ である（大域観測可能（globally observable)）ことである。
これまで，大域観測可能な問題のうち一部の問題に関しては Regret が Θ(log T ) であることが知られている [2]。
本研究では，大域観測可能なすべての問題に関して Regret が Ω(log T ) であることを示す。さらに，その log T の定数を明らかにし，また漸近的に最適性能のアルゴリズムを提案する。全ての証明や一部の定理，実験詳細などは省略するが，これらに興味がある方は著者らの論文 [3] を参照されたい。
上記のように単純な問題設定であるが，オンライン広告のクリック確率の最適化や実験計画の最適化のモデルである確率的バンディット問題，検索エンジンのランキング最適化のモデルである効用ベース比較バンディット問
題，繰り返しオークションにおける価格最適化問題など，実用上広いクラスの問題が部分モニタリング問題に属することが分かっており，この問題に関する効率的なアルゴリズムには大きな価値があると考える。

6000
R(t): regret

4000
3000
2000

R(t): regret

RandomFeedExp3
CBP
BPM-LEAST
BPM-TS
PM-DMED
LB

5000

1000
0 0
10

101

102

103 104
t: round

105

106

350
300
250
200
150
100
50
0 0
10

(a) 繰り返しオークション問題

RandomFeedExp3
CBP
BPM-LEAST
BPM-TS
PM-DMED
LB-Theory

101

102

103 104
t: round105106(b) 確率的バンディット問題

図 2: Regret-round 片対数プロット．LB は Regret 漸近下界，Random は一様ランダムなアクション選択，FeedExp3,

CBP, BPM-LEAST, BPM-TS は既存手法である．

2

Regret 下限
任意の結果分布 p ∈ PM と a > 0 に関して E[Regret(T )] = o(T a ) を満たすアルゴリズムを強一致なアルゴリズ

ムと定義する。強一致なアルゴリズムは，全ての真に最適ではないアームが実は最適であったという仮説を各ラ
ウンド t に 1/t 程度の有意水準で棄却する必要がある。

Lemma 1. 強一致なアルゴリズムに関して以下が成立する。
∑∀q∈C1cE[Ni (T )]D(p∗i ∥Si q) ≥ log T − o(log T ),i∈[N ]

ここで，p∗i ，Si q は，真の p∗ および仮説 q における，アクション i を選択したときのシグナルの確率分布である。
∑
また，D(p∥q) = i (p)i log ((p)i /(q)i ) は２つの分布間の KL ダイバージェンスである。

Regret の最小化解は本稿では省略するが，上記の制約を満たしつつ Regret(T ) を最小化するようなアームの選択方法から，C ∗ log T （C ∗ は最適係数）の形で導出できる．

3

PM-DMED アルゴリズム
アームを引く回数に関して，Lemma 1 の制約を満たしつつ Regret(T ) を満たすという問題は，真の分布が分かっていれば半無限計画（LSIP）を解くことによって求められる．しかし，真の確率分布はもちろん学習者にとって未知である．本研究では，この問題を確率分布に関しての経験推定値をプラグインすることで解く partial monitoring

deterministic minimum empirical divergence (PM-DMED) アルゴリズムを提案する．詳細は [3] を参照されたい．
推定精度を上げるための項や，推定失敗からの復帰の項などいくつかの理論保証のための項を含めることによって，このアルゴリズムは T に関して下界と漸近一致する Regret を持つ．
繰り返しオークションおよび確率的バンディット問題におけるアルゴリズムの計算機シミュレーションでの経験的な性能を図 3 に示す．図から明らかであるように，既存手法の Regret が漸近下界よりかなり大きいのに対し，提案手法である PM-RMED は漸近に下界と一致する Regret の傾きを持つ．

参考文献[1] Gábor Bartók, Dean P. Foster, Dávid Pál, Alexander Rakhlin, and Csaba Szepesvári. Partial monitoring classification, regret bounds, and algorithms. Math. Oper. Res., 39(4):967–997, 2014.
[2] Gábor Bartók, Navid Zolghadr, and Csaba Szepesvári. An adaptive algorithm for finite stochastic partial
monitoring. In ICML, 2012.
[3] Junpei Komiyama, Junya Honda, and Hiroshi Nakagawa. Regret lower bound and optimal algorithm infinite stochastic partial monitoring. In NIPS, 2015.

マルコフ連鎖モンテカルロ法の高次元漸近論
大阪大学大学院基礎工学研究科，数理・データ科学教育研究センター

鎌谷研吾

マルコフ連鎖モンテカルロ (MCMC) 法のより良いデザインのための高次元漸近論の枠組みと，高次元で有効な MCMC 法を紹介する．MCMC 法の高次元漸近論は 1980 年台後半から始まり，ベイズ統計学においても理論だけでなくその実用性が注目され，今も活発に研究がなされている．

Pd (dx) = pd (x)dx は Rd 上の確率分布とする．最も基本的な MCMC 法であるランダムウォーク型メトロ
ポリス (RWM) 法は次の手続きで定義された．ある x ∈ Rd を出発点として以下を繰り返す：

{

x∗
x

←
←

x + w, w ∼ Nd (0, σ 2 Id ),x∗ with probabilty α(x, x∗ ){}pd (x∗ )α(x, x ) = min 1,.
pd (x)

ただし，∗

d}m が定義される．このマルコフ連鎖は Pd (dx) を不変分布とこのステップにより d 次元のマルコフ連鎖 {Xm

して持つように出来ていて，例えば pd (x) > 0 (x ∈ Rd ) であればエルゴード的になり，下記右辺の積分が存在する限り大数の法則

1 ∑
d
f (Xm
)→
M m=1
M

∫
f (x)Pd (dx) a.s. (M → ∞)(1)が成り立つ [7]．したがって右辺の積分の近似値として左辺の有限和を用いることが正当化されるわけである．
次元が高くなるに連れ，(1) の右辺の数値積分は極めて困難になり，MCMC 法の意義が出る．しかし，高次d}m の振る舞いを見
元での計算負荷の増大は MCMC 法でも不可避である．これは d → ∞ としたときの {Xm

ることでその一端を説明できる．このような次元の極限の議論の枠組みを与えるのが高次元漸近論である．高次元漸近論で様々なことがわかる．

1. Pd が正規分布に十分近い時は，採択率がおおよそ 23% 程度であるように σ > 0 を調整するのが良いと示唆される [5]．

2. Pd の裾の形状に依存するが，RWM 法よりも収束レートの良い MCMC 法が存在する [6]．
3. Pd の裾が軽ければ，MCMC の繰り返し回数は O(d) 回必要であり，裾が重い場合は O(d2 ) 回必要になる [3]．
本発表では高次元漸近論を使って，高次元でも有効な MCMC 法のデザインを考える．上述のように高次元では MCMC 法の繰り返し回数を長くしなければならないが，加えて一回の pd (x) の値の計算量も d に依存する．従って高次元で MCMC 法を適用する際は，計算量を抑える工夫が不可欠である．

RWM 法が困難なほどの高次元積分では，ラプラス変換を援用する方法や，並列化を用いる方法など様々な方法が存在する．本発表では少ない繰り返し回数で収束する MCMC 法を考える．近年，エルゴード的な提案
カーネルを使った preconditioned Crank-Nicolson (pCN; [1]) が理論的に注目されている．本発表では pCN法を改良した mixed preconditioned Crank-Nicolson (MpCN; [2]) 法を紹介する．MpCN 法は高次元漸近論に基づいてデザインされた有効な手法である：ρ ∈ (0, 1) とする．


 r
x∗

x

∼ Gamma(d/2, ∥x∥2 /2),
← ρ1/2 x + (1 − ρ)1/2 r−1/2 w, w ∼ Nd (0, σ 2 Id ),← x∗ with probabilty α(x, x∗ )

ただし，

}
{
pd (y)∥x∥−dα(x, y) = min 1,.
pd (x)∥y∥−d

この手法がなぜ有効であるか，どれくらい従来手法を改善するかを数値実験と高次元漸近論を用いて考察する．高次元漸近論から示唆される，MCMC 法デザインの一般的な注意も述べる．MpCN 法を実際に使う際の注意についても述べたい．本発表は論文 [2, 4] にもとづいている．
本研究は JSPS 科研費 24740062 及び CREST, JST の助成を受けている．

参考文献[1] Beskos, A., Roberts, G., Stuart, A., and Voss, J. MCMC methods for diﬀusion bridges. Stoch.
Dyn. 8, 3 (2008), 319–350.
[2] Kamatani, K. Eﬃcient strategy for the Markov chain Monte Carlo in high-dimension with heavytailed target probability distribution. Arxiv (2014).
[3] Kamatani, K. Rate optimality of Random walk Metropolis algorithm in high-dimension with heavytailed target distribution. Arxiv (2014).
[4] Kamatani, K., and Uchida, M. Hybrid multi-step estimators for stochastic diﬀerential equations
based on sampled data. Statistical Inference for Stochastic Processes 18, 2 (2015), 177–204.
[5] Roberts, G. O., Gelman, A., and Gilks, W. R. Weak convergence and optimal scaling ofrandom walk Metropolis algorithms. Ann. Appl. Probab. 7, 1 (1997), 110–120.
[6] Roberts, G. O., and Rosenthal, J. S. Optimal scaling for various Metropolis-Hastings algorithms.
Statist. Sci. 16, 4 (2001), 351–367.
[7] Tierney, L. Markov chains for exploring posterior distributions. Ann. Statist. 22, 4 (1994), 1701–
1762. With discussion and a rejoinder by the author.

変分鞍点近似によるベイズ推論とモデル選択
国立情報学研究所，JST ERATO 河原林巨大グラフプロジェクト
林浩平

January 21, 2016
高次元データに対してなんらかの構造を仮定し，低次元の特徴を取り出すことはデータ解析においてしばしば行われている．例えば多数の文書から単語のク
ラスタリングを行うトピック抽出は，単語の生成確率がそれぞれの文書で異なる重み付けを持つ (しかし基底となる分布は全文書で共有される) 混合分布に従うという仮定のもとでモデリングし，パラメータ推定を行う．このような例は次元圧
縮，ネットワーク分割，特徴抽出などさまざまな機械学習の問題に頻出する．
このような問題を統一的に扱う枠組みの 1 つとして，隠れ変数モデルを用いたアプローチが存在する．N 個のサンプルからなるデータを X = (x1 , . . . , xN ) とする．隠れ変数モデルは個々のサンプル xn を隠れ変数 zn とパラメータ Θ で表現する．さきほどのトピック抽出の例を用いると xn が文書 n の bag-of-words 表
現，Θ が基底となる単語分布，zn が文書 n が従う単語分布を構成する混合分布の重みとなる．ベイズ統計では隠れ変数モデルは同時確率分布

p(X, Z, Θ | K) = p(X | Z, Θ, K)p(Z | Θ, K)p(Θ | K)

(1)として表現され (K は Z の次元)，隠れ変数やパラメータはベイズ推論によって推定される．
隠れ変数モデルを扱う上で問題となるのが• ベイズ推論の計算困難性• K をどう決めるかの 2 点である．まずベイズ推論に関しては同時確率分布 (1) の積分計算を解く必要があるが，この計算は一般的には解析解を持たないためなんらかの近似計算を行う必要がある．また K の決めかたとして，周辺尤度}∫ {∫p(X | K) =p(X, Z | Θ)p(Θ)dΘ dZ(2)が最大となる K を選択するベイズモデル選択が存在するが，ここの積分計算も一般に計算困難である．
この 2 つの問題点を解決する新しい方法として，因子化漸近ベイズ推論法が近年提案された [1]．因子化漸近ベイズ推論の基本的なアイデアは変分法と鞍点近似を組み合わせた点にある．まず変分法を使うと周辺尤度はlog p(X | K) = Eq [log p(X, Z | K)] + H(q) + KL(qkp∗ )(3)

と変形できる．ここで q(Z) は Z の任意の分布，H(q) はエントロピー，KL(qkp∗ )∗は∫ KL 距離，p (Z) = p(Z|X, K) は周辺事後分布である．さらに p(X, Z | K) =p(X, Z, Θ | K)dΘ に対して鞍点近似を用いることで，ある仮定を満たすとき，
以下のように精度のよい近似が得られる：

log p(X, Z | K) = log p(X, Z | Θ̂, K) + log p(Θ̂)
−

N
1
1
|Θ|
log det |F| −log+ O( ).
2
2
2π
N

(4)

ここで Θ̂ は最尤推定量，F は Θ̂ まわりでの − log p(X, Z | Θ̂, K) のヘッセ行列，|Θ| はパラメータの次元を表すとする．この結果 (4) を (3) に代入することで周
辺尤度の近似が得られる．この近似はもとの周辺尤度と異なり計算困難な積分計
算を持たないか，あるいは持っていたとしてもある程度タイトな下限を構築することができ，反復法による推定アルゴリズムで解くことができる，
因子化漸近ベイズ推論法はもともと混合分布 [1] に対して提案されたが，その後隠れマルコフモデル [2]，混合エキスパート [3, 4]，スパース特徴選択 [5]，関係データモデル [6, 7, 8] などさまざまなモデルに対して拡張された．本発表ではその基本的なアイデアと応用事例について紹介する．

References[1] Ryohei Fujimaki and Satoshi Morinaga. Factorized asymptotic bayesian inference
for mixture modeling. In AISTATS, 2012.
[2] Ryohei Fujimaki and Kohei Hayashi. Factorized asymptotic bayesian hidden
markov model. In International Conference on Machine Learning (ICML), 2012.
[3] Riki Eto, Ryohei Fujimaki, Satoshi Morinaga, and Hiroshi Tamano.
automatic bayesian piecewise sparse linear models. In AISTATS, 2014.

Fully-[4] Jialei Wang, Ryohei Fujimaki, and Yosuke Motohashi. Trading interpretability
for accuracy: Oblique treed sparse additive models. In Proceedings of the 21th
ACM SIGKDD International Conference on Knowledge Discovery and Data Mining (KDD), 2015.
[5] Yohei Kondo, Kohei Hayashi, and Shin ichi Maeda. Bayesian masking: Sparsebayesian estimation with weaker shrinkage bias. In 7th Asian Conference on Machine Learning (ACML), 2015.
[6] Kohei Hayashi and Ryohei Fujimaki. Factorized asymptotic bayesian inference forlatent feature models. In 27th Annual Conference on Neural Information Processing Systems (NIPS), 2013.
[7] Kohei Hayashi, Shin-ishi Maeda, and Ryohei Fujimaki. Rebuilding factorized information criterion: Asymptotically accurate marginal likelihood. In International
Conference on Machine Learning (ICML), 2015.
[8] Chunchen Liu, Lu Feng, Ryohei Fujimaki, and Yusuke Muraoka. Scalable model
selection for large-scale factorial relational models. In International Conference on
Machine Learning (ICML), 2015.

グラフ表現を用いた異常変数同定における一致推定量 ∗†原聡‡

IBM 東京基礎研究所

1 グラフ表現を用いた異常変数同定
D 次元の確率変数を X = (X1 , X2 , . . . , XD ) ∈ RD とする．異常変数同定とは，2 つのデータセット
i.i.d.

i.i.d.

(m) M
P := {x(n) }N}m=1 ∼ q(X) とを比較し，分布に変化のある異常な確率変数の
n=1 ∼ p(X) と Q := {y集合 S ⊆ {1, 2, . . . , D} を探す問題である．異常変数同定に対して複数の先行研究があるか，どれも実験的な性能評価に留まっており理論的な性能保証は行われていなかった [1, 2, 3, 4]．本研究では異常変数同定に対して解の一致性を保証する手法を提案する．なお，問題としては井手ら [1, 2] のグラフ表現を用いた定式化を対象とする．ここで，グラフ表現とは確率変数間の関係性をグラフで表現したものである．例えば，井手らの研究ではグラフ表現として確率変数間の相関や，ガウシアン・グラフィカルモデルが用いられている [1, 2]．
問題 1（グラフ表現を用いた異常変数同定） Λp , Λq ∈ RD×D をそれぞれ分布 p, q のグラフ表現の重みつき隣接行列とする．このとき，以下の (i), (ii) を満たす集合 S ⊆ {1, 2, . . . , D} を見つける．ただし S c は S の補集合である．

(i) Λpdd′ = Λqdd′ , ∀d, d′ ∈ S c , (ii) Λpdd′ 6= Λqdd′ , ∀d ∈ S, ∃d′ ∈ {1, 2, . . . , D} \ {d}.
条件 (i) は異常変数集合 S の補集合，つまり正常変数においてグラフ表現に変化がないことを意味する．他
方，条件 (ii) は異常変数においては隣接する少なくとも 1 つの辺においてその重み，つまり確率変数間の関係に変化があることを意味する．これらは近傍構造保存仮定として井手ら [1, 2] により提唱された条件である．

2 k-疎部分グラフ問題による定式化，解法及び一致性
データセット P, Q から推定されたグラフ表現をそれぞれ Λ̂p , Λ̂q とする．また，行列 Γ̂ ∈ RD×D をΓ̂dd′ := |Λ̂pdd′ − Λ̂qdd′ | により定義する．このとき，推定量 Λ̂p , Λ̂q が真の隣接行列 Λp , Λq に十分近ければ，条p

件 (i), (ii) から以下が言える：(i’) Γ̂dd′ ≈ 0, ∀d, d′ ∈ S c , (ii’) Γ̂dd′ > 0, ∀d ∈ S, ∃d′ ∈ {1, 2, . . . , D} \ {d}．

(i’), (ii’) は集合 S c （及びその補集合 S ）が，行列 Γ̂ の値の小さい部分行列を探すことで同定できることを示している．異常変数集合 S のサイズが k だとすると，これは以下の最適化問題として定式化できる．

Ŝ =argminXΓ̂dd′ , s.t. |S| = k.

(1)

S⊆{1,2,...,D} d,d′ ∈S c

∗ 本稿の執筆にあたり JST CREST の助成を受けた．
† 本研究は IBM 東京基礎研究所の森村哲郎氏，高橋俊博氏，柳沢弘輝氏，及び東工大/さきがけの鈴木大慈氏との共同研究である．
‡ satohara@jp.ibm.com

問題 (1) は Γ̂ を隣接行列とするグラフの k-疎部分グラフ問題 [5] であり，一般には NP 困難である．k の値が既知であれば問題 (1) の解法としては，厳密解法や貪欲法による近似解の導出などが考えられる．特に，厳密解については以下の一致性定理が成り立つ．
定理 1（Ŝ の一致性） 推定量 Λ̂p , Λ̂q が真の隣接行列 Λp , Λq の一致推定量ならば，問題 (1) の解 Ŝ は集合 Sの一致推定量である．つまり，limN,M→∞ P (Ŝ =6 S) = 0 が成り立つ.
集合 S のサイズ k が未知の場合は，問題 (1) の凸近似から導かれる凸二次計画を解くことで S の推定量を得ることができる．

S̃ = {d | s̃d = 0}, s̃ = argmin s⊤ Aµ s, s.t.
s∈RD

D
X

sd = 1, sd ≥ 0, ∀d.

d=1

ただし，行列 Aµ ∈ RD×D は Aµ := Γ̂ + µID （ID は D 次元単位行列）であり，定数 µ は Aµ が正定値になる十分大きな実数とする．このとき，推定量 S̃ について以下の一致性定理が成り立つ．
定理 2（S̃ の一致性） 推定量 Λ̂p , Λ̂q が真の隣接行列 Λp , Λq の一致推定量であり，かつ以下を満たす δ > 0が存在するならば S̃ は集合 S の一致推定量である．

min
d∈S

X

|Λpdd′ − Λqdd′ | ≥ (1 + δ)µ.

d′ ∈S c定理 2 は k の値が未知でも一定条件下で集合 S の一致推定量が得られることを示している．また，S̃ は凸二次計画を解くことで得られるため，必要な計算時間を多項式時間で抑えることができる．これは NP 困難な問題 (1) を直接解くよりも計算時間の面で有利である．

参考文献[1] T. Idé, S. Papadimitriou, and M. Vlachos. Computing correlation anomaly scores using stochastic
nearest neighbors. Proceedings of the 7th IEEE International Conference on Data Mining, pages
523–528, 2007.
[2] T. Idé, A. C. Lozano, N. Abe, and Y. Liu. Proximity-based anomaly detection using sparse structurelearning. Proceedings of the 2009 SIAM International Conference on Data Mining, pages 97–108,
2009.
[3] S. Hirose, K. Yamanishi, T. Nakata, and R. Fujimaki. Network anomaly detection based on eigen
equation compression. Proceedings of the 15th ACM SIGKDD International Conference on Knowledge
Discovery and Data Mining, pages 1185–1194, 2009.
[4] R. Jiang, H. Fei, and J. Huan. Anomaly localization for network data streams with graph joint sparse
PCA. Proceedings of the 17th ACM SIGKDD International Conference on Knowledge Discovery andData Mining, pages 886–894, 2011.
[5] R. Watrigant, M. Bougeret, and R. Giroudeau. Approximating the sparsest k-subgraph in chordal
graphs. Approximation and Online Algorithms, Lecture Notes in Computer Science, 8447:73–84,
2014.

Non-convex Compressed Sensing
with the Sum-of-Squares Method
相馬 輔 (東京大学)
𠮷田 悠一 (国立情報学研究所・Preferred Infrastructure)

1 概要
本研究では，ℓq ノルム (0 < q ≤ 1) に対する stable signal recovery を考える．この問題では，未知の信号 x ∈ Rn に対して，行列 A ∈ Rm×n を用いたベクトル y = Ax が観測できるものとする．目標は，観測
行列 A と 復元アルゴリズム ∆ : Rm → Rn であって，∥x − ∆(Ax)∥q ≤ O(1) · min{∥x − x∗ ∥q : ∥x∗ ∥0 ≤ s}となるものを設計することである．既存の ℓq -stable signal recovery に対するアルゴリズム [1] は ℓq 最小化

min ∥z∥qq

subject to Az = y(1)を解く必要があるが，これは非凸最適化問題であるため，多項式時間の ℓq -stable signal recovery アルゴリ
ズムの存在は分かっていなかった．本研究では，二乗和多項式法 (Sum-of-Squares Method) を利用し，Rademacher 行列では，定数 0 < q < 1 に対して (ほぼ)ℓq -stable signal recovery が多項式時間で可能である
ことを示す．
定理 1. q = 2−k (k ∈ Z+ ), s, n ∈ Z+ とする．A ∈ Rm×n を m = O(s2/q log n) を満たす正規化

Rademacher 行列とすると，以下が高確率で成り立つ: 任意の x ∈ [−1, 1]n と ϵ > 0 に対して
∥x − ∆(Ax)∥q ≤ O(1) · min{∥x − x∗ ∥q : ∥x∗ ∥0 ≤ s} + ϵを満たす多項式時間アルゴリズム ∆ が存在する．

表 1 ℓq -stable signal recovery の既存手法と提案手法[2, 3]
[1]
本研究

q
q=1
0<q≤1
0 < q ≤ 1(定数)観測行列 A

サンプル数 mアルゴリズム ∆∆ の計算量Subgaussian
Subgaussian
RademacherO(s log(n/s))
C1 (q)s log(n/s) + C2 (q)sO(s2/q log n)ℓ1 最小化 (LP)ℓq 最小化

多項式時間

二乗和多項式

多項式時間

2 多項式最適化と二乗和多項式法
多項式最適化は，多項式 f, g1 , . . . , gi ∈ R[z] (z = [z(1), . . . , z(n)]⊤ ) に対して

min
z

f (z)

sub. to gi (z) = 0
1(i = 1, . . . , m)—で定義される最適化問題である．多項式最適化は一般に非凸最適化問題であり，整数計画問題等を含むためNP 困難である．
多項式最適化に対する強力な凸緩和法として二乗和多項式法 [4, 5, 6, 7] がある．二乗和多項式法においては，次数と呼ばれるパラメータ d を指定する．そして，元の変数 z のかわりに pseudoexpectation と呼ば

e : R[z]d → R を変数とした最適化問題を解く (R[z]d は次数 d 以下の多項式の集合)．
れる線形作用素 Emine (z)]
E[fsub. to

e =1
E[1]

e
E

2
e
E[p(z)]≥0(p ∈ R[z] : deg(p) ≤ d/2)

e i (z)p(z)] = 0 (p ∈ R[z] : deg(gi p) ≤ d, i = 1, . . . , m)E[g二乗和多項式法は nO(d) サイズの行列変数を持つ半正定値計画 (SDP) に帰着可能であることが知られている．特に，d が定数であれば多項式時間で (任意の精度で) 解くことが可能である．

3 提案アルゴリズムℓq 最小化は多項式最適化として定式化可能であるが，愚直な定式化のままでは緩和ギャップが大きすぎるため，ℓq 最小化の最適解の情報を得ることができない．そこで提案アルゴリズムでは，冗長な変数と制約を追加した多項式最適化に対して二乗和多項式緩和を適用することで，緩和ギャップを減らし ℓq 最小化の最適解の情報を得られるようにした．追加する制約は，ℓq 最小化の stability の証明に必要な ℓqq 擬距離における三角e から元信号の推定量であるベクトルを抽出するアルゴリズムも構成した．
不等式などである．また，緩和解 E
詳細は [8] を参照されたい．

参考文献[1] R. Chartrand and V. Staneva, “Restricted isometry properties and nonconvex compressive sensing,”
Inverse Problems, vol. 24, no. 3, 2008.
[2] E. J. Candès, “The restricted isometry property and its implications for compressed sensing,” Comptes
Rendus Mathematique, vol. 346, no. 9-10, pp. 589–592, 2008.
[3] E. J. Candès, J. K. Romberg, and T. Tao, “Stable signal recovery from incomplete and inaccurate
measurements,” Communications on Pure and Applied Mathematics, vol. 59, no. 8, pp. 1207–1223,
2006.
[4] J. B. Lasserre, “Global optimization with polynomials and the problem of moments,” SIAM Journal
on Optimization, vol. 11, no. 3, pp. 796–817, 2006.
[5] Y. Nesterov, “Squared functional systems and optimization problems,” in High Performance Optimization. Springer, 2000, pp. 405–440.
[6] P. A. Parrilo, “Structured semidefinite programs and semialgebraic geometry methods in robustness
and optimization,” Ph.D. dissertation, California Institute of Technology, 2000.
[7] N. Z. Shor, “An approach to obtaining global extremums in polynomial mathematical programming
problems,” Cybernetics, vol. 23, no. 5, pp. 695–700, 1987.
[8] T. Soma and Y. Yoshida, “Non-convex compressed sensing with the sum-of-squares method,” inProceedings of the 26th Annual ACM-SIAM Symposium on Discrete Algorithms, 2016, pp. 570–579.

2教師付き学習における Barron and Cover 理論
BARRON AND COVERS’ THEORY IN SUPERVISED LEARNING
川喜田雅則
Masanori Kawakita
九州大学システム情報科学研究院
Information Science and Electrical Engineering, Kyushu University

1

まえがき

機械学習，統計学，情報理論において，様々な推定法の精度を理論的に評価，あるいは保証することは重要であるが，これまで従来の理論の多くは標本数，あるいは共変量の数を無限大と仮定したり，確率変数の有界性，
あるいはモーメントに関する様々な技術的な条件を必要とする．これらの条件は現実に推定を行うときには成り立っていないものが少なくない．そこで，本稿ではこのような仮定をほとんど必要とせず性能評価を行う方法を考えてみたい．そのための土台となる理論が MDL(最小記述長) 原理 [1, 2, 3, 4] に関する代表的な成果の一つである Barron and Cover 理論 (BC 理論)[5] である．MDL
原理とはデータを可能な限り短く記述することが，データを生成した確率分布に関する推定の精度を高めることにつながるという原理である．MDL 原理に基づくモデ
ル選択規準として最も有名なものの一つが Rissanen によって提案された二段階符号に基づく情報量基準 [1] である．ここでいう二段階符号とは，まずデータを記述するのに用いる確率モデルを記述し，次にそのモデルを用
いてデータを記述するという符号化を指し，このモデル
記述長とデータ記述長の和を最小にする推定量を MDL
推定量と呼ぶ．BC 理論は MDL 推定量に対して，次のような理想的な性質を満たすリスクの上界評価を与える理論である．まず BC 理論は確率変数の有界性やモーメント条件などの技術的な条件をほとんど必要としない．
また有限の標本数で成立し，（経験的に）タイトなリスクの上界評価を与える．このような理論は著者の知る限り他にない．しかし BC 理論を機械学習に教師付き学習に応用しようとすると，これまで近似的かつ限定的にしか可能ではなかった．この問題の解決には BC 理論にとって本質的な困難さが伴う．最近我々はこの問題を一定程度解決し，BC 理論を教師付き学習に適用するための一つの方法を提案し，実際に lasso の乱数計画におけるリスク上界を導出することに成功した [6, 7, 8, 9] ので，それを本稿で紹介する．

2

教師付き学習における MDL 推定量

あ る 分 布 p∗ (x, y)=p∗ (x)p∗ (y|x) か ら 独立に生成された n 個の教師付きデータ D:={(x1 , y1 ), (x2 , y2 ), · · · , (xn , yn )} が 与 え ら れ た と する．以降では xn := x1 x2 · · · xn ，y n := y1 y2 · · · yn と書く．教師付き学習とは，このデータセット D を用いて，p∗ (y|x) を推定する問題である．ここでは p∗ (y|x)
を推定するためにパラメトリックモデル pθ (y|x) を用いる．ただし θ ∈ Θ(= ℜp ) は p 次元のパラメータベクト
ルとする．MDL 推定量を定義するためにはパラメータ
空間 Θ を量子化しなくてはならない．量子化近似したパラメータ空間を Θ(xn )(⊂ Θ) と書く．また Θ(xn ) に含まれる各モデルの記述長を L(θ̄|xn ) と定義する．モ
デル記述長は Kraft の不等式を満たすとする．モデルpθ (y|x) に関する MDL 推定量 θ̈(xn , y n ) は{}θ̈(xn , y n ) = arg min − log pθ̄ (y n |xn ) + βL(θ̄|xn )θ̄∈Θ(xn )と定義される．ただし，β は 1 より大きい実数とする．この記述長は，共変量 xn を与えられた下での，目的変数 y nの条件付き記述長であることに注意する．最小化の目的
関数の第 1 項をデータ記述長，第 2 項をモデル記述長と呼ぶ．また，それらの和を全記述長と呼ぶ．なお，本稿では log で自然対数を表し，符号長の単位に nat を用いる．
この二段階符号化によって達成できる記述長の最小値をLβ2-p (y n |xn ) で表す．このとき，Lβ2-p (y n |xn ) は y n について Kraft の不等式を満たすため，語頭符号の符号長とみなせる．従って，pβ2-p (y n |xn ) := exp(−Lβ2-p (y n |xn )) =pθ̈ (y n |xn ) exp(−βL(θ̈|xn )) とすれば，pβ2-p は条件付きの劣確率分布となる．
教師付き学習における BC 理論
オリジナルの BC 理論 [5] は，MDL 推定量についての理論であり，推定量は量子化近似された空間上で定義されていなければならなかった．それに対して本節では量子化近似をしない推定量3

θ̂(xn , y n ) := arg min {− log pθ (y n |xn ) + Pen(θ|xn )}θ∈Θを考える．BC 理論を教師付き学習について，このように量子化近似をしない推定量に適用することは本質的に難しかったが，[6, 10, 7, 8, 9] で提案された解決法について紹介する．まず BC 理論では推定量の精度を測るための損失関数に Rényi ダイバージェンス( q(Y |X) )1−λ
1log Ep∗ (x)p(y|x)
dλ (p, q) = −1−λp(Y |X)を用いる．ここで λ ∈ (0, 1) は Rényi ダイバージェンスの次数と呼ばれる．また真の分布 p∗ (x) についての適当
(n)
(n)
(n)
な典型集合 Aϵ を導入し，Pϵ = Pr{xn ∈ Aϵ } と書
く．この典型集合を用いて次の定義を導入する．
Definition 1. ある β > 1 について，各 p∗ (xn ) ∈ PX
に対して，量子化空間 Θ(p∗ ) ⊂ Θ と，Kraft の不等式を満たすモデル記述長 L∗n : Θ(p∗ ) → [0, ∞) が存在し，それらが共変量 xn に依存せず，n∀xn ∈ A(n)ϵ , ∀y ,)(p∗ (y n |xn )
∗
−
βL
(
θ̄)
max ndλ (p∗ , pθ̄ ) − lognpθ̄ (y n |xn )θ̄∈Θ(p∗X)()
p∗ (y n |xn )n≥ max ndλ (p∗ , pθ ) − log−
Pen(θ|x)θ∈Θpθ (y n |xn )(n)が成り立つとき，Pen は (PX , β, Aϵ ) に対して ϵ リス
ク妥当であるという．

ϵ リスク妥当である罰則関数について，次の定理を示せる．
(n)定理 2. Pen(·|xn ) が (PX , β, Aϵ ) について ϵ リスク
妥当であるとき，任意の λ ∈ (0, 1 − β −1 ] について

E1 dλ (p∗ , pθ̂ ) ≤ E1

1p∗ (y |x )β
1log
+ log (n)nnnpβ2-p (y |x ) nPϵnnが成り立つ．ただし，pβ2-p の定義の中で用いられる推(n)
定量は θ̂ とし，E1 は xn ∈ Aϵ の下での条件付き期待値を表わす．
証明は例えば [9] を参照されたい．リスク上界の第二項は典型集合の性質から n に応じて小さくなる項である．従ってリスク上界の主要な項は第一項であり，これはオリジナルの BC 理論と同様に二段階符号の（条件付き）冗長度の形をしている．従って冗長度が小さくなるように二段階符号を設計することがリスクの良い保証（上界）を得ることにつながるという意味で，「BC 理論による MDL 原理の正当化」を教師付き学習の場合に拡張していることになる．ただし主要項が厳密に冗長度とみなせるためには，一般には罰則関数に別途条件が必要であることに注意する．

lasso のリスク上界
本節では lasso について，前節の ϵ リスク妥当性を用いて乱数計画における lasso の具体的なリスク上界を導出する．まず教師付きデータ {(xi , yi )|i = 1, 2, · · · , n} は回帰モデル yi = xTi θ∗ + ϵi に従うとする．パラメータ θの次元を p とし，n より大きくても良いとする．ノイズ
ϵi は平均 0, 分散 σ 2 を持つ任意の確率変数とする．また
PX := {p(xn ) = Πni=1 N (xi ; 0, Σ)| 任意の正則なΣ} とする．ただし N (x; µ, Σ) は平均ベクトル µ，分散共分散行
列 Σ の正規分布を表すとする．正規性を仮定しても ϵ リ
スク妥当な罰則関数を求めるのは容易ではないが，著者
らが文献 [6, 9] において導出した ϵ リスク妥当な重み付(n)
き ℓ1 ノルムを紹介する．まず典型集合 Aϵ には∑n}{(1/n) i=1 x2ij
(n)
n
≤1+ϵ
Aϵ := x ∀j, 1 − ϵ ≤Σjj4を採用する．ただし ϵ ∈ (0, 1) とする．このとき以下の定理が成り立つ．

定理 3. 次数 λ の Rényi ダイバージェンスに対して，(n)次の条件を満たす重み付き ℓ1 ノルムは，(PX , β, Aϵ )についての ϵ リスク妥当性を満たす罰則関数となる．

Pen(θ|xn ) :=
µ1

≥

µ2

≥

µ1 ∥θ∥w,1 + µ2 ,
v
(
)
u
√
u 2n log 4p
λ
1+ϵ
t
+
,
σ2 β
4(1 − ϵ)
1−ϵ
log 2
.
β

lasso 推定量は定理 2 に示されたリスク上界を持つ．
(n)
また，文献
2 に示されているように Pϵ ≥)( [6] の補題
2
1 − 2p exp − nϵ7 が成り立つ．従って定理 2 のリスク上界の第二項は指数的に小さくなる項であることがわかる．
参考文献[1] J. Rissanen, “Modeling by shortest data description,” Automatica, vol. 14, no. 5, pp. 465–471,
1978.

[2] A. R. Barron, J. Rissanen, and B. Yu, “Theminimum description length principle in coding
and modeling,” IEEE Transactions on Information Theory, vol. 44, no. 6, pp. 2743–2760, 1998.
[3] P. Grünwald, The Minimum Description Length
Principle. MIT Press, 2005.
[4] P. Grünwald, I. Myung, and M. Pitt, Advances in
Minimum Description Length: Theory and Applications, ser. Neural information processing series.
Bradford Book, 2005.
[5] A. R. Barron and T. M. Cover, “Minimum complexity density estimation,” IEEE Transactions
on Information Theory, vol. 37, no. 4, pp. 1034–
1054, 1991.
[6] 川喜田雅則, 豊暉原侑心, 竹内純一, “MDL 理論による lasso のリスク上界,” 信学技法 IBISML2015-16,
vol. 115, no. 112, pp. 101–107, 2015.
[7] 竹内純一, 川喜田雅則, “教師付き学習における MDL
推定量の設計とその収束速度,” 信学技法 IT2015-26,
vol. 115, no. 137, pp. 53–58, 2015.
[8] ——, “教師付き学習における MDL 推定,” 第 9 回
シャノン理論ワークショップ予稿集, 2015.
[9] 川喜田雅則, 竹内純一, “Barron and cover 理論の教
師付き学習への拡張及び lasso への応用.” 第 18 回
情報論的学習理論ワークショップ（IBIS 2015), 11
2015.
[10] ——, “教師付き学習における MDL 推定のリスク
バウンドについて,” 第 38 回情報理論とその応用シンポジウム (SITA 2015), pp. 78–83, 11 2015.

状態空間モデルを用いた時系列の周期成分分解と位相推定
東京大・情報理工松田 孟留

東京大・情報理工，理研・脳センター

駒木 文保

多くの時系列は，さまざまな周波数の周期成分が重ね合わさったものと見なすことができる．たとえば脳波や Local Field Potential などの神経電位時系列は，その周期成分の位相が生理学的に重要な意味をもつことが近年明らかにされてきている [2]．従来はバンドパスフィルタリングとヒルベルト変換によって位相を推定している．すなわち，バンドパスフィルタをかけた時系列 y(t) を∫ ∞yH (t) =

y(τ )
dτ,
t
−∞ − τ

z(t) = y(t) + iyH (t)によって解析信号 z(t) に変換し，その偏角を位相としている．しかし，このような推定法は観測ノイズやフィ
ルタリング方法によって結果が不安定になってしまう [4]．
本研究では，状態空間モデルに基づいた時系列の周期成分分解と位相推定の手法を開発した．Wiener は，脳波のスペクトルにおいてアルファ帯域（10Hz 付近）にピークができる原理について，確率微分方程式のモデルを用いて考察した [5]．このモデルは，振動子の周波数がランダムにゆらぐ（「時計の針が震える」）という発想から導入されたものである．提案手法では，このランダム周波数変調の発想をもとに線形ガウス型の状態
空間モデル

(

(k)

xt+1,1(k)
xt+1,2)((cos(2πfk ∆t)∼ N aksin(2πfk ∆t)− sin(2πfk ∆t)cos(2πfk ∆t)yt ∼ N(K∑) ( (k) ) ( 2
xt,1
σk
(k) ,
0
xt,2

0
σk2

))
(k = 1, · · · , K),)(k)xt,1 , τ 2

k=1を用いる．ここで ∆t は yt のサンプリング間隔である．このモデルでは，時系列 yt には K 個の（ランダム周波数変調あり）振動子の活動が重ね合わさっていると仮定し，k 番目の振動子の周波数を fk ，座標を(xt,1 , xt,2 )⊤ としている．よって，yt にカルマンフィルタ・スムーザを適用することで，各振動子の振幅と位
(k)

(k)

相が推定できる．これは，時系列をトレンド成分や季節成分に分解する [3] の季節調整法と同じ発想である．
2
状態空間モデルのパラメータ a1 , · · · , aK , f1 , · · · , fK , σ12 , · · · , σK
, τ 2 は経験ベイズ法によって，また振動子の個数 K は情報量規準 ABIC[1] の最小化によってあらかじめ決定しておく．これにより，時系列を自然な形で（ランダム周波数変調ありの）周期成分に分解することができる．

参考文献[1] H. Akaike. Likelihood and the Bayes procedure. In Bayesian Statistics, Eds. J. M. Bernardo,
M. H. DeGroot, D. V. Lindley and A. F. M. Smith, pp. 1–13. Valencia, Spain, 1980.
[2] G. Buzsáki. Rhythms of the Brain. Oxford University Press, Oxford, 2011.
[3] G. Kitagawa and W. Gersch. A smoothness priors-state space modeling of time series with trend
and seasonality. Journal of American Statistical Association, Vol. 79, 378–389, 1984.
[4] A. G. Siapas, E. V. Lubenov and M. A. Wilson. Prefrontal phase locking to hippocampal theta
oscillations. Neuron, Vol. 46, pp. 141–151, 2005.
[5] N. Wiener. Nonlinear problems in random theory. MIT Press, Cambridge, 1966.

5

0

10

20

30

40

50

60

70

80

90

100

10

20

30

40

50

60

70

80

90

100

10

20

30

40

50

60

70

80

90

100

10

20

30

40

50

60

70

80

90

100

20

30

40

50

60

70

80

90

100

5

0

5
0
-5

0.5
0
-0.5
-1

5

x 10

-8

0
-5

図1

10

提案手法による Wolfer sunspot data の分解結果．1 行目：元データ（対数変換した黒点時系列），2〜4 行目：周期成分，5 行目：観測ノイズ

関数予測における漸近ミニマックス予測分布
東京大・情報理工 矢野 恵佑
東京大・情報理工, 理研・脳科学総合研究センター

駒木 文保

１ 概要
本発表では，ノンパラメトリックな関数予測における漸近ミニマックス予測分布の構成を議論する．まず，真の関数が楕円型の制約をもち我々がその制約値を知っている場
合の漸近ミニマックス予測分布を構成する．次に，真の関数が楕円型の制約をもち我々がその制約値を知らない場合の漸近ミニマックス予測分布を構成する．最後に，数値実験を通して我々の構成した漸近ミニマックス予測分布の性質を議論する．

２ 関数予測における漸近ミニマックス予測分布
未知の関数 θ(·) にガウス過程 W (·) が加法された曲線 X(·) を観測し，θ(·) に W (·) とf (·) が加法された曲線 Y (·) を予測する状況を考える．X(·) は既知の独立なガウス過程 W
ε > 0 および [0, 1] 上のガウス過程 W (·) を利用してX(t) = θ(t) + εW (t) for t ∈ [0, 1](1)f (·) を利用してと与えられ，Y (·) は既知の ε̃ > 0 および [0, 1] 上のガウス過程 W
f (t) for t ∈ [0, 1]
Y (t) = θ(t) + ε̃W(2)と与えられるとする．
我々は上記の予測を次のようなガウス列モデルでの予測に帰着して考える．すなわち，
f を R∞ 上のガウス分布 ⊗∞ N (0, 1) に従う独立な確率変数とし，ε と ε̃ を既知
W とW
i=1
f に従う Y の分布をの正数とした時に，X = θ + εW に従う X を観測し，Y = θ + ε̃W予測する．X の従う分布を Pθ と書き，Y の従う分布を Qθ と書く．
本発表では，パラメータ空間が楕円型である場合の Kullback–Leibler リスクのもとでの漸近ミニマックスな予測分布を構成する．予測分布とは観測 X をもとに構成した R∞
上の確率分布のことである．楕円型のパラメータ空間は非零な増加列 a = (a1 , a2 , · · · ) お
よび正数 C を利用して∞

Θ(a, C) = {θ ∈ l2 : Σ a2i θi2 ≤ C}

(3)

i=1

と与えられる．Kullback–Leibler リスクのもとで漸近ミニマックスな予測分布とは∫dQθ(y)dQθ (y)dPθ (x)
inf
sup
log
dQ̂(·; X = x)Q̂(·;X)∈A θ∈Θ(a,C)∫
dQθ
= (1 + o(1)) suplog(y)dQθ (y)dPθ (x) as ε → 0
dQ̂minimax (·; X = x)θ∈Θ(a,C)(4)を満たす予測分布 Q̂minimax (·; X) である．ここで，A は R∞ 上の確率分布全体であり，ε̃は 1/ε̃ = O(1/ε) と与えられているとする．
まず，増加列 a および C が既知である場合にある τ ∗ (ε, a, C) ∈ l2 によるガウス分布∗2Gτ ∗ (ε,a,C) = ⊗∞i=1 N (0, τi (ε, a, C)) を事前分布に用いたベイズ予測分布∫
QGτ ∗ (ε,a,C) (·|X = x) :=Qθ (·)dGτ ∗ (ε,a,C) (θ|X = x)

(5)が漸近ミニマックスであることを示す．ここで，Gτ ∗ (ε,a,C) (·|X = x) は Gτ ∗ (ε,a,C) の事後分布である．[1] はパラメータの次元数が ε に応じて増加するという漸近論を利用して楕
円型のパラメータ制約下での Kullback–Leibler リスクのもとでの漸近ミニマックス予測分布を求めている．本結果は [1] の結果の無限次元のパラメータ空間への拡張である．
次に，増加列 a および C が未知の場合での漸近ミニマックスを達成する予測分布の構成法について紹介する．我々の構成法は有限次元の Stein の事前分布を利用する．d 次元∑
ユークリッド空間における Stein の事前分布とは密度が h(θ(d) ) := ( di=1 θi2 )(2−d)/2 と表される事前分布である．
最後に，数値実験を通して有限の ε に対する漸近ミニマックス予測分布の性質を議論
する．

参考文献
[1] Xu, X. and Liang, F. (2010), Asymptotic minimax risk of predictive density estimation for non-parametric regression. Bernoulli 16, pp. 543–560.

低次元多様体上の Data Sharpening
島根大学総合理工学研究科

工藤 雅紀

島根大学総合理工学研究科

内藤 貫太

1．はじめに ノンパラメトリック核型回帰において，バイアス縮小を導く方法として Data Sharpening が知られている．Data Sharpening という手法は [2] で提案され，局所線形 (以下，LL) 推定量に残差の LL 推定量を加えることで Data Sharpening(以下，DS) 推定量が得られる．そして，p 変量の説明変数に対する DS 推定量の漸近的振る舞いは，バンド幅 h → 0 と標本数 n → ∞, nhp → ∞ のもとで，漸近バイアスが Op (h4 ) であり，漸近分散は Op (1/(nhp )) である．他方，LL 推定量の漸近的振る舞いは [4] で議論され，その漸近バイ
アスは Op (h2 ) であり，その漸近分散は Op (1/(nhp )) であることが知られている．したがって，LL 推定量よりも，DS 推定量の方が漸近的には小さいバイアスを持つことがわかる．
一方，説明変数がある低次元の多様体に埋め込まれているという仮定を含む回帰モデルが議論されてきた．
そこでのモデルは 2 種類あり，説明変数が埋め込まれている多様体を既知とするモデルと，未知とするモデルに大別される．多様体を既知とするモデルの代表例は方向統計学であり，多様体として円周や球面を想定しており，これまで盛んに研究されている．しかし，多様体が未知であるとする研究は，既知である場合の研究に比べるとあまり多くないように思われ，例えば [1] では，LL 推定量を未知の d 次元多様体上に構成し漸近的振る舞いが調べられ，漸近バイアスが Op (h2 ) であり，その漸近分散は Op (1/(nhd )) であることを導いている．
本発表では，[1] と [2] の結果を伴わせて，説明変数が未知の低次元多様体に埋め込まれている設定で DS 推
定量を構成し，その理論的・数値的検証について報告する．

2. Data Sharpening 推定量 p 変量の説明変数を X ，1 変量の目的変数を Y とし，その独立な n 個の観測値を {(Xi , Yi )|1 ≤ i ≤ n} とする．ε1 , · · · , εn を i.i.d の誤差とし，回帰モデル

Yi = m(Xi ) + v(Xi )1/2 εi , E[εi ] = 0, V ar[εi ] = 1, i = 1, · · · , nを想定する．推定のターゲットは m である．[1] での設定を踏まえ，説明変数が d 次元多様体 X (d ≤ p, X ⊂Rp ) に埋め込まれているとする．すなわち，Xi ∈ X を仮定する．このとき，x ∈ X に対して，Local chart√pqd→ X ∩ Bx,νが存在する．ここで，By,r= {z ∈ Rq | (z − y)T (z − y) < r} は中と言われる写像 φ : B0,µ
心が y で半径が r である q 次元の球．さらに非退化な d 次元の密度関数 f (z) により，f (φ−1 (x)) を X の密度関数とする．Rp ∋ x における m(x) の初期推定量 mb LL (x, h) は，局所線形推定量で構成する．ここで，mb LL (x, h) = eT1 (XxT Wx Xx )T XxT Wx Y であり，
1 ..
Xx =  .

 (X1 − x)TY1 .. 
..
 , Y =  .  , Wx = diag(Kh (X1 − x), · · · , Kh (Xn − x)),.

1 (Xn − x)TYn

K はカーネル関数，h > 0 はバンド幅，Kh (U ) = K (h−p U ) , e1 = [1, 0, · · · , 0]T である．
残差 ri := Yi − mb LL (Xi , h) から，{(Xi , ri )|1 ≤ i ≤ n} による局所線形推定量として rbLL (x, h) =eT1 (XxT Wx Xx )T XxT Wx [r1 , · · · , rn ]T を構成する．Data Sharpening 推定量はf)m
b DS (x, h) := mb LL (x, h) + rbLL (x, h) = eT1 (XxT Wx Xx )T XxT Wx (2Y − M
f = [mとして得られる．ここで，M
b LL (X1 , h) · · · mb LL (Xn , h)]T である．
3. 記号 関数 g : Rq → R (q = d or q = p) に対して, ∇2 g(x) はヘッセ行列を表す．φ(z) = [φ1 (z) · · · φp (z)]T1

と z = (z1 , · · · , zd ) により，a における φ のヤコビ行列を[J (a|φ) =]∂φi (z)∂zj1≤i≤p,1≤j≤dz=aとし，また G : Rp → R に対して，−1∫

B(x|g, φ) = C(x|φ)
uT J (φ−1 (x)|φ)T ∇2 g(x)J (φ−1 (x)|φ)uK(J (φ−1 (x)|φ)u)du,Rd∫
R(x|G, φ) =G(J (φ−1 (x)|φ)u)2 du
Rd

とおく．ここで∫C(x|φ) =K(J (φ−1 (x)|φ)u)duRd

であり，さらに，G∗ (y) = 2G(y) − C(x|φ)−1 G ∗ G(y), G ∗ G(y) =∫RdG(J (φ−1 (x)|φ)v)G(y − J (φ−1 (x)|φ)v))dv.

以上の記号を用いると，[1] の Theorem 2.1 における mb LL (x, h) のバイアスと分散，それぞれのリーディング
タームは J1 (x) = B(x|m, φ) と J2 (x) = v(x)R(x|K, φ)/{f (φ−1 (x))C(x|φ)2 } である．

4. 漸近挙動 適当な条件のもと，x ∈ X , n → ∞, h → 0 のとき
h4
B(x|J1 , φ) + op (h4 ),4
1 R(x|K ∗ , φ)J2 (x)(1 + op (1)).
V ar[mb DS (x, h)|X1 , · · · , Xn ] =nhd R(x|K, φ)

Bias[mb DS (x, h)|X1 , · · · , Xn ] = −[1] で議論された LL 推定量は漸近バイアスが Op (h2 )，漸近分散は Op (1/(nhd )) であった．一方，DS 推定
量の漸近バイアスは Op (h4 )，漸近分散は Op (1/(nhd )) であり，バイアス縮小が生じていることがわかる．

5. DS 推定量の漸近正規性 3 節の条件のもと，E[ε4 ] < ∞ かつ x ∈ X で h = κn−1/(d+8) のとき
(
)/
√
h4
D
V ar[mb DS (x, h)|X1 , · · · , Xn ] −→ N (0, 1)m
b DS (x, h) − m(x) − B(x|J1 , φ)4Dが成り立つ．ここで，κ はある正の定数で， −→ は分布収束を表し，N (0, 1) は標準正規分布である．

6. バイアス縮小の検証 mb DS (x, h) のバイアス縮小が実際に生じていることを，シミュレーションと実データへの適用を通して確認する．結果は本発表にて報告する．

参考文献[1] Bickel, P. and Li, B. (2007). Local polynomial regression on unknown manifolds, Complex Datasets andInverse Problems: Tomography, Networks and Beyond. Institute of Mathematical Statistics Lecture
Notes-Monograph Series, 54, 177-186.
[2] Choi, E., Hall, P. and Rousson, V. (2000). Data sharpening methods for bias reduction in nonparametric regression, The Annals of Statistics, 28, 1339-1355.
[3] Naito, K. and Yoshizaki, M. (2009). Bandwidth selection for a data sharpening estimator in nonparametric regression, Journal of Multivariate Analysis, 100, 1465-1486.
[4] Ruppert, D. and Wand, M. P. (1994). Multivariate locally weighted least squares regression, The
Annals of Statistics, 22, 1346-1370.
[5] Yoshizaki, M. and Naito, K. (2004). Practical aspects of bias reducing estimators in nonparametric
regression, Japanese Journal of Applied Statistics, 33, 131-155 (in Japanese).

2確率的交互方向乗数法とマルチクラスグラフ型正則化学習への応用
鈴木大慈 † , 吉田悠一 ‡†‡東京工業大学，JST さきがけ (s-taiji@is.titech.ac.jp)

国立情報学研究所，Preferred Infrastructure, Inc., (yyoshida@nii.ac.jp)

1 はじめに
本発表では，「多クラス連結正則化」と呼ばれる新しい正則化を導入し，その効率的確率的最適化手法を提案する．
グラフに沿った正則化は半教師あり学習や画像のデノイジング，セグメンテーション，遺伝子データ解析など，多くの場面で現れる．本研究では，そのような正則化の中でもスパース性を有する正則化に焦点を当てる．既存の手法はグラフの各頂点に当てられたラベルが一つであるか，複数であってもそれらが独立に扱われていたり，スパース性を有していないかった．本研究では，k-劣モジュラ関数に発想を得て，新しい正則化を提案する．これまでの研究で一般化 Fused Lasso 正則化 (Tibshirani et al., 2005) はグラフカットの Lovász 拡張となることが知られていた (Bach, 2013)．
グラフカットは劣モジュラ関数であるが，劣モジュラ関数は {0, 1}N 変数を入力として受け取るため，複数ラベルは扱えてない．そのため，我々はグラフカットは異なる {0, 1, 2, . . . , k}N 上の離散関数の線形拡張を正則化として導入する．提案する離散関数はある特殊な状況で k-劣モジュラ関数と呼ばれる関数になる (Huber & Kolmogorov, 2012;

Kolmogorov, 2011)．また，提案した正則化学習問題を効率的に解く確率的最適化手法を提案する．提案する最適化手
法は，双対問題を確率的交互方向乗数法 (Suzuki, 2014) を用いて解くものである．双対問題の実効定義域の多面体的性質を利用することで，提案手法は線形収束することが示される．

2 多クラス連結正則化
G = (V, E, ℓ) を部分的にラベルの付けられたグラフとする．すなわち，V は頂点の集合，E は辺の集合であり，ℓ ∈ {0, 1}S は頂点の部分集合 S ⊆ V の上に与えられたラベルを表している．
我々の目的は部分的に与えられたラベル ℓ から V 全体の本当のラベルを当てることである．なお，ℓ は必ずしも正しい値であるとは限らず，ノイズが乗っていることも想定する．この問題を解くために，辺で連結された頂点のラベルは “近い” と仮定して，隣接頂点のラベルを近づけるような正則化をかける．そのような正則化の代表的な手法として一般化 Fused Lasso 正則化がある．これは次のような問題で記述される:∑
1∥xS − ℓS ∥2 + γ|xu − xv |,x∈RV 2min(1)(u,v)∈E

ただし，γ は正則化パラメータである．上の問題は連続最適化問題であるが，x を丸めることで離散的な解を得る．
上の問題はラベルが {0, 1} に二値だったが，k ≥ 3 値の問題も解けるよう，正則化を拡張する．再び G = (V, E, ℓ)を部分的にラベルの与えられたグラフとする．ただし，ℓ ∈ [k]S とする．この時，次の正則化学習問題を考え，これを多クラス連結正則化と呼ぶ:

γ1∥xS − ℓS ∥2 +2x∈RV ×[k] 2
min

k
∑ ∑

|xu,i − xv,i | +(u,v)∈E i=1

1

k
k
∑
γ(a − 1) ∑xu,i −
xv,i .
2
i=1
i=1

(2)

k = 1, a = 1 のとき，問題 (2) は問題 (1) と一致する．第三項のため，あるラベルに対応する座標が大きな値を取れば，他のラベルに対応する座標の値が小さくならなくてはいけない．このため，複数のラベルのうち，ただひとつのラベルが大きな値を取るような作用が働く．

3 双対問題とアルゴリズム
問題 (2) の目的関数は微分できず，単純な勾配法を当てはめることはできない．そこで，我々は双対問題を考え，その緩和問題を確率的交互方向乗数法で解く方法を提案する．F ∈ RE×V を隣接行列，すなわち (F x)e = xu,i − xv,iが x ∈ RV ×[k] に対して成り立つ行列とする．すると，問題 (2) の双対問題は

miny∈RE×[k]1γψ ∗ (y/γ) + ∥FS⊤ y − ℓS ∥2
2

(3a)

⊤
subject to F:,vy = 0 (v ∈/ S),

で与えられる．ここで，FC = [F:,v ]v∈C (C ⊆ V ) で，ψ ∗ はある関数 f ∗ を用いて ψ ∗ (y) =

(3b)
∑

e∈E f

∗

(ye,: ) と表され∗る. ここで，f (ye,: ) は集合 C1 と C2 を用いて次のように与えられる:{0f (y) =∞(y ∈ C1 + C2 ),(otherwise).

∗

ただし，C1 = {y ∈ Rk | ∥y∥∞ ≤ 1/2}, C2 = {y ∈ Rk | y = α1, |α| ≤ (a − 1)/2} である. この双対問題を
直接勾配法で解くことは難しいが，f ∗ が二つの単純な凸集合 C1 と C2 のミンコフスキー和であることを利用し，ye,: = y1,e + y2,e (y1,e ∈ C1 , y2,e ∈ C2 ) と分解することで，各変数 y1,e , y2,e に関しては f ∗ への射影が容易に計算できるようになる．このことを用いることで，確率的交互方向乗数法に基づいた線形収束する方法を構築することができる．

Theorem 3.1. x(t) を確率的交互方向乗数法の t ステップ目の解，x∗ を最適解とする．すると，ある ν > 0 と C > 0が存在して，x(t) は次のように R-線形収束する:

E[∥x

(t)

(

)
ν
− x ∥ ] ≤C exp −
t ,
|E|
∗ 2

収束レートの定数 ν は最適解を含むアフィン空間と実効定義域の角度で決まる．

参考文献
Bach, Francis. Learning with submodular functions: A convex optimization perspective. Foundations and Trends in
Machine Learning, 6(2-3):145–373, 2013.
Huber, Anna and Kolmogorov, Vladimir. Towards minimizing k-submodular functions. In Combinatorial Optimization,
pp. 451–462. Springer, 2012.
Kolmogorov, Vladimir. Submodularity on a tree: Unifying L♮ -convex and bisubmodular functions. In MFCS, pp. 400–411.
Springer, 2011.
Suzuki, Taiji. Stochastic dual coordinate ascent with alternating direction method of multipliers. In ICML, pp. 736–744,
2014.
Tibshirani, Robert, Saunders, Michael, Rosset, Saharon, Zhu, Ji, and Knight, Keith. Sparsity and smoothness via the fusedlasso. Journal of the Royal Statistical Society: Series B, 67(1):91–108, 2005.

2余震活動の確率予測：ベイズ統計を用いたアプローチ	 	
 東京大学生産技術研究所	
 
近江崇宏	
 
	
 
概 要 ： 大きな地震（本震）が起こると、その後に多くの地震（余震）が引き起こされる。大きな余震は被災地に追加の被害を及ぼしうることから、日本では気象庁が余震活動の確率予報を業務として行っている。その一方で、本震後に迅速に余震活動の予測を行うには、様々な技術的な問題が存在する。例えば、本震の直後では、多くの特に小さな余震が観測から抜け落ちてしまっている。そのため、本震直後ではこのような不完全な観測データに基づいて予測をしなければならない。また通常、観測データに基づき予測モデルのパラメータ値を一つだけ決めて、それに基づき予測を行っているが、短期間のデータからの推定には大きな不確定性が伴うため、そのような予測方式では予測が実際の観測から大きく外れてしまうことがある。我々はこのような問題を統計的な手法を用いて解決し、予測をより有用なものとするための研究をこれまで行ってきた。本講演
では、私たちのこれまでの研究を紹介し、余震活動の予測において、統計モデリングがどのように用いられ、いかに重要な役割を果たしているかを議論する。	
 	 余 震 活 動 の 統 計 モ デ ル と 確 率 予 測 ： 本震後の時刻 t におけるマグニチュード M の余震の発生率λ(t,M)は次のようなモデルで与えられる。

Kλ(t ,M) =e− βM .
p(t + c)!
ここで!K /(t + c)p は余震活動のべき的な時間減衰（大森—宇津則）を、!e

− βM

は地震頻

度のマグニチュードに対する指数減衰（グーテンベルク・リヒター則）をそれぞれ表している。
	  原理的にはこのモデルのパラメータ !{K ,p,c, β } を決めてやれば将来の余震活動の予測を立てることができる。その一方で、パラメータ値を決定するためには幾つかの難点がある。まず、余震活動はそれぞれの本震に応じて非常に個性が強く、たとえ本震のマグニチュードが同一だとしても、活動強度が全く異なることもある。そのため事前にパラメータ値を決めることが難しく、その時点で得られている対象となる余震活動の観測データからパラメータ値を直接推定する必要がある。しかしながら、特に本震直後の観測データは多くの観測漏れが含まれていることが知られており、不完全な観測データからパラメータ値の推定を行わなくてはならない。このことがこれまで、余震活動の迅速を難しくしていた。
	  このような問題に対して、私たちはデータの欠損を統計的に特徴づけるために、次のような時刻 t におけるマグニチュード M の余震の検出率関数!Φ(t ,M)M
1(x − µ(t ))2Φ(t ,M) = ∫ dx
exp[−
]
−∞
2
2σ 2
2
πσ
!
を導入した。ここで!µ(t ) は検出率が 50％に対応するマグニチュードを表し、時間変動すると仮定する。ここではベイズ平滑化により!µ(t ) を推定する[1]。この検出率関数と予測モデルを組み合わせることにより、不完全なデータからのパ
ラメータの推定が可能になり、本震直後からより精度の高い予測が行えるようになる[1,2]。
	  さらに予測行う際にはベイズ型の予測を導入した。これまでは最適なパラメータを一つ選び予測を行う（プラグイン予測）というようなことが行われてきたが、短期の余震観測データを用いた場合には、大きな推定の不確定性を伴う。
そのためプラグイン予測はしばしば観測値から大きく外れることがある。そこ
で、よりロバストな予測を行うため、MCMC によって事後分布からサンプルされた多数のパラメータ値からの予測を組み合わせるベイズ型予測を用いた。その結果、ベイズ型予測によって予測精度が優位に改善することが明らかになった[3]。
参考文献：	
 [1] T. Omi, Y. Ogata, Y. Hirata and K. Aihara, "Forecasting large aftershocks within one day
after the main shock", Scientific Reports 3, 2218 (2013).

[2] T. Omi, Y. Ogata, Y. Hirata and K. Aihara, "Estimating the ETAS model from an early
aftershock sequence", Geophysical Research Letters 41, 850 (2014).

[3] T. Omi, Y. Ogata, Y. Hirata and K. Aihara, "Intermediate-term forecasting of aftershocks
from an early aftershock sequence: Bayesian and ensemble forecasting approaches", Journal ofGeophysical Research: Solid Earth 120, 2561 (2015).

点過程モデルによる地震活動異常の解析と群発地震への応⽤
統計数理研究所

熊澤貴雄

統計数理研究所、東京⼤学

尾形良彦

気象研究所

⽊村⼀洋

気象研究所

前⽥憲⼆

気象研究所 ⼩林昭夫

地震活動の通常パターンからの異常はモデルのパラメータの時間的な変化として表れ、その変化は地震断層の状態が変化したことを反映する。断層の状態が変化する原因として近隣の地震性・⾮地震性の断層すべりに伴う応⼒場変化の他に、間隙流体圧の増加に伴う断層強度の低下等が報告されている。⾮地震性の断層すべりや間隙流体圧の増減は直接的な観測が困難だが、⼀定の震源域において地震が多発する群発性地震を引き起こし、通常の地震活動パターンから
統計的に外れる活動を⾒せる。これらの異常は短期的或いは断続的であり、時々刻々と変化する状態を記述するためには変化点的解析では不⼗分な場合が多く、より連続的にパラメータの時間変化を扱うモデルが必要となる。
報告者はETASモデル(Ogata, 1988)のパラメータの時間変化をベイズ平滑化により推定するベイズ型⾮定常ETASモデル（１）を導⼊し、(1)

λθ (t | H t ) = µ qµ (t) +∑K 0qK (ti )e ( i{i: S<ti <t}α M −M z )(t − t + c)pi$" q*,i+1 − q*,i$& N
q* (t) = ∑ I (t ,t ) (t) #(t − ti ) + qµ ,i ' = ∑ q*,i Fi (t) 、i i+1$% ti+1 − ti$( i=1i=1
N

ただし

2011年東北地⽅太平洋沖地震に誘発された内陸地震群を解析した。誘発地震の⼀部は本震の断層活動からは説明できないものであったが、誘発地震系列のパラメータの時間変化を推定することでそれらが間隙流耐圧の増加による断層強度減少に起因することを⽰すことがでた。
また、⽕⼭活動と地震活動の相関をモデル化する為に、⾮定常ETASモデルで推定したパラメータの時間変化と体積歪み計測量の関係を⽤いた予測モデルを

提案した。伊⾖地⽅東部でのマグマ貫⼊に伴う地殻の体積歪みは群発地震発⽣数と⾼い相関を⽰すことが報告されているが、モデルの或るパラメータの時間
変化がより⾼い相関を⽰すことが分かってきた。この関係をETASモデルに組み込むことで群発地震の短期的予測が可能となる。
より具体的には、ベイズ型⾮定常 ETAS モデル (１)でこれらの地震活動からインバージョンで求めた背景地震活動度µ (t) が体積歪の変動に約半⽇ほど遅れた相互相関を⽰すことが分かった。そこで先⾏時間の体積歪の変動から常時地震活動への線形回帰で予測すると、回帰係数は遡る先⾏時間に関して指数関数的に減少することが分かった。これより、先⾏する体積歪の変動時系列と指数関数の線形畳み込み（linear convolution）から背景地震活動を与える⾮定常
ETAS モデルの最尤推定値を求めると、どの群発地震でも殆ど同⼀な減衰係数となることが解った。しかし他⽅で、指数関数の切⽚パラメータは群発地震の開始位置から歪計までの距離に依存していることが分かった。すなわち群発開始位置を与えれば、歪データから正確な背景地震発⽣率を得ることができる。
結局⼆つの独⽴した解析結果が正当であることを裏付けられた。すなわち（1）ベイズ型⾮定常ETASモデルを使⽤した群発地震の背景活動度のインバージョン結果の精度が⾼いこと（2）群発地震の背景率とマグマ貫⼊の場所に歪観測点からの距離に依存する体積応⼒の時間変化との間に物理的に⾃然な因果関係がみられることである。マグマ貫⼊などの流体の作⽤によるストレス変化の定量化のメカニズムを理解する有⼒なヒントになりうると考える。

統計的言語モデルを用いた化学構造の特徴抽出と分子設計への応用総合研究大学院大学 池端 久貴 北陸先端科学技術大学院大学 本郷 研太(株)地球快適化インスティテュート 磯村 哲
北陸先端科学技術大学院大学 前園 涼 統計数理研究所 吉田 亮

1 背景
医薬品や材料を含む多く製品の開発現場において､計算機を利用した分子設計は非常に有用な手法である．データベースに含まれる既存の化学構造からの特性予測はQSAR(Quantitative Structure Activity Relationship)モデルと呼ばれ，これまでに多くの特徴抽出方法の提案，回帰モデルを始めとしたさまざまな統計的手法が適用されている[A. Varnek , and I. Baskin].しかしながら, 特性から化学構造の予測であるInverse-QSAR問題は，QSAR問
題に比べると研究成果の報告件数はかなり少ない．その理由としては化合物空間の広大さが考えられ，例えば，元素をC，N，O，Sに限定し，非水素原子の総数を30までとした場合でも，可能な有機分子の組み合わせ総数は1060個と膨れ上がる．
これまでInverse-QSAR問題に対して提案されている手法は２種類に大別できる．前者は特徴空間上のQSARモデルから得られた領域の内点に対応する化学構造を全列挙する方法で，後者は環構造を維持したまま化学構造を断片化して得られたフラグメントを構成単位とし，遺伝的アルゴリズムを用いて目的関数を最小化する化学構造を探索するといった手法である．
本研究では，統計的言語モデルと逐次モンテカルロ法を組み合わせた分子設計法を提案する．既存化合物から学習されたパターンを基に単語単位の生成モデルを構成できるため，これまでのフラグメントを用いた手法でカ
バーしきれなかった多様な化学構造を得ることができる．また，逐次モンテカルロ法を用いることで，遺伝的アルゴリズムや全列挙による手法ではできなかった確率的な取り扱いが可能となる．

2 手法
2.1

n-gramモデルによる化学構造のモデリング

SMILES形式は文字列で化学構造を表現するために広く用いられている形式の一つである.各元素は標準的な略記形が用いられ，有機分子骨格を形成するB, C, N, O, P, S, F, Cl, Br，I以外は[Fe]のように角括弧を用いる.環構造は始点と終点を数字で,側鎖部分は丸括弧で囲むことで表現される．このSMILES形式で表現された化学構造をn-gramモデルを仮定し，学習することで，化学構造に特徴的なパターンを抽出し，化学構造のランダム
サンプリングが可能となる．ただし，化合物の環構造，側鎖などの分岐構造を文字列で表現する際は離れた文字間の依存関係が現れるため，以下のような工夫を施す．


開始している環（ペアが現れていない数字）の個数，開始している側鎖の有無の条件毎に学習を行う
閉じた側鎖の情報を削減するオペレータ𝛹の導入 𝛹(CCC(CCCC)O) = CCC(C)O

以上を考慮した上で，SMILES文字列を𝒔 = (𝑠1 , ⋯ , 𝑠|𝒔| )，𝑎文字目から𝑏文字目までの部分文字列を𝒔[𝑎:𝑏] とすると
き，改良n-gramモデルは以下のように表すことができる．
|𝒔|

(1)

|𝐴|

𝑝𝑛 (𝒔) = ∏ ∑ 𝑝𝑛 (𝑘) (𝑠𝑖 | 𝛹(𝒔[𝑖−𝑛+1:𝑖−1] )) 𝐼(𝒔[1:𝑖−1] ∈ 𝐴𝑘 ) ,𝑖=1 𝑘=1

ここで𝑝𝑛 (𝑘) はグループ𝐴𝑘 における文字生成確率を表し，|𝐴|個あるグループへの割り当ては生成文字以前の部分
文字列の非ペアの数字の個数，閉じていない側鎖の有無から非確率的に決定される．パラメータ推定は既存化学構造からKneser-Ney smoothing [S. F. Chen et.al] などの推定方法を用いる．

2.2

Inverse-QSARモデリングある目標の特性𝑦 ∗ が与えられたときに，その構造𝒔についての条件付き分布𝑝(𝒔 | 𝑦 ∗ )を得たい.ベイズの定理を用いると𝑝(𝒔 | 𝑦 ∗ ) ∝ 𝑝(𝑦 ∗ | 𝒔)𝑝(𝒔)となり，右辺第一項はQSARモデルとなり,第二項は構造に関する事前知識となる．
これは式(1)で表される構造モデルに対応する． QSARモデルを構築する際には，化学構造を記述子と呼ばれるベクトルで表現されたものを入力とすることが多い．本稿では適当な記述子 𝜙: 𝐬 → 𝑥 を用いて𝑝次元のベクトル𝒙が得られたとする．目的の特性値𝑦は1次元実数値とし，未知の係数ベクトル𝒘を用いて，𝑦 = 𝒘𝑇 𝒙 + 𝜖 2 のように𝒙との間に線形の関係を仮定する．ここで𝜖 2 ~𝑖𝑖𝑑 𝑁(0, 𝜎 2 )とする．更に事前分布として，𝒘~𝑁(𝟎, 𝑉0 −1 )，𝜎 2 ~𝐼𝐺(𝑎0 , 𝑏0 )を仮定する．訓練データ𝐷 = {𝑋 = (𝒙𝑻1 , ⋯ , 𝒙𝑻 𝑁 )𝑇 , 𝒚 = (𝑦1 , ⋯ , 𝑦𝑁 )}が得られたとき,新しい入力𝒙∗ に対する特性値𝑦 ∗は次のようなnon-standardized Student's t-distributionに従う．
(2)𝑝(𝑌 = 𝑦 ∗ |𝒙∗ , 𝐷) = 𝑇 (𝒘𝑁 𝑇 𝒙∗ ,−1𝑏𝑁
(1 + 𝒙∗ 𝑇 𝑉𝑁 𝒙∗ ), 2𝑎𝑁 ),
𝑎𝑁

𝑁

1

2

2

ここで，𝒘𝑁 = 𝑉𝑁 𝑋 𝑇 𝒚， 𝑉𝑁 = (𝑉0 −1 + 𝑋 𝑇 𝑋) , 𝑎𝑁 = 𝑎0 + , 𝑏𝑁 = 𝑏0 + (𝒚𝑇 𝒚 − 𝒘𝑁 𝑇 𝑉𝑁 −1 𝒘𝑁 )である．
通常、目標の特性は[𝑦1 , 𝑦2 ] のように区間を取ることが多い．このとき𝒔の事後分布は，以下のように表され
る．
(3)

𝑦2

𝑝( 𝒔 | 𝑦1 ≤ 𝑌 ≤ 𝑦2 ) ∝ ∫ 𝑝(𝑌 = 𝑦 ∗ |𝒔, 𝐷)𝑑𝑦 ∗ 𝑝𝑛 (𝒔)𝑦1

2.3

逐次モンテカルロ法得られた𝒔の事後分布を目標分布とし，逐次モンテカルロ法 [P. Del Moral et.al] を適用することで式(3)に対応した化学構造の事後分布を得る．ステップ t における目標分布𝛾𝑡 を
𝛾𝑡 = {𝑞(𝒔)}𝜅(𝑡) 𝑝𝑛 (𝒔).
(4)𝑦とする．ここで，𝑞(𝒔) = ∫𝑦 2 𝑝(𝑌 = 𝑦 ∗ |𝒔, 𝐷)𝑑𝑦 ∗ で，𝜅(𝑡)は正の単調関数で少なくとも最終ステップにおいては
1

𝜅(𝑡) = 1 となる関数である．逐次モンテカルロ法の手続きとしては以下のようになる．
1. initialize
2. at time t
3. Local move from 𝒔𝑡−1 (𝑖) to 𝒔𝑡 (𝑖) with the following procedure (𝑖 = 1, ⋯ 𝑁)
 obtain 𝑢~𝐵𝑖𝑛𝑜𝑚(𝐿, 0.5) for prefixed 𝐿.
 reduce u letters from the right edge and generate 𝐿 − 𝑢 letters from the right edge with probability inthe learned probability from the equation(1).
4. update weights as 𝑤𝑡 (𝑖) ≔ 𝑤𝑡−1 (𝑖)
1

2 −1

(𝑖)
5. if (( ∑𝑁)𝑖=1 𝑊𝑡𝑁{𝑞(𝒔𝑡 (𝑖) )}𝜅(𝑡){𝑞(𝒔𝑡−1 (𝑖) )}𝜅(𝑡−1), for 𝑖 = 1, ⋯ 𝑁, and normalize them to obtain 𝑊𝑡 (𝑖)≤ 𝑁⁄2), then resample with weights as the probability, and 𝑊𝑡 (𝑖) =1
𝑁

6. Reordering the randomly chosen atom being at the first position by some software such as OpenBabel.
7. n:=n+1, then go back to 2.

3 検証
材料開発への適用例を発表にて報告する．

参考文献[1] A. Varnek , and I. Baskin, Machine Learning Methods for Property Prediction in
Chemoinformatics: Quo Vadis, Journal of chemical information and modeling, 52, 1413-1437, 2012.
[2] S. F. Chen and J. Goodman, An Empirical Study of Smoothing Techniques for Language Modeling, Computer Speechand Language, 13, 359-394, 1999.
[3] P. Del Moral et.al, Sequential Monte Carlo samplers, J. R. Statist. Soc. B, 68, 411-436, 2006.

Bay-Bridge : local linear approximated bridge and its modelbuilding procedure via the Bayesian model
中央大学大学院理工学研究科

保科

架風

1 はじめに
近年の様々な技術の発展に伴い, 多くの場面で大規模かつ複雑な構造を有するデータが獲得・蓄積されている. そのようなデータから線形回帰モデリングによって現象のメカニズムの特定や予測を行う際, 最尤・最小2 乗推定やモデル選択基準に基づく変数選択が適用されてきた. しかし, これらの手法では推定の精度や計算コストの面で問題が生じることがある. これに対し, モデルの推定と変数選択を同時に行うスパース回帰モデ
リングでは, このようなデータを効率的に取り扱うことが可能であり, この 20 年間盛んに研究が進められ, 飛躍的な発展を遂げている.
スパース回帰モデリング手法の 1 つである bridge 回帰 (Frank and Friedman, 1993) は 2 つの調整パ
ラメータの値を適切に選ぶことで回帰係数ベクトルの幾つかの成分を厳密に 0 と推定することが可能であり, 様々なデータから有益な情報を抽出することができる. これに対し, Kawano (2014) では GBIC 基準(Konishi et al., 2004) によって bridge 回帰の調整パラメータを選択する手法を提案している. しかし, 複数の調整パラメータの存在はモデル選択のコストを増大させる. また, bridge 回帰の推定では非凸最適化問題を扱う必要がある.
これに対し本研究では, 1. bridge 回帰の回帰係数ベクトル, 誤差分散, 1 つの調整パラメータをベイズモデルに基づく MAP 推定を行うことでモデル選択の次元を下げ, 2. 残りの調整パラメータの値を選ぶためのベイズモデルに基づくモデル選択基準の導出を行った. 一般的にベイズモデルに基づくモデル選択基準はサンプ
ルサイズに依存するラプラス近似を使用しており, 高次元データに対してしばしば有効に機能しなくなる. これに対し本研究におけるモデル選択基準の導出では近似にモンテカルロ積分を適用し, また, 高次元モデルでの有効性が知られている EBIC (Chen and Chen, 2008) 基準の考え方を取り入れた. これにより, 提案手法
は高次元データのモデリングにおいても有効に機能することが期待できる.

2 提案手法
目的変数 y, p 次元説明変数 x = (x1 , . . . , xp )T に対するサンプルサイズ n の線形回帰モデル

y = β0 1n + Xβ + ε(1)を考える. ただし, n 次元目的変数ベクトル y = (y1 , . . . , yn )T , 切片項 β0 , すべての成分が 1 の n 次元ベク
トル 1n , n × p 計画行列 X = (x1 , . . . , xn )T (xi = (xi1 , . . . , xip )T , i = 1, . . . , n), p 次元回帰係数ベクトル
i.i.d.

β = (β1 , . . . , βp )T であり, ε = (ε1 , . . . , εn )T を n 次元誤差ベクトルとし, εi ∼ N (0, σ 2 ), σ 2 を未知の誤差
∑n
∑n
∑n
2
分散とする. また, 一般性を失うことなく i=1 yi = 0,i=1 xij = 0,i=1 xij = n, (j = 1, . . . , p) と目的
変数と説明変数を基準化する.
このとき, bridge 回帰 (Frank and Friedman, 1993) は次の Lq 正則化法として定義が与えられる:β̂ lasso := argmin ∥y − Xβ∥2 + λ
β

p
∑
j=1

1


|βj |q  .

(2)

q ≤ 1 のとき, bridge 回帰は lasso と同様にスパース性を持つが, q < 1 のときの推定では非凸最適化問題を解く必要がある. また, bridge のモデル選択では λ と q の最適な組み合わせを選ぶ必要があり, 候補のモデルの数が増大する.
これに対し, 本研究では bridge にベイズモデルの視点から β, δ = 1/σ 2 , λ2 を MAP 推定する以下の

Bay-Bridge (Bayesian MAP estimation with local linear approximated bridge regression) を提案した:


p
∑
√
n−2
δ
p + νλ − 2ηλ(β̂, δ̂, λ̂2 ) := argmin −log δ + ∥y − Xβ∥2 + λq δ
log λ2 + λ2  . (3)wj |βj | −2222β,δ,λ2j=1(0)

ただし, wj = |βj |q (j = 1, . . . , p), (νλ /2, ηλ /2): λ2 のガンマ事前分布のパラメータである. また, この

MAP 推定では bridge に local linear approximation (Zou and Li, 2008) を適用したものを基にしている.
さらに本研究では Bay-Bridge に対するモデル選択基準として EBIC 基準の周辺尤度の近似にモンテカル
ロ積分を適用し, 以下のモデル選択基準を導出した:

(

)

(ν )
ηλ
λ
+ 2 log Γ
2
2
( )
(
νλ )
1
∗
− 2 log Γ p +− log |W | + log γ + Γ+ 2 log M22}−(n−p∗ )/2
{
M
∑
1 T
−1
T
−1/2
− 2 logy (In − XAq A(m) XAq )y|A(m) |
2
m=1
 


−p∗ +(νλ /2)
( )
1 ∑

∑
p
2
2




γτj(m) + 2α ∗ .
τj(m) + ηλ
exp
2

p
∗

EBIC(q) =n log 2π + 2p log 2 − log Γn − p∗2− νλ log(4)β̂j ̸=0β̂j ̸=0

ただし, p∗ = |Aq |, Aq = {j | β̂j ̸= 0, j = 1, . . . , p}, XAq は Aq に対応する説明変数で構成される計画行列,−1T
A(m) = XA
XAq + W D(m), W = diag(ξj2 |j ∈ Aq ), ξj = qwj (j = 1, . . . , p), D = diag(τ12 , . . . , τp2 ), D(m) :q∏22
対角成分に {τj(m)| j ∈ Aq } を持つ対角行列, {τj(m)| j ∈ Aq , m = 1, . . . , M }: j∈Aq Gamma(τj2 |1/2, γ)からのサイズ M のランダムサンプルである. これにより, 高次元データに対するモデリングにおいて有効に機能するモデル選択基準を導出することができた.

参考文献[1] Chen, J. and Chen, Z. (2008). Extended Bayesian information criteria for model selection with large
model spaces. Biometrika, 95, 759–771.
[2] Frank, I. and Friedman, J. H. (1993). A statistical view of some chemometrics regression tools (withdiscussion). Technometrics, 35, 109–135.
[3] Kawano, S (2014). Selection of tuning parameters in bridge regression models via Bayesian information criterion. Statistical Paper, 55, 1207–1223.
[4] Konishi, S., Ando, T., and Imoto, S. (2004). Bayesian information criteria and smoothing parameterselection in radial basis function networks. Biometrika, 91, 27–43.
[5] Zou, H. and Li, R. (2008). One-step sparse estimates in nonconcave penalized likelihood models.
The Annals of Statistics, 36, 1509–1553.

2

潜在的歩数増加パターンと体組成指標との統計的関連性についての分析大橋洸太郎（立教大学社会情報教育研究センター・慶應義塾大学システムデザインマ
ネジメント研究所）

小熊祐子（慶應義塾大学大学院健康マネジメント研究科・スポーツ医学研究センター）渡辺美智子（慶應義塾大学大学院健康マネジメント研究科）
1．目的
本研究は，生体ログデータの計量化のために混合成長曲線モデルを適用する一連の統計的分析処理技法を提案し，その適用例として分析対象を年代や性別によって層別した上で，「歩数」特に「累積歩数」に関する潜在的なパターンの相違が体組成に影響するプロセスモデルを例示することを目的としている。厚生労働省(2013)では，がん，循環器疾患，糖尿病や慢性閉塞性肺疾患を生活習慣病と位置づけ，日本国民のこれらの疾患への対策として，食生活の改善や運動習慣の定着に根差した一次予防と呼ばれる分野に重点を置いた対策が望まれていることを述べている。このような中で今回我々は，混合成長曲線モデル(Growth Mixture Model，GMM)の応用と，その結果に共分散分析を組み合わせることが，運動に伴う生体ログデータの解析に有効ではないかと考えた。また適用例として，独立行
政法人情報通信研究機構(National Institute of Information and Communications
Technology，NICT)の採択課題「課題 A ソーシャル・ビッグデータ利活用アプリケーショ
ンの研究開発」における「ヘルスリテラシー向上のための生体ログデータ分析に基づく健康情報フィードバック」プロジェクトで収集が行われている「歩行」と「体組成指標」のデータを用いた。

2．方法
プロジェクトへの協力者の内，性別と年齢の情報と，4 週間分 28 日以上の歩行記録が揃っていた 2830 名(男性 1757 名，平均 48。38 歳(SD=11。86)，女性 1073 名，平均 41。59 歳(SD=12。
29))を分析対象とした。また協力者を人数の多かった層を中心に「30 代男性」「30 代女性」「40代男性」「40 代女性」「50 代男性」「50 代女性」の 6 層に分けた。続いてこれらの各層の協力
者の 1 週目から 4 週目までの「累積歩数」に関して GMM を適用し，各層が何群に分けられるのかを検討した。そして分類による効果を吟味するために，分類された各「群」と「経過週数」を説明変数とする共分散分析を行った。検定の目的変数となる体組成指標には，株式会社タニタの測定器による「内臓脂肪レベル」，野村・渡辺・小熊(2015)で提案された「体型指標」「隠れ肥満指標」の 3 つを用いた。「経過週数」は 1 週目の累積歩行数と 2 週目の体組成の指標の平均値の組合せを「1 週経過」，2 週目の累積歩行数と 3 週目の体組成の指標の平均値を「2 週経過」，
3 週目の累積歩行数と 4 週目の体組成の組合せを「3 週経過」として 3 時点のデータを作成し，前の週までの累積歩行数が次の週の体組成にどれほど影響しているのかを検討した。

3．結果

GMM の結果については，50 代男性と 50 代女性について Table1 に掲載した。AIC による群数の比較の結果，50 代男性，50 代女性共に 5 つの群に分類することが適切であることが分かった。

Table 1：GMMによる群数
50代男性 50代女性
群数
AIC
AIC
4
30208.51 11458.79
529901.77 11306.53
629907.77 11312.53

Table 2：50代男性の「体型指標」の群毎の平均値
50代男性
群番号
体型指標 隠れ肥満指標 内蔵脂肪レベル
5
1.403
0.743
9.73
4
1.912
0.160
11.1113
3.101-0.04512.7892
3.298-0.271
12.8231
4.193-0.82513.991

Table2 に 50 代男性の「体型指標」「隠れ肥満指標」「内臓脂肪レベル」の群毎の平均値を示した。群番号が大きい程多く歩いていた群を指す。これらの指標は値が小さい程体組成が良いと判断されるものであり，この結果から多くの歩数を記録した群程体組成指標もよいとわかった。
共分散分析の結果については，50 代男性の「体型指標」の結果を Table 3 に掲載した。Table 3
における身長，体重，基礎代謝，筋肉量は共変量である。これらの変数を統制した上で 50 代男
性は，「群」「経過週数」による体型指標の変化に有意な差がみられた。またこの他 50 代男性には，「隠れ肥満指標」「内臓脂肪レベル」に「群」による有意差がみられていた。
Df
群
経過週数
身長
体重
基礎代謝
筋肉量
群×経過週数
残差

4
2
1
1
1
1
8
570

Table 3 ：50代男性「体型指標」
Sum SqMean Sq
F value
421.0
105
1.579E+03
1.0
18.121E+00
0.0
02.750E-01
6993.06993
1.049E+05
1.0
12.027E+011.0
11.704E+010.0
04.500E-02
38.0
0

Pr(>F)
2.000E-16 ***3.330E-04 ***5.990E-01
2.000E-16 ***8.170E-06 ***4.210E-05 ***9.990E-01

4．考察
4 週間，1 ヶ月程度の歩行記録をつけ，歩いていくことで「体型指標」は有意に変化していく
ことが分かった。また，群による違いは「体型指標」「隠れ肥満指標」「内臓脂肪レベル」にみられ，より多く歩いた記録を持つ群の方が，そうでない群よりも体組成指標が改善する様子が統計的に確認された。ただし今回分析対象としたデータは，記録のあった協力者の測定結果であり，
この結果の中には測定機器の非装着による歩行数の記録の未計上分が含まれていると考えられる．例えば日に 10 時間以上の装着時間のあるものを有効データとするといった，データの峻別を視野に入れ，今後の分析を行いたい
5．文献
・厚生労働省，2013，二十一世紀における第二次国民健康づくり運動 (健康 21(第二次)) ，厚生
労働省 HP<http://www1.mhlw.go.jp/topics/kenko21_11/b2.html#A22>(最終閲覧日，2016 年 1 月 7 日).
・野村俊一・渡辺美智子・小熊裕子，2015，体組成指標の推移パターンの解析と適正な減量の判
別，第 6 回横幹連合コンファレンス論文集，pp223-226.
・タニタからだカルテ HP <http://www.karadakarute.jp/tanita/column/columndetail.do?columnId=172>(最終閲覧日 2016 年 1 月 7 日)

生体ログデータを使用した歩行パターンと体組成指標との関連を検討する分析方法の提案
渡辺 真弓 1
1

渡辺 美智子 1

慶應義塾大学健康マネジメント研究科

2

小熊 祐子 1,2
スポーツ医学研究センター

E-mail: yamam00@sfc.keio.ac.jp

1.

はじめに

身体活動の増加は，将来的な疾病の予防という効果に留まらずメンタルヘルス・自己効力感の向上，腰痛・膝痛軽減，上気道感染予防等に寄与することが知られている．一方，身体活動の不足は生活習慣病発生の危険因子である 1)．近年，スマートフォン等の普及等によって，日常生活中の情報をデジタルデータとして蓄積できる環境が整ってきた 2)．本研究の目的は，歩数計や体重計において蓄積された生体ログデータの活用方法の提案である．身体活動に関する情報を用いて歩行パターンを導出し，健康状態との関連を検証する一連の過程を提示する．
2. 方法
歩数データには，活動量計で記録された 232300 件，3812 人分の 2013 年 4 月 1 日から 2014 年 3 月 31
日までの一年間の歩行データが時間単位で保存されている．まずは各データの ID 平均からの偏差を求め，
この偏差を第 1～3 四分位点で 4 分割した．ただし，0 時から 5 時の間は 2 分割した．このデータを用いて，歩行パターン把握のために潜在クラス分析を行った．クラス数は AIC，BIC の情報量基準を参考にして決定した．分析には，MplusVersion7.11 を使用した．
体組成データは，歩行データと同じく 2013 年 4 月～2014 年 3 月の 1 年間に体組成計にて測定された，
男性 361772 件（5771 人），女性 213102 件（5354 人）のデータである．体重，体脂肪率，筋肉量，基礎
代謝量，体内年齢のデータを使用し男女別に主成分分析を行った．分析には，JMP12.0.1 を使用した。
3. 結果
3.1 歩行パターン
AIC,BIC の値と解釈妥当性を考慮して，
4 クラスを選択した．
各クラスの構成割合は，
クラス 1 が 22.1%，
クラス 2 が 21.6%，クラス 3 が 28.4%，クラス 4 が 27.6%だった． 応答確率から，クラス 1 を「休日型アクティブ日」，クラス 2 を「休日型非アクティブ日」，クラス 3 は「通勤型アクティブ日」，クラス 4は「通勤型非アクティブ日」とした．各クラスへの所属確率の ID 平均を算出し，最も高い所属確率を示すクラスをその者のクラス（歩行パターン）とした．また，平均歩数を 3 分割し，歩行パターンと平均歩数を掛け合わせて，12 類型を求めた． 12 類型毎に，4 つの歩行パターンに対する所属確率割合の月次推
移を Fig.1 に示した．
休日型アクテ 休日型非アク 通勤型アクテ 通勤型非アク
3.2 体組成パターン
ティブ日多
ィブ日多
ティブ日多
ィブ日多
体組成データの主成分分析の結果,主成分
歩
1 は総合的な身体の大きさ，中でも太さに強数
少く関連する指標であると考えられた．主成分
2 は，身体がどの程度引き締まっているかを示す指標である考えられた．
歩
数
3.3 歩行パターンと主成分パターンの結合中主成分1 と主成分 2 の主成分得点を算出し，歩行データと結合させた．この結果，男性 690
歩人，女性 290 人のデータを得た．歩行パター
数
多ン 12 類型毎に主成分 1 と主成分 2 の月次推移を示したのが Fig.2 である．男性においては,平均歩行数が「少」で「休日型アクティ「休日型アクティブ日多」所属確率
ブ日多」「通勤型非アクティブ日多」の者は「休日型非アクティブ日多」所属確率
5 人以下だったため，分析からは除外した．
「通勤型アクティブ日多」所属確率
総合的な身体の大きさ（太さ）は，平均歩数「通勤型非アクティブ日多」所属確率が増加するほど下がり，痩せる傾向にあるとFig.1 Posterior probabilities for each walking
言える．しかし，平均歩数が同程度であってpattern by average walk steps

も，どの歩行パターンが多いかということで主成分 1 の得点は異なり，通勤型非アクティブ日が多い者の得点が低く，痩せている傾向にあると言える．身体の引き締まり度は，やはり平均的な歩数が多くなるほど低くなる傾向にあり，身体が引き締まっていることが分かる．同程度の平均歩数であっても，特に「休日型アクティブ日」が多い者と「通勤型アクティブ日」が多い者で低く，身体が引き締まっている傾向にある．主成分 1 の得点が低く，総合的には痩せている傾向にある平均歩数「多」で「通勤型非アクティブ
日多」の者は，引き締まり度合いはあまり高くないという結果であった．
次に女性においては，平均歩数「少」で「休日型アクティブ日多」「通勤型非アクティブ日多」の者，
平均歩数が「多」で「通勤型アクティブ日多」の者は 1-2 人と非常に少数だったため，分析からは除外した．総合的な身体の大きさは，平均歩数が多くなるほど低くなり，痩せる傾向にある．同程度の平均歩数で比較すると，歩行パターンによる大きな違いは見られない．身体の引き締まり度は，平均歩数が少ない
者で低く，身体が引き締まっている傾向にある．これは，一般的な認識とは逆の傾向であるが，データ数の少なさに起因することも考えられる．
Male

Female

休日型アクテ 休日型非アク 通勤型アクテ 通勤型非アク

休日型アクテ 休日型非アク 通勤型アクテ 通勤型非アク

ィブ日多

ティブ日多

ィブ日多

ティブ日多

歩
数
少

ィブ日多

ティブ日多

ィブ日多

ティブ日多

歩
数
少
歩
数
中

歩
数
中

歩
数
多

歩
数
多

主成分 1：身体の大きさ（太さ）主成分 1：身体の大きさ（太さ）主成分 2：引き締まり度

主成分 2：引き締まり度

Fig.2 principal component scores by each walking patterns

4.

考察

平均的な歩数が多くなるほど，身体は痩せ，引き締まる傾向にあった．これは，歩行を行うほど筋肉が増加し健康に良いと言った一般的な認識を裏付けるものであると言える．ただし，男性においては同程度の平均歩数でも，歩行パターンによって体組成は異なる傾向を見せた．「休日型」「通勤型」の双方で，アクティブな日を多く過ごす者では，そうでない者と比べてより引き締まり度合いが高まり，健康的な身体状態であると考えられた．平均的に同程度の歩数であっても，歩行パターンが異なることで身体への影
響も異なることが示唆される．また，外見的にはスリムでも，筋肉量が少なく不健全な痩せ方をしている者の存在が示唆され，今後の検討に値する．ただし，女性においてはこのような傾向は見られなかった．
むしろ平均歩数が少ない程身体が引き締まるという，解釈が困難な結果となったが，この点に関してはデータ数の少なさが大きく影響している可能性を否定できず,今後データ数を増加した上で再検討する必要がある．
謝辞
本研究は，独立行政法人情報通信研究機構（NICT）の委託研究「ソーシャル・ビッグデータ利活用・基盤研究の研究開発」における採択課題「ヘルスリテラシー向上のための生体ログデータ分析に基づく健康情報フィードバック」の支援を受けたものである．
参考文献
1) 厚生労働省，健康づくりのための身体活動基準 2013（2013）2) 西山勇穀他，ライフログデータを用いたチームの行動変容促進，情報処理学会論文誌，56 巻 1 号，349/361(2015)

体型を示す新たな体組成指標の提案とそのダイナミクス
野村 俊一 1,2 渡辺 美智子 3 小熊 祐子 3,41

東京工業大学情報理工学研究科 2 慶應義塾大学システムデザイン・マネジメント研究科
3

慶應義塾大学健康マネジメント研究科 4 慶應義塾大学スポーツ医学研究センター
E-mail: nomura@is.titech.ac.jp

１

はじめに
本稿では，成人男女の体重および体脂肪率の分布に基づいて構成した，体型および隠れ肥満度を判定するための２つの指標を提案する．通常の肥満には該当しない体重であっても，提案指標により体脂肪量の多寡を基に隠れ肥満の傾向にあるかどうかを判断することができる．さらに，１年間にわたる個々人の体重および体脂肪率の計測記録を分析し，２つの提案指標のダイナミクスの違いを明らかにすることで，両指標を改善する減量ほど長期的に継続することが示された．また，提案指標について都道府県別に集計した結果を示し，メタボリックシン
ドロームの該当率と比較しながら地域的な体型の傾向を考察する． なお，本稿における解析に用いるデータは，2013 年 4 月 1 日から 2014 年 3 月 31 日までの 1 年間にわたって継続的に計測された，個人別の体重および体脂肪率のログデータである．記録蓄積された体組成計測定データを分析して見出された体組成推移パターンに共通する性質を活用することで，個人ごとの健康状態や目的に応じた精度の高い健康指導情報の提供を目指す．
２

体型判定への新たな指標の提案
現在，肥満度の判定基準として最も利用されている指標に，BMI (Body Mass Index) ＝体重÷身長 2 (kg/m2)がある．また，BMI と同様に，体脂肪率も肥満の判定に広く用いられている．BMI が健常な範囲内であっても体脂肪率が高いと，内臓脂肪が過剰な隠れ肥満の兆候があり，逆に体脂肪率が低いとアスリートのような筋肉質な体型といえる． BMI と体脂肪率は、両者の物理的単位が異なることから相関関係が非線形で捉えづらい．そこで，BMI を次式のように体脂肪にかかる部分（Fat Mass Index; FMI＝BMI × 体脂肪率）と体脂肪以外にかかる部分（Lean Mass Index; LMI＝BMI × (1－体脂肪率)）に分解する．FMI と LMI は同じ物理的単位 (kg/m2)
をもっており，体重に占める体脂肪量とその他の重量を，BMI のように身長に関して標準化した数値といえる．
ここで，FMI と LMI の変数の組に対して，男女ごとに主成分分析を行うと，そ
Male
Male

れぞれ Fig. 2 の矢印のようにばらつきの大きい第一主成分と，それに直交するば（Body Size Index）および隠れ肥満指標（Hidden Obesity Index）を，体型分類のための新たな指標として提案する．
体型指標（男性）= 1.10 × FMI + 0.85 × LMI－19.65

2020
BMI(bodyfat)
BMI(bodyfat)
5 5 1010 1515

主成分をスケール変換および平行移動することにより得られる下式の体型指標

Body size
Body
size
index
index

Hidden obesity
Hiddenobesity
index
index

0
0

00

らつきの小さい第二主成分とに分けることができる．この第一主成分および第二1010

隠れ肥満指標（男性）= 0.85 × FMI－1.10 × LMI + 16.70

15
20
25
15 BMI(others)
20
25
BMI(others)

FemaleFemale肥満指標が一定の下，体型指標が 1 単位増えることは BMI が 1 単位増えること

BMI(bodyfat)fat)BMI(body
5 5 1010 1515 2020に相当している．従って，体型指標の符号は体重が標準に比べて過剰であるか否0 0

体型指標（女性）= 1.20 × FMI + 0.35 × LMI－13.10
隠れ肥満指標（女性）= 0.35 × FMI－1.20 × LMI + 16.50
Fig. 1 の矢印は上式で提案した指標の座標軸を表しており，座標軸の原点，すなわち両指標が 0 である状態は，BMI の標準値 22 に対応している．また，隠れかを示しており，体型指標の正負から「痩せ型」か「肥満型」に体型を大きく分類

30
30Body size index
Body size index

HiddenHiddenobesityobesity
index
index

0
0
10
10

15
20
25
15 BMI(others)
20
25
BMI(others)

30
30

することができる．一方，隠れ肥満指標は，体重の内訳を重視した指標であり，体Fig. 1: Scatter plots of

重に占める体脂肪の比率が高くなるほど隠れ肥満指標は大きくなる．隠れ肥満指FMI and LMI with axes of

標は測定者の集団で平均 0 となるよう中心化されており，隠れ肥満指標の正負に

proposed indices.

より「筋肉質型」と「隠れ肥満型」への分類が可能となる．隠れ肥満は外観からわかりづらいため，本人が自覚しないまま進行することが多く，それゆえ，この隠れ肥満指標は，隠れ肥満の進行度を表すバロメータとして有用であると考えられる．
３

提案指標の地域的傾向上で提案した 2 つの指標について，計測者の居住都道府県別・男女別・年代別に集計を行い，提案指標の地域的な傾向を探る．Fig. 2 では 40 歳以上 59 歳以下男性の体型指標および隠れ肥満指標の平均値に加え，さらに，比較対象として厚生労働省の資料より集計したメタボリックシンドローム該当者の受診者に対する割合および該当者＋予備群該当者の割合を示している．全体的に見て，東北地方～北海道と，四国地方～九州地方にかけて，各指標の平均値およびメタボリックシンドローム該当割合が高い傾向にあることが分かる．男性の両指標の大きさとメタボリックシンドローム該当割合の高低は概ね同じ傾向を示している以上から，提案指標はメタボリック
シンドロームと一定の整合性が取れた指標と見ることができるが，２つの指標の地域性には違いも見られ，体型の類型により細かい地域性があるように考えられる．
４

提案指標の月別推移
ここでは，上で提案した 2 つFig. 2: Cartograms of proposed indices and rates of Metabolic syndrome (%).

の指標について，年間を通した
各月における指標の推移を分析し，指標ごとの推移の仕方に現れる特徴を探る．分析するデータはここまでと同
じ 2013 年 4 月から 2014 年 3 月の 1 年間の体重および体脂肪率の計測記録であるが，データの精度を確保するため，毎月 10 日以上にわたり計測が行われた男性 273 名，女性 97 名のデータに対象を絞った．体型指標と隠れ肥満指標それぞれの月次推移に有意水準 1%の単位根検定（Harris-Tsavalis test [1]）を行うと，男女ともに隠れ肥満指標の推移は定常性をもつという結論が得られる．つまり，隠れ肥満指標は，一時的に変化しても，長期的には各人の生活習慣に見合った水準へと落ち着くという平均回帰性をもっているといえる．一方，体型指標について同じ単位根検定を行っても定常性をもつという結論は得られず，一度増減した体型指標は自然と元には戻るものではなく，戻すための働きかけが必要になることを示唆している．
体重・体脂肪率の計測記録の平滑化
不定期かつ誤差を含む体組成計測記録から，精度の高い体重・体脂肪率の状態推移を得る手法として，ここではカルマンフィルターによる計測値の平滑化を提案する．Fig. 3 のように，体重および体脂肪率の実

Weight (kg)
59 60 61 62 63 64

５0

マンフィルターによる平滑化推定値は安定した推移を示しており，計測ミスの影響も抑えられている．
謝辞

100

200
Days

300

Body fat percentage (%)
16
18
20
22

計測値は，大きく上下に振れているのに対して，カルFig. 3: Time series of weight with Kalman smoother.

本研究は，独立行政法人情報通信研究機構（NICT）の委託研究「ソーシャル・ビッグデータ利活用・基盤研究
0
100
200
300
の研究開発」における採択課題「ヘルスリテラシー向上のための生体ログデータ分析に基づく健康情報フィード
Days

バック」の支援を受けたものである．
参考文献[1] Harris, R. D. F., and Tsavalis, E. Inference for unit roots in dynamic panels where the time dimension isfixed, Journal of Econometrics, 91, 201--226 (1999).