関数データに基づく半教師あり識別・判別問題

1

九州大学大学院数理学府

川野秀一

九州大学大学院数理学研究院

小西貞則

はじめに

識別・判別問題において，どの群に属するか既に分かっているデータ (ラベルありデータ) と群の所属が未知のデータ (ラベルなしデータ) が混在している状況の下，ラベルありデータのみに基づいて判別モデルを構成するのではなく，ラベルなしデータを効果的に用いることにより判別モデルを構成し，予測の精度を向上させる手法は，半教師あり学習法と呼ばれている (Chapelle et al., 2006)．近年この学習法は大きな注目を集めており，
学習データの取得が高コストである分野，特にバイオインフォマティックスやテキスト分
類などの分野で応用されつつある (例えば，Nigam et al., 2000; Chapelle and Zien, 2005;
Dean et al., 2006; Pan et al., 2006 参照)．
一方，離散時点で経時的に測定・観測されたデータに関数化処理を施し，得られる関数
データ集合を解析対象とする手法は関数データ解析と呼ばれ，諸科学のさまざまな分野でその有用性・重要性が報告されている (Ramsay and Silverman, 2002; 2005)．しかし，群の所属が既知の関数データ (ラベルあり関数データ) と未知の関数データ (ラベルなし関
数データ) の両方の情報を用いた判別モデルに関する研究，つまり，関数データ集合に対する半教師あり学習法の研究についてはいまだ十分に確立していない．
本報告では，各個体の観測データに対して関数化処理を施し，関数データ集合を対象とした半教師ありロジスティック識別・判別法を提案する．モデルに含まれるパラメータの推定には，モデルの安定化を図るために正則化法を適用し，推定アルゴリズムに EM アル
ゴリズムを用いる．また，モデルに含まれるチューニングパラメータの値を客観的に選択するために，ベイズアプローチに基づくモデル評価基準を導出する．さらに，実データへの適用を通して提案手法の有効性を検証する．

2

関数データ

いま，n 個の個体に対して特性 X(∈ R) について変数 T (∈ R) をともない観測された離散的なデータ集合{(tαj , xαj ); j = 1, . . . , Nα , tαj ∈ T (⊂ R)},(α = 1, . . . , n),(1)が与えられたとする．ただし，xαj は α 番目の個体の時点 tαj での観測データとする．この離散データに対して非線形回帰モデルを当てはめることにより，関数データ集合 {xα (t); t ∈T } (α = 1, . . . , n) を求める．
観測データ集合 {(tαj , xαj ); j = 1, . . . , Nα } は，次のような回帰モデルに従うと仮定する．

xαj = uα (tαj ) + εαj ,1

j = 1, . . . , Nα .

(2)

2.5
2.0
1.5
1.0
0.5
0.0
0.0

0.2

0.4

0.6

0.8

1.0

図 1: 基底展開法による離散データの関数化の例．丸印は離散データ，図下部の曲線は基底関数，図上部の曲線は関数データを表す．
ここで，uα (·) は未知の平均構造，εαj は平均 0，分散 σx2α の正規分布に従うものとする．未知の平均構造 uα (·) は，基底関数と呼ばれる m (< Nα ) 個の既知の関数 φi (t) (i = 1, . . . , m)の線形結合により，

xαj =

m

i=1

cαi φi (tαj ) + εαj = cTα φ(tαj ) + εαj(3)と表されると仮定する．ただし，cα = (cα1 , . . . , cαm )T は m 次元係数ベクトル，φ(t) =(φ1 (t), . . . , φm (t))T は m 次元基底関数ベクトルとする．
基底関数としては，これまでデータの特性・分析目的に応じて，フーリエ級数，B-ス
プライン (Imoto and Konishi, 2003)，ウェーブレット (Fujii and Konishi, 2006) などさまざまなものが用いられてきた．本報告では，基底関数としてガウス型基底関数を用い，特に次で与えられるガウス型基底関数を用いることとする (Kawano and Konishi, 2007)．


(t − τi+2 )2φi (t) = exp −
,
2h2

h=

τi+2 − τi.
3(4)

ここで，τi は昇順の数列で，τ4 ，τm+1 がそれぞれ全観測時点の始点，終点となるように等間隔に配置された節点とする．したがって，(2) 式と (4) 式よりガウス型基底関数に基づ
く非線形回帰モデルは次式で表される．

f (xαj |tαj ; cα , σx2α ) = 12πσx2α{xαj − cTα φ(tαj )}2
exp −
,
2σx2α

j = 1, . . . , Nα .

(5)係数パラメータ cα と誤差分散 σx2α は，正則化法に基づいて推定する．すなわち，モデルの対数尤度関数に正則化項を付与した正則化対数尤度関数

ζα (cα , σx2α ) =

Nα

j=1

log f (xαj |tαj ; cα , σx2α ) −2

N α ζα T
c Ωcα
2 α

=−

Nα
1
Nα ζα Tlog(2πσx2α ) − 2 (xα − Φα cα )T (xα − Φα cα ) −c Ωcα (6)
2
2σxα
2 α

の最大化により推定する．ここで，Φα = (φ(tα1 ), . . . , φ(tαNα ))T ，xα = (xα1 , . . . , xαNα )Tであり，ζα (> 0) は罰則の度合いを調整する正則化パラメータ，Ω は m × m の既知の非負値定符号行列とする．特に，本報告では D2 を 2 次の差分行列としたとき，非負値定符号
行列として Ω = D2T D2 を用いる．(6) 式の最大化によって得られる正則化推定量 ĉα , σ̂x2αは，それぞれ次式で与えられる．

ĉα = (ΦTα Φα + Nα ζα σ̂x2α Ω)−1 ΦTα xα ,σ̂x2α =

1(xα − Φα ĉα )T (xα − Φα ĉα ).
Nα(7)

正則化パラメータ ζα および基底関数の個数 m の選択方法については，Konishi and Kitagawa (2008) を参照されたい．図 1 は，ガウス型基底関数を用いた非線形回帰モデルによる離散データの関数化の例を表している．ここで，丸印は離散データ，図下部の山形の曲線群が基底関数，図上部の曲線が得られた関数データである．
したがって，uα (t) の推定量 ûα (t) = ĉTα φ(t) を関数データ xα (t) として，離散データ集
合 {(tαj , xαj ); j = 1, . . . , Nα , α = 1, . . . , n} に替えて，関数データ集合{xα (t); t ∈ T , α = 1, . . . , n}(8)をデータとして用いる．

3

半教師あり関数データ判別

3.1

関数ロジスティックモデル

いま，n1 組のラベルあり関数データ {(xα (t), gα ); α = 1, . . . , n1 } と (n − n1 ) 組のラベルなし関数データ {xα (t); α = n1 + 1, . . . , n} に基づいて L 群半教師あり関数ロジスティック
判別方式を構成する．群の所属が既知の関数データ xα (t) が観測されたとき，それが第 k群からのものである事後確率をPr(gα = k|xα ),k = 1, . . . , L,とすると，関数データ集合に基づくロジスティックモデルはPr(gα = k|xα )log
Pr(gα = L|xα )= βkf +xα (t)βk (t)dt = βkT zα ,k = 1, . . . , L − 1,(9)として与えられる (Araki et al., 2009)．ただし，βk (t) は (3) 式の基底関数を用いること
により，

βk (t) =

m

j=1

3

βkj φj (t)(10)

と基底展開できるものと仮定し，βk = (βkf , βk1 , . . . , βkm )T ，zα = (1, wαT J)T とする．また，J は基底関数交差積行列 J = φ(t)φT (t)dt を表し，各 (i, j) 成分は以下で与えられる．

Jij =

√





(τi+2 − τj+2 )2πh2 exp −
.
4h2

(11)

一方，(9) 式を変形すると

exp{βkT zα },T
1 + L−1j=1 exp{βj zα }Pr(gα = k|xα ) =Pr(gα = L|xα ) = 1 −L−1k = 1, . . . , L − 1,

(12)

1

Pr(gα = k|xα ) =

1+

k=1

,
L−1
T
j=1 exp{βj zα }(13)Tが得られる．確率 Pr(gα = k|xα ) はパラメータ β = (β1T , . . . , βL−1
)T に依存しているため，Pr(gα = k|xα ) = πk (xα ; β) と記述するものとする．
次に，群のラベルを表す (L − 1) 次元ラベル変数 Y = (Y1 , . . . , YL−1 )T を

⎧
⎪
⎨

(k)
(α)
(α) T
(0, . . . , 0, 1 , 0, . . . , 0)Tyα = (y1 , . . . , yL−1 ) =⎪T⎩(k = 1, . . . , L − 1),

if gα = k,(0, . . . , 0)

if gα = L,と定義する．このとき，ラベル変数ベクトル yα は多項分布

f (yα |xα ; β) =

L−1


(α)

πk (xα ; β)yk {πL (xα ; β)}1−L−1 (α)
y
j=1 j

(14)

k=1

に従う．また，(n − n1 ) 個のラベルなし関数データ {xα (t); α = n1 + 1, . . . , n} に対して，どの群から得られたかを表す観測変数
⎧
⎪
⎨

(k)
(α)
(α) T
(0, . . . , 0, 1 , 0, . . . , 0)Ttα = (t1 , . . . , tL−1 ) =⎪T⎩

if xα (t) が第 k 群に属する,
if xα (t) が第 L 群に属する,(0, . . . , 0)を別に得たとする．このとき，tα は (14) 式と同じ確率 πk (xα ; β) を持つ多項分布に従うと仮定することができる．したがって，n1 個のラベルあり関数データと (n − n1 ) 個のラ
ベルなし関数データに基づいて，対数尤度関数を

(β) =

n1 L−1


α=1

+

k=1


(α)
yk βkT zα − logn
α=n1 +1

L−1
 (α)
k=1

tk

1+

L−1

l=1



exp{βlT zα }

βkT zα − log

1+

L−1

l=1



exp{βlT zα }(15)として構成することができる．
また，将来のデータ xα (t) は事後確率 Pr(gα = k|xα ) (k = 1, . . . , L − 1) が最大となる群に判別される．

4

3.2

モデルの推定T
モデルのパラメータ β = (β1T , . . . , βL−1
)T は，対数尤度関数に正則化項を付与した次の正則化対数尤度関数の最大化に基づいて推定する．

λ (β) = (β) −n1 λ L−1
β T Kβk .
2 k=1 k(16)

ここで，λ (> 0) は正則化パラメータであり，モデルの複雑さを調整する働きを持つとともに，モデルの安定化に寄与する機能も持つ．K は (m + 1) × (m + 1) 型の行列であり，具体的には，
⎛
⎞
T
0
0
⎠
K=⎝
(17)
0 K∗
で与えられる．ここで，0 は m 次元 0 ベクトル，K ∗ は m × m 型の既知の非負値定符号行列である．K ∗ としては，単位行列 Im ，2 次差分行列 D2 による D2T D2 などさまざまものが用いられるが (詳細については，例えば，Konishi and Kitagawa, 2008 参照)，本報告では単位行列 Im を用いる．
(α)
しかし，(16) 式の正則化対数尤度関数は未知の量 tk を含んでいるため，この正則化対(α)
数尤度関数の最大化は一般に困難である．そこで，変数 tk を潜在変数と捉え，EM アル
ゴリズム (Dempster et al., 1977) を適用することを考える．
EM アルゴリズムは不完全な状態で観測されたデータについて，最尤法 (正則化法) に基づいた推測を行うための推定アルゴリズムであり，統計学のみならず工学・医学・経済学などのさまざまな分野で用いられている (詳細については，例えば，McLachlan andPeel, 2000; 小西他, 2008 参照)．EM アルゴリズムを今回の半教師あり関数ロジスティッ
クモデルに適用することにより，(16) 式の最大値を与える推定値 β̂ は以下のアルゴリズ
ムに従い求めることができる．

[推定アルゴリズム]
Step1 ラベルあり関数データ {(xα (t), gα ); α = 1, . . . , n1 } のみに基づいて，パラメータβ を推定する．また，この関数の最大化には数値的最適化法の 1 つである NewtonRaphson 法を用いる (詳細については，例えば，Green and Silverman, 1994 参照)．
Step2 得られた推定値 β̂ に基づいて L 群判別方式 P̂r(gα = k|xα ) = πk (xα ; β̂) を構成する．
Step3 Step2 で得た判別方式を用いることにより，ラベルなし関数データ xα (t) (α = n1 +1, . . . , n) に対して事後確率 P̂r(gα = k|xα ) (k = 1, . . . , L) を求める．この事後確率に従い，tα を各 α に対して，次のように推定する．
(α)

(α)

t̂α = (t̂1 , . . . , t̂L−1 )T = (P̂r(gα = 1|xα ), . . . , P̂r(gα = L − 1|xα ))T .
(α)

(α)

(18)

Step4 (16) 式に含まれる tk を t̂k に置き換えた関数の最大化に基づいて推定値 β̂ を得る．また，この関数の最大化には Newton-Raphson 法を用いる．
5

Step5 収束条件が満たされるまで Step2 から Step4 を繰り返す．
ここで，上記のアルゴリズムの収束条件は，j 番目，(j + 1) 番目の更新における正則化対
(j)
(j+1)
(j+1)
(j)
数尤度関数をそれぞれ λ ，λと定めるとき，|λ− λ | < ε とする (ε は十分小さ
い正の数とする)．
(16) 式の最大化に基づいて得られる推定値 β̂ を (14) 式に代入することにより，統計モ
デル f (yα |xα ; β̂) が得られる．しかし，この統計モデルは正則化パラメータ λ の値に依存しているため，このパラメータの値を適切に選択する必要がある．次節では，このパラメータの値を適切に選択するために，統計モデル f (yα |xα ; β̂) を評価する基準をベイズ理論の観点から導出する．

4

モデル評価基準

Schwarz (1978) は，ベイズアプローチに基づいてモデル評価基準 BIC を提案した．BIC
は，最尤法によって推定されたモデルを評価するための基準であったが，Konishi et al.
(2004) はパラメータ β の事前分布を (退化した) 正規分布と考えることで，正則化法により推定されたモデルの事後確率のラプラス近似から得られるモデル評価基準を導出した．
GBIC = −2

n1

α=1

logf (yα |xα ; β̂) + n1 λ

L−1

k=1

β̂kT K β̂k + log|R(β̂)|




2π
.
− (L − 1)log|K|+ − (L − 1)(m + 1 − d)logλ − (L − 1)d log
n1

(19)

ここで，|K|+ は行列 K の 0 でない固有値の積であり，d は行列 K の階数，R(β̂) は次で与えられる (L − 1)(m + 1) × (L − 1)(m + 1) 型行列である．

R(β̂) = −1
1
(H  G)T (H  G) + I + λJ.
n1
n1

(20)

ここで， は行列の成分毎の積を表す．G は G = (Z, . . . , Z) となる n1 × (L − 1)(m + 1)型行列であり (Z = (z1 , . . . , zn1 )T )，H = (π1 1Tm+1 , . . . , πL−1 1Tm+1 ) (πk = (Pr(g1 =k|x1 ), . . . , Pr(gn1 = k|xn1 )))，I = diag{Z T diag(π1 )Z, . . . , Z T diag(πL−1 )Z}，J = diag{K,. . . , K} である．
モデル評価基準 GBIC の値を最小にする正則化パラメータ λ の値を選択し，対応するモデルを最適なモデルとする．

5

適用例

Spellman et al. (1998) による遺伝子発現プロファイルデータへの適用を通して提案手法の有効性を検証する．このデータは，77 時点における 6178 個の酵母菌の遺伝子発現量を測定したものであり，今回は 24 時点における細胞周期が同定された欠損値を除く 632

6

−3

−2

1

2

5

0

−1.5

−1.0

−1

10
15

5

7
0.5

1.0

2

1.5

15

0.0

1

10

−0.5

expression level0

expression level

5−1

expression level

−2

−1.5

−2

−1.0

−1

−0.5

0.0

0.5

expression level0

expression level
1

1.0

1.5

2

M/G1

G1

time
20

5

time
20
5

10

15

10

time

S
S/G2

10

15
20

time
15
20

G2/M

time

20

図 2: 各群における関数データ集合の例．各群にはそれぞれ 10 本の関数データが描かれ
ている．

0.35
0.30
0.25

Prediction errors

15

20

25

30

35

40

45

50

Percentages of classified data sets

図 3: 10 回の反復に対する予測誤差の平均の推移 (実線：提案手法，破線：教師あり関数
ロジスティック判別手法)
個の遺伝子データを用いた．また，細胞周期は M/G1，G1，S，S/G2，G2/M の 5 群に分類されている．
まず，各遺伝子に対して，得られた 24 時点の離散データに非線形回帰モデルをあてはめることにより，図 2 のような関数データ集合を得た．次に，得られた 632 個の関数デー
タを無作為に 300 個の学習データと 332 個のテストデータに分けた．また，半教師あり学習を行うために，学習データ数の 15%, 20%, 30%, 40%, 50%をラベルあり関数データと無作為に割り当て，各ラベルあり関数データの割合に対して残りをラベルなし関数データと設定した．
適用に際して，得られた関数データに基づいて半教師あり関数ロジスティックモデルを構築し，332 個のテストデータを用いて予測誤差を計算した．また，正則化パラメータ λの値は 4 節で導出したモデル評価基準 GBIC の最小化に基づいて選択した．比較対象として，Araki et al. (2009) により提案された教師あり関数ロジスティック判別モデルを用いた．ここで，上記の教師あり関数ロジスティックモデルは，ラベルあり関数データのみに基づいて推定されることに注意されたい．
図 3 は 10 回の反復により得られた予測誤差の平均の推移を表している．ここで，実線が提案手法，破線が教師あり関数ロジスティック判別手法 (Araki et al., 2009) を示す．図3 から，提案手法は，学習データ中のラベルあり関数データの割合が少ないときに予測誤差を小さくしているのが見て取れる．これより，少量のラベルあり関数データと大量のラベルなし関数データが得られたという状況においては，提案手法は従来の教師あり判別手法より有効に働くものと考えられる．

8

参考文献
[1] Araki, Y., Konishi, S., Kawano, S. and Matsui, H. (2009): Functional logistic discrimination via regularized basis expansions. Communications in Statistics - Theory
and Methods (to appear).
[2] Chapelle, O., Schölkopf, B. and Zien, A. (2006): Semi-Supervised Learning. Cambridge, MA: MIT Press.
[3] Chapelle, O. and Zien, A. (2005): Semi-supervised classiﬁcation by low density separation. Proceedings of the 10th International Workshop on Artiﬁcial Intelligence andStatistics, 57–64.
[4] Dean, N., Murphy, T.B. and Downey, G. (2006): Using unlabelled data to updateclassiﬁcation rules with applications in food authenticity studies. Journal of the Royal
Statistical Society C, 55, 1–14.
[5] Dempster, A.P., Laird, N.M. and Rubin, D.B. (1977): Maximum likelihood fromincomplete data via the EM algorithm (with discussion). Journal of the Royal Statistical Society B, 39, 1–38.
[6] Fujii, T. and Konishi, S. (2006): Nonlinear regression modeling via regularizedwavelets and smoothing parameter selection. Journal of Multivariate Analysis, 97,
2023–2033.
[7] Green, P.J. and Silverman, B.W. (1994): Nonparametric Regression and Generalized
Linear Models. London: Chapman & Hall.
[8] Imoto, S. and Konishi, S. (2003): Selection of smoothing parameters in B-spline
nonparametric regression models using information criteria. Annals of the Instituteof Statistical Mathematics, 55, 671–687.
[9] Kawano, S. and Konishi, S. (2007): Nonlinear regression modeling via regularized
Gaussian basis functions. Bulletin of Informatics and Cybernetics, 39, 83–96.
[10] Konishi, S., Ando, T. and Imoto, S. (2004): Bayesian information criteria andsmoothing parameter selection in radial basis function networks. Biometrika, 91,
27–43.
[11] Konishi, S. and Kitagawa, G. (2008): Information Criteria and Statistical Modeling.
New York: Springer.
[12] 小西貞則, 越智義道, 大森裕浩 (2008): 計算統計学の方法. 朝倉書店.
[13] McLachlan, G.J. and Peel, D. (2000): Finite Mixture Models. New York: Wiley.

9[14] Nigam, K., McCallum, A., Thrun, S. and Mitchell, T. (2000): Text classiﬁcationfrom labeled and unlabeled documents using EM. Machine Learning, 39, 103–134.
[15] Pan, W., Shen, X., Jiang, A. and Hebbel, R.P. (2006): Semi-supervised learning
via penalized mixture model with application to microarray sample classiﬁcation.
Bioinformatics, 22, 2388–2395.
[16] Ramsay, J.O. and Silverman, B.W. (2002): Applied Functional Data Analysis. New
York: Springer.
[17] Ramsay, J.O. and Silverman, B.W. (2005): Functional Data Analysis. Second Edition. New York: Springer.
[18] Schwarz, G. (1978): Estimating the dimension of a model. Annals of Statistics, 6,461–464.
[19] Spellman, P. T., Sherlock, G., Zhang, M. Q., Iyer, V. R., Anders, K., Eisen, M.
B., Brown, P. O., Bostein, D., Futcher, B. (1998): Comprehensive identiﬁcationof cell cycle-regulated genes of the yeast Saccharomyces cerevisiae by microarrayhybridization. Molecular Biology of the Cell, 9, 3273–3297.

10