GAMをもう少し理解したい④
前回の続きです。過去記事はこちらから。
GAMの実装
bakfit
さて bakfit ですが、このサブルーチンはバックフィッティングというアルゴリズムを定義しています。このアルゴリズムは、簡単に言えば「関心のある平滑化対象変数を順番に残差に対してフィッティングしていく」方法です。Wikipediaのバックフィッティングのページからアルゴリズムの部分について抜粋しましょう：
ポイントは４行目の (a) のところで、
とあるように、目的変数から定数部（線形部分）および関心のある平滑化対象変数"以外"（非線形部分）を引いています。すなわち「平滑化対象変数で説明される部分」と「確率的な誤差」で構成される部分のみを残し、平滑化を行なっていることがわかります。これを非線形としたい変数全てについて収束するまで順次繰り返す、というのがバックフィッティングの流れです。
ちなみにこのバックフィッティング、Wikipediaのページでは以下のような説明があります：
the backfitting algorithm is equivalent to the Gauss–Seidel method algorithm for solving a certain linear system of equations.
さらにGauss–Seidel法について調べると（こちら）、これはn元の連立一次方程式を反復的に解くための手法であることが解説されているのですが、その中で以下の記述があります：
The Gauss–Seidel method now solves the left hand side of this expression for x, using previous value for x on the right hand side.
で、これは以下のように解くことができます：
ここで括弧の中に着目すると、$b$から$ax$が順次引かれている（二項目と三項目）のですが、総和記号のインデックスには$i$が含まれていません。つまり更新の対象である$i$番目の$x$以外を順番に引いており、バックフィッティングと同じような計算になっていることがわかります。なお括弧の中の三項目は添え字が (k) ですが、二項目では (k+1) になっており、更新後の$x$を直接使っていることがわかります。さらにちなみに、ここで逐次更新するのではなく一度全ての解を求めてから一斉に更新する方法をヤコビ法といいます。
では上記のアルゴリズムが bakfit でどのように実装されているのかを確認してみましょう。と言いつつ実は bakfit 自体は基本的に以下のように backf1 を呼び出しているだけです。
subroutine bakfit(x,npetc,y,w,which,spar,dof,match,nef, etal,s,eta,beta,var,tol, qr,qraux,qpivot,effect,work) implicit double precision(a-h,o-z) logical ifvar integer npetc(7),iter integer n,p,q,which(*),match(*),nef(*),nit,maxit,qrank,qpivot(*) double precision x(*),y(*),w(*),spar(*),dof(*), etal(*),s(*),eta(*),beta(*),var(*),tol, qr(*),qraux(*),effect(*),work(*) n=npetc(1) p=npetc(2) q=npetc(3) ifvar=.false. if(npetc(4)==1)ifvar=.true. maxit=npetc(6) qrank=npetc(7) do i=1,q{work(i)=dof(i)} ### ここで backf1 を呼び出している call backf1(x,n,p,y,w,q,which,spar,dof,match,nef, etal,s,eta,beta,var,ifvar,tol,nit,maxit, qr,qraux,qrank,qpivot,effect,work(q+1),work(q+n+1), work(q+2*n+1),work(q+3*n+1),work(q+4*n+1)) npetc(7)=qrank return end
なので backf1 を確認することにしましょう。
backf1
前回同様まずは backf1 の処理全体をさっと眺めて見通しを良くします。 backf1 ではざっくりと
- 前処理
- バックフィッティング
etaを再定義
という処理が行われます。
subroutine backf1(x,n,p,y,w,q,which,spar,dof,match,nef, etal,s,eta,beta,var,ifvar,tol,nit,maxit, qr,qraux,qrank,qpivot,effect,z,old,sqwt,sqwti,work) implicit double precision(a-h,o-z) logical ifvar integer n,p,q,which(q),match(n,q),nef(q),nit,maxit,qrank,qpivot(p) double precision x(n,p),y(n),w(n),spar(q),dof(q), etal(n),s(n,q),eta(n),beta(p),var(n,q),tol, qr(n,p),qraux(p),effect(n),work(*) double precision z(*),old(*),dwrss,ratio double precision sqwt(n),sqwti(n) logical anyzwt double precision deltaf, normf,onedm7 integer job,info onedm7=1d-7 job=1101;info=1 if(q==0)maxit=1 ratio=1d0 ## 1. 前処理 ### データに重みをつける anyzwt=.false. do i=1,n{ if(w(i)>0d0){ sqwt(i)=dsqrt(w(i)) sqwti(i)=1d0/sqwt(i) } else{ sqwt(i)=0d0 sqwti(i)=0d0 anyzwt=.true. } } ### QR分解する if(qrank==0){ do i=1,n{ do j=1,p{ qr(i,j)=x(i,j)*sqwt(i) } } do j=1,p{qpivot(j)=j} call dqrdca(qr,n,n,p,qraux,qpivot,work,qrank,onedm7) } ### eta に 非線形項 s(i,j) を加算する do i=1,n{ eta(i)=0d0 for(j=1;j<=q;j=j+1){ eta(i)=eta(i)+s(i,j) } } ## 2. バックフィッティング nit=0 while ((ratio > tol )&(nit < maxit)){ deltaf=0d0 nit=nit+1 ### まずは Linear part + error 部分に最小二乗法を当てはめる do i=1,n{ z(i)=(y(i)-eta(i))*sqwt(i) old(i)=etal(i) } call dqrsl(qr,n,n,qrank,qraux,z,work(1),effect(1),beta, work(1),etal,job,info) ### 重みの逆数で戻し、 etal を更新 do i=1,n{ etal(i)=etal(i)*sqwti(i) } ### 次に平滑化対象の変数についてループ for(k=1;k<=q;k=k+1){ j=which(k) do i=1,n{ old(i)=s(i,k) z(i)=y(i)-etal(i)-eta(i)+old(i) } if(nit>1){dof(k)=0d0} ### splsmで平滑化 call splsm(x(1,j),z,w,n,match(1,k),nef(k),spar(k), dof(k),s(1,k),s0,var(1,k),ifvar,work) do i=1,n{ eta(i)=eta(i)+s(i,k)-old(i) etal(i)=etal(i)+s0 } deltaf=deltaf+dwrss(n,old,s(1,k),w) } ### whileループの打ち切り判定 normf=0d0 do i=1,n{ normf=normf+w(i)*eta(i)*eta(i) } if(normf>0d0){ ratio=dsqrt(deltaf/normf) } else {ratio = 0d0} } ### 3. eta を再定義 do j=1,p {work(j)=beta(j)} do j=1,p {beta(qpivot(j))=work(j)} if(anyzwt){ do i=1,n { if(w(i) <= 0d0){ etal(i)=0d0 do j=1,p{ etal(i)=etal(i)+beta(j)*x(i,j) } } } } do i=1,n eta(i)=eta(i)+etal(i) do j=1,q { call unpck(n,nef(j),match(1,j),var(1,j),old) do i=1,n {var(i,j)=old(i)} } return end
３つ目の工程で eta の再定義とあるのですが、これまでの処理（〜 s.wam）において eta は「 y をリンク関数で変換したもの」として扱われてきました。一方 backf1 の中では 「平滑化対象変数によって構成される部分」（非線形/平滑化部分）として扱われており、混乱しやすいので注意が必要です。つまり backf1 の中では y を「線形部分」と「非線形（平滑化）部分」に分割し、前者を etal 、後者を eta として計算が行われます。
では一つずつ見ていきましょう。
1. 前処理
まずは前処理です。ここでは
という処理が行われます。
### weight の平方根をとって逆数にする anyzwt=.false. do i=1,n{ if(w(i)>0d0){ sqwt(i)=dsqrt(w(i)) ## Fortranの組み込み関数。sqrt() sqwti(i)=1d0/sqwt(i) ## weight の逆数 } else{ sqwt(i)=0d0 sqwti(i)=0d0 anyzwt=.true. } } ### qrank が 0 のときは、x に重みを加えた上で QR 分解する if(qrank==0){ do i=1,n{ do j=1,p{ ### 重みを乗じた x を生成 qr(i,j)=x(i,j)*sqwt(i) # qr は n*p の行列 } } do j=1,p{qpivot(j)=j} ### qvivot は長さ p の integer。1~p の数字を格納 ### QR分解。dqrdca のソースは以下 ### https://github.com/cran/gam/blob/master/inst/ratfor/linear.r call dqrdca(qr,n,n,p,qraux,qpivot,work,qrank,onedm7) } ### eta に 非線形項 s(i,j) を加算する ### ただしここで eta は fitted value ではなく 0 スタート ### s は線形部分とは独立な平滑化部分で、バックフィッティングで更新する do i=1,n{ eta(i)=0d0 # 0.0 にリセット for(j=1;j<=q;j=j+1){ # 平滑化対象の変数について eta(i)=eta(i)+s(i,j) # 平滑化部分を順次加算していく } }
上のブロックでQR分解を行なっていますが、これは次の工程で線形部分に対して通常の回帰を実行するためです。以前の記事で glm がどのように実装されているかを調べた際、 glm では、
dqrdc2というサブルーチンを用いてQR分解dqrslというサブルーチンで最小二乗法を当てはめ
という処理で解が得られることを示しましたが、 glm 同様に gam においても線形部分のフィッティングではQR分解によって最小二乗解を求める、という流れになっているようです。
2. バックフィッティング
では続いてバックフィッティングの実装を見ていきましょう。ここが、ある意味で gam の本体とも言えるルーチンになります。
### while ループ開始 nit=0 while ((ratio > tol )&(nit < maxit)){ # ratio のデフォルトは 1.0 。 tol は多分 0.0005 # first the linear fit deltaf=0d0 nit=nit+1 # イテレータに +1 ### まずは Linear part + error 部分に最小二乗法を当てはめる ### y から非線形項を減じたもの（Linear part + error）に重みを乗じて z を更新(これを y として dqrsl に渡す) # etal を old に格納。ループ１回目の時点では etal は Null？ do i=1,n{ z(i)=(y(i)-eta(i))*sqwt(i) old(i)=etal(i) } ### https://github.com/wch/r-source/blob/trunk/src/appl/dqrsl.f ### dqrsl の6個目の引数として z を渡しているが、これは y になる ### etal は11番目の引数として渡しているが、これは xb 。 call dqrsl(qr,n,n,qrank,qraux,z,work(1),effect(1),beta, work(1),etal,job,info) ### 重みの逆数で戻し、 etal を更新 ### etal は y の予測値の線形部分 do i=1,n{ etal(i)=etal(i)*sqwti(i) }
まずは上のブロックで線形部分のフィッティングが行われます。 dqrsl の直前の処理に注目すると、 y から eta を減じたものを z として定義し、これを dqrsl の引数として渡しています。すなわち目的変数から eta （非線形部分）を減じた、線形部分ですね。これを、先ほど解説したように dqrsl というサブルーチンに渡し、最小二乗法による解を得ます。引数には eta ではなく etal が入っていることに注意しましょう。この etal は直後の処理で重みを乗じた上で後続の処理に引き渡されます。
続いて平滑化対象変数の処理に入ります。
### ここで平滑化対象の変数についてループ for(k=1;k<=q;k=k+1){ # 平滑化対象の変数ごとに j=which(k) do i=1,n{ ### old に k 番目の平滑化対象変数の値を入れる old(i)=s(i,k) ### y から etal (Linear part)と eta (Nonlinear part) を引き、old (k 番目の Nonlinear part) を足している ### 残差 + k 番目の平滑化変数の値になっている z(i)=y(i)-etal(i)-eta(i)+old(i) } ### df は 0 にリセットされてしまう if(nit>1){dof(k)=0d0} ### splsm を呼びだす ### 切片と k 番目の平滑化変数を x として渡す ### 9番目の引数 s が splsm 内での平滑化変数 smo なので s が更新される call splsm(x(1,j),z,w,n,match(1,k),nef(k),spar(k), dof(k),s(1,k),s0,var(1,k),ifvar,work) do i=1,n{ ### eta を更新 ### 古い eta にk番目の非線形部分を加算して、古い非線形部分を減じる) eta(i)=eta(i)+s(i,k)-old(i) ### etal を更新(古い etal に s0 を乗じる。) ### s0 は weighted mean of y etal(i)=etal(i)+s0 } ### deltaf を更新（while ループの打ち切り判定） deltaf=deltaf+dwrss(n,old,s(1,k),w) } # ここまで for ループ
細かな処理はコメントとして書いてありますが、このブロックでは平滑化の対象となる変数について平滑化を実行して s を更新します。ポイントは
z(i)=y(i)-etal(i)-eta(i)+old(i)
のところで、この記事の冒頭でバックフィッティングのアルゴリズムを解説した際に「関心のある平滑化対象変数を順番に残差に対してフィッティングしていく」と述べましたが、そのために残差を求めているのが上の式になります。 etal は線形部分、 eta は非線形部分なので両者を引くと確率的な誤差しか残りませんが、そこに old(i) つまり s(i, k) を加算することで、k番目の平滑化対象変数 + 確率的な誤差のみが残ります。
続いてこの誤差に対して平滑化を行なっているのが
splsm(x(1,j),z,w,n,match(1,k),nef(k),spar(k), dof(k),s(1,k),s0,var(1,k),ifvar,work)
の部分です。この splsm は内部でさらに sbart というCのサブルーチンを呼び出しているのですが、この解説を始めるとキリがないのでまた別の機会に取っておきましょう。
平滑化により s が更新されたら、 eta にもk番目の平滑化対象変数の更新を反映させ、最後に while ループの打ち切り判定のための指標を計算します。
### normf の更新（while ループの打ち切り判定） normf=0d0 do i=1,n{ normf=normf+w(i)*eta(i)*eta(i) } ### ratio（while ループの判定に使われる）を計算 if(normf>0d0){ ratio=dsqrt(deltaf/normf) } else {ratio = 0d0} # call DBLEPR("ratio",-1,ratio,1) }
3. eta の再定義
以上の工程で十分に収束したと判断されたらバックフィッティングアルゴリズムは終了し、 s.wam に値を返す処理に入ります。冒頭で解説したように backf1 内では eta は非線形部分として扱われていたので、これを線形部分と加算して元の eta に戻します。なお beta は線形部分の各変数の回帰係数です。
### etal を再定義 do j=1,p {work(j)=beta(j)} do j=1,p {beta(qpivot(j))=work(j)} if(anyzwt){ do i=1,n { if(w(i) <= 0d0){ etal(i)=0d0 do j=1,p{ etal(i)=etal(i)+beta(j)*x(i,j) } } } } ### 非線形項 + 線形項として eta を再定義 do i=1,n eta(i)=eta(i)+etal(i) ### 結果をまとめる do j=1,q { call unpck(n,nef(j),match(1,j),var(1,j),old) do i=1,n {var(i,j)=old(i)} } return end
さて、これで backf1 および bakfit の処理は終了したので s.wam に戻ることになります。 s.wam の残りの工程では、前回記事で確認したように結果を格納したオブジェクトである rl を生成するのでした。それと元々のデータを順位で置き換える形で生成された smooth.frame をまとめて返します。
さらに s.wam を抜けると gam.fit 内でのループに戻ってくるわけですが、これは前々回記事で解説した通り、残差からデビアンスを計算した上でループの打ち切りを判定します。
gam.fit
収束が十分であれば、 gam.fit のループを抜けます。残りの工程は基本的に結果をまとめあげて gam に返すだけのようです。
### gam.fit に戻ってきた fitqr <- fit$qr xxnames <- xnames[fitqr$pivot] nr <- min(sum(good), nvars) if (nr < nvars) { Rmat <- diag(nvars) Rmat[1:nr, 1:nvars] <- fitqr$qr[1:nr, 1:nvars] } else Rmat <- fitqr$qr[1:nvars, 1:nvars] Rmat <- as.matrix(Rmat) Rmat[row(Rmat) > col(Rmat)] <- 0 dimnames(Rmat) <- list(xxnames, xxnames) names(fit$residuals) <- ynames names(mu) <- ynames names(eta) <- ynames ### eta and mu fit$additive.predictors <- eta fit$fitted.values <- mu names(fit$weights) <- ynames names(fit$effects) <- c(xxnames[seq(len = fitqr$rank)], rep.int("", sum(good) - fitqr$rank)) if (length(fit$smooth) > 0) fit$smooth.frame <- smooth.frame[smooth.labels] wtdmu <- if (a$intercept) sum(weights * y)/sum(weights) else linkinv(offset) nulldev <- sum(dev.resids(y, wtdmu, weights)) n.ok <- nobs - sum(weights == 0) nulldf <- n.ok - as.integer(a$intercept) rank <- n.ok - fit$df.residual aic.model <- aic(y, nobs, mu, weights, new.dev) + 2 * rank if (!is.null(fit$smooth)) { nonzeroWt <- (wz > 0) nl.chisq <- gam.nlchisq(fit$qr, fit$residuals, wz, fit$smooth) } else nl.chisq <- NULL fit <- c(fit, list(R = Rmat, rank = fitqr$rank, family = family, deviance = new.dev, aic = aic.model, null.deviance = nulldev, iter = iter, prior.weights = weights, y = y, df.null = nulldf, nl.chisq = nl.chisq)) fit }
gam
ようやく gam に帰ってきました。 gam.fit を呼び出した以降の処理は以下のブロックだけです。特定の条件下においてNull Deviance（切片のみ考慮されたモデルによる逸脱度）を計算し、結果をまとめあげて返します。
### offset が指定されており intercept 項がある場合 if (length(offset) && attr(mt, "intercept") > 0) { fit$null.dev <- glm.fit(x = X[, "(Intercept)", drop = FALSE], y = Y, weights = weights, offset = offset, family = family, control = control[c("epsilon", "maxit", "trace")], intercept = TRUE)$deviance } if (model) fit$model <- mf fit$na.action <- attr(mf, "na.action") if (x) fit$x <- X if (!y) fit$y <- NULL fit <- c(fit, list(call = call, formula = formula, terms = mt, data = data, offset = offset, control = control, method = method, contrasts = attr(X, "contrasts"), xlevels = .getXlevels(mt, mf))) class(fit) <- c("Gam", "glm", "lm") if (!is.null(fit$df.residual) && !(fit$df.residual > 0)) warning("Residual degrees of freedom are negative or zero. This occurs when the sum of the parametric and nonparametric degrees of freedom exceeds the number of observations. The model is probably too complex for the amount of data available.") fit }
終わりに
以上で終わりです！長くなりましたが、ようやく gam の内部で何が行われているかを確認することができました。振り返ると、 gam という関数のコアとなる部分では：
- バックフィッティングによる非線形項のフィッティング（
bakfitおよびbackf1） - バックフィッティングに渡すためにデータを順位に置き換える（
s.wam） {general, s, lo}.wamのいずれのエンジンに渡すかを特定し、bf.callを作成（gam.fit）gam.fitに渡すための諸々の設定（gam）
といった処理が行われているようでした。またバックフィッティングというアルゴリズムは、端的には「関心のある変数による影響のみを残してフィッティングすることを繰り返す」というものであることがわかりました。
一方で、今回の一連の記事では平滑化の仕組みを追いかけることがまだ出来ていません。また平滑化の関数にしても s.wam しか見ておらず、 lo.wam や general.wam の挙動を確認していません。その辺りもいずれ調べたいですね。
樋口先生退任記念シンポジウム参加記録
11/5に開催された元統計数理研究所の所長である樋口先生の退任記念シンポジウムに参加してきましたので、そのメモを共有しておきます。ご参考まで。やっぱり長らく研究をされてきた方のお話というのは面白いですね。
ご挨拶
統数研 椿広計 氏
バイオインフォマティクスとベイズモデリング
東京大学 井本清哉 氏
- 樋口先生との共同研究
- ベイジアンネットワークによる遺伝子ネットワークの推定
- 遺伝子間の因果関係（制御関係）
- グラフ構造の推定
- ヒトゲノム、21000個の遺伝子
- それらが複雑に関係し合う
- たかだか有限のデータでは推定できない
- ネットワーク構造の事後確率
- 生物学的知識、他のデータで事前確率を構成
- 結合配列
- タンパク質の相互作用
- 進化情報
- データベース
- 遺伝子ノックダウン
- バイオインフォマティクスはベイズモデリングの格好の応用対象であり、問題の宝庫
- マイクロバイオームを含んで「私」
- 健康や病気をデータ化できる時代が来た
- データ同化による細胞内分子シミュレーションモデルの構築
- インフルエンザ感染拡大エージェントベースシミュレーション
- インフルエンザの拡大を抑えるためには会社員のワクチン接種が有効
- ベイジアンネットワークによる遺伝子ネットワークの推定
- 予測だけを目的としたら、AI（機械学習技術）は人を超えたかも
- Kaggle、AutoMLがコンペで2位になった
- ゲノム医療、理由が説明できないと医療に入っていくことは難しい現状
- 医療分野で AI は Augumented Intelligence
- 人の能力を向上させるもの
- 医科学研究所、2900万件の医学論文を学習
- 思考のリミッターを外しなさい by 樋口先生
消費者行動理解のためのベイジアンモデリング
筑波大学 佐藤忠彦 氏
- マーケティングに関わる領域も最近は樋口先生の表芸になってきている
- 能書き
- マーケティング、現実の市場で発生する現象
- 問題を解決するために以下の問いに答える必要がある
- マーケティングを取り囲む環境は？
- 市場はどこにある？
- 敵は誰？
- 顧客はどこに？
- 顧客の嗜好は？
- 消費者の価格反応は？
- 広告、販促への反応は？
- ものやサービスを選ぶ基準は？
- 上記は経験や感覚で評価できるか？
- 科学的アプローチが必要
- 消費者を理解する ＝ 仮説を発見する
- 消費者を理解する ≠ 真の姿を解明する
- 消費者理解のプロセス
- 複数のモデルを立て
- 全てのモデルを推定・比較し
- この時点でもっとも妥当と考えられる構造の推定
- 情報の変換
- 結果の活用
- 現在持ちうる知識で記述能力・汎化能力が高いモデルを選ぶ
- 消費者理解のキーワード
- 消費者異質性
- 時間的異質性
- 潜在変数
- パラメータがいっぱいで普通のモデルでは対応できない
- 事後分布が関心の対象
- 樋口流ベイジアンモデリング
- マーケティングの事例①：佐藤、樋口（2009）
- スーパーで牛乳を購買する確率が、潜在変数である家庭内在庫量/消費量の推移に基づき変化する構造
- 一般状態空間
- 店に来て牛乳を買う状況
- n時点での在庫量 = 前日開始時点での在庫量 + 購買量 - 消費量
- 消費量 = 残っている量 * 消費率
- 店に来る効用、牛乳を買う効用
- 動的に変化する
- 在庫量が変化すると購入確率がどのように変化するか
- マーケティングの事例②：宮津、佐藤（2015）
- まとめ
固体地球科学とデータ同化
- 北川源四郎 時系列解析プログラミング
- データ同化
- 数値シミュレーションと観測データの融合
- 大規模な状態空間モデルと大規模データ
- 四次元変分法
- データ同化入門
樋口先生とブリジストンと私
ブリジストン 花塚泰史 氏
- 樋口先生と10年以上の関係がある
- タイヤ事業が8割、多角化事業が2割
- 売り切り型からソリューションプロバイダーとして
- 状態空間モデルとカーネルマシンを併用した・・・
- ブリジストンでの事例
- タイヤ、ほぼ唯一の電子化されていない部品
- CAISによる安全・安心なモビリティ社会への貢献
- インテリジェント化したタイヤ
- 空気圧、温度を拾うのに加え、タイヤ内のセンサーから
- 磨耗状態推定
- 路面状態推定 → 自動車の制御へフィードバック
- 加速度
- 路面状態に応じて波形が特徴的に変化
- ウェット、凍結
- 周波数領域において特徴的な波形
- Windowをシフトしながら
- カーネルSVM
- （速度によって）系列長が異なるものを比較
- グローバルアライメントカーネル
- 計算量がでかい
- リアルタイムでやりたい
- 窓同士の類似度
- ドライとウェットを分類する判別問題
- GAカーネルの高速化
- 類似度を計算する窓を制限する
- タイヤの波形は特徴的な２つのピークが出るので・・・
- 元のGA、93.9%で判別できるがタイヤ100回転分の時間が必要
- 工夫により96.9％、タイヤ一回転分ぐらいで判断できるようになった
- ７種類に判別
- 冬季道路管理の最適化への適用
- 路面のメンテナンス
- 巡回車に乗っている人の目で判断
- 北海道で実用化されている
- ISCOS
- 路面のメンテナンス
- 次世代のデータサイエンティスト育成
- 社内研修
- 将来的な発展
- センシング技術
- 解析、予測
- デジタルサービス
ベイズモデリングと歩んだ30年
中央大学 樋口知之 氏
- 平成元年に入所、平成31年に退任
- 200MBのデータ解析に大型計算機が必要だった時代
- 人間の未来予測は当たらない
- 1968年に2001年を予測した映画
- 1982年に2019年を予測した映画（ブレードランナー）
- モデルとモデリング
- 何が違うのか
- 常にモデルを改良する姿勢
- ものの見方、捉え方を柔軟に変化
- 思考のプロセス自体を科学する
- モデルを比較するための羅針盤としての情報量基準
- 何が違うのか
- 赤池先生との出会い
- 多変量自己回帰を応用したセメント工場、火力制御
- 船のオートパイロット（北川先生）
- 当時
- 支配方程式がない、計算コストが大
- 今
- データから直接モデリング
- 生成モデルと識別モデル
- 生成モデルはデータが生成される過程
- 識別モデルはデータやパラメータを与えたときにどのように識別されるか
- 後者が簡単（機械学習コミュニティで盛ん）
- 統計の人たちはデータがどんな機序で生まれたのかに強い関心がある
- 頻度主義 vs ベイズ主義
- ベイズの定理がなぜ役立つか？
- コンピュータの性能向上
- 周辺尤度の計算
- 高精度センサーのコモディティ化
- 尤度関数
- ストレージの廉価化
- 細かい情報を保存
- コンピュータの性能向上
- 1980年代の統数研
- 統数研
- アカデミアから産業界への人材流出は欧米でも同じ
- クロスアポイント
- メゾスコピックモデリング
- 数理表現された法則
- 第一原理法則
- 演繹と帰納
- 逐次データ同化技術
- KFとPFの間にある
- 状態ベクトルの次元数、アンサンブルメンバー数
- 統計的推論法のパラダイムシフトを促す外的要因
- ビッグデータ
- データ駆動サービス社会
- 因果から相関へ
- 目的特化型計算機
- 分析者の技術から計算パワーへ
- 力学的ダウンスケーリング
- エミュレーション、仮想計測
- やっていること
- GANを使ってエミュレータを作る
- 深層学習と状態空間モデルの融合
- 高等教育
- 機械学習と統計
- Data driven vs Science driven
- Predictive model vs Descriptive model
- Correlation vs Causality
- 認識科学から設計科学
- 対象理解から帰納の最適化
- 「真のモデルなど存在しない」という見方は、当面の問題のイメージ能力の欠如
AI社会での数理モデリングと現場主義の意義と価値
東京大学 合原一幸 氏 東北大学 早瀬敏幸 氏 大阪大学 鷲尾 隆 氏 ブリジストン 花塚泰史 氏 中央大学 樋口知之 氏
- 分野横断的なアプローチで豊富な共同研究実績の持ち主
- 合原氏
- 数理工学
- 統計とは兄弟みたいな
- 暮らしを変える驚きの数理工学
- ニューラルコンピュータ
- アルファ碁解体新書
- AIと人間の共同作業、協働
- Pair Go
- ファッション
- GANでデザイナーっぽいデザイン
- 数理工学
- 早瀬氏
- 流体力学、制御工学、生体工学
- データ同化
- 鷲尾氏
- 三者は生成モデルにこだわりがある
- ３つのお題（６つの事前質問から）
- 現場主義
- １．現場主義とKKDの違いは？ 現場主義に重要性について
- ２．学生や研究者に現場主義をどう教育していくべきか？
- 樋口
- 現場に出向く研究者魂
- 現場の人と厳しく議論できるコミュニケーション力、分野跳躍力
- 本質的な課題を数理の問題にしたてるセンス
- しつこさ
- 鷲尾
- 向こうの分野の研究者に統計や機械学習を教えた方が良いという経験
- 現場を知っている人の方がよく使える
- 合原
- 中小企業において、社長が現場を見に行く
- そうすると階層構造にショートカットが生まれる
- 中小企業において、社長が現場を見に行く
- 早瀬
- 研究時間が減っていて、研究に割ける時間がない
- 花塚
- 樋口
- コミュニケーションが取りづらい相手は？
- 無茶振りする（合原）
- 専念しても一年かかるようなもの
- データ化されていないものを大量に持ってくる
- 割ける時間の中でうまく成果が出るテーマ選び
- 他の人に無茶振りしてみる
- その人だったら無茶ではないかもしれない
- 人脈づくり、人的ネットワーク
- 有名な人を知っているだけではダメ
- あの人はいま時間がある
- 分析・テーマの振り返りが繰り返されるので、じっくり取り組む必要がある（鷲尾）
- ボキャブラリーが通じないとこからスタート
- 無茶振りする（合原）
- 人とAI
- 教育システム
- ５．分野横断的な学術体系をどのようにして教育システムに組み込むか
- ６．１０年後に国立大学はどのような姿になっているか
- 東北大学では全員にデータサイエンス系の講義を受けさせることになる（早瀬）
- オーストラリアの某大学は数学科を三等分している（合原）
- 若い人を早めに独立させる（合原）
- 教授、准教授、助教みたいな階層構造をやめる
- 数学が苦手な人たちにいかに脱落させないか（鷲尾）
- OJTとして授業より前に使ってもらう
- 古典的製造業は縦割りがピシッとしている（花塚）
- 物売りからこと売り、ソリューションやるなら横串が必要
- サプライチェーンの初めから終わりまでを担う
- それがデータサイエンス組織
- コミュニケーションを密にして話をしにいく
- 現場ごとにデータサイエンスを理解している人が必要
- CoEとしての機能
- 育成
- CoE
- 高度な手法を理解できる
- 現場
- 話を理解できる人を多く育てる
- 裾野を広げる
- CoE
- 知っている人をいかにつなげていくか
- 統数研に期待すること、文科省に頑張ってほしいこと
- 予算（合原）
- 日本のトップの研究所として世界をリードする研究（早瀬）
- 地道に学問を究めるところが魅力（鷲尾）
- 外にアピールして外部資金を獲得
- 基礎研究を地道にコツコツと続けられる
- リカレント教育（花塚）
- DSを外から取ってくるのは至難
- 社内だけでの教育も困難
- アカデミアとの連携
- 現場主義
GAMをもう少し理解したい③
前回の続きです。過去記事はこちらから。
ushi-goroshi.hatenablog.com ushi-goroshi.hatenablog.com
GAMの実装
s.wam
s.wam は gam.fit の中で bf.call として定義されたフィッティングのための関数で、 s.wam の wam はweighted additive modelの意味のようです。名前からして gam の本体のような気がしてきますが、先に言ってしまうとこの関数内で bakfit というFortranの関数を呼び出しており、少なくともRの「 gam というライブラリの本体」という意味では合っているのではないでしょうか。
s.wam の説明に入る前に、見通しを良くするために関数全体をさっと眺めてみましょう。 s.wam はざっくりと以下の４つの処理に分けることができそうです：
- 平滑化の対象変数について元のデータを順位に置き換え、smooth.frameを再定義する
- 後でFortranに渡すために必要な指定を行う
- Fortran で書かれた bakfit を呼び出す（本体）
- 後処理
### gam::s.wam function (x, y, w, s, which, smooth.frame, maxit = 30, tol = 1e-07, trace = FALSE, se = TRUE, ...) { ### 1. 平滑化の対象変数について元のデータを順位に置き換え、smooth.frameを再定義する if (is.data.frame(smooth.frame)) { first <- TRUE data <- smooth.frame[, names(which), drop = FALSE] smooth.frame <- gam.match(data) dx <- as.integer(dim(x)) smooth.frame$n <- dx[1] smooth.frame$p <- dx[2] oldClass(data) <- NULL smooth.frame$spar <- unlist(lapply(data, attr, "spar")) smooth.frame$df <- unlist(lapply(data, attr, "df")) } else first <- FALSE ### 2. 後でFortranに渡すために必要な指定を行う storage.mode(tol) <- "double" storage.mode(maxit) <- "integer" which <- unlist(which) storage.mode(which) <- "integer" storage.mode(y) <- "double" storage.mode(w) <- "double" p <- smooth.frame$p n <- smooth.frame$n for (ich in which) x[, ich] = signif(x[, ich], 6) ### 3. Fortran で書かれた bakfit を呼び出す（本体） fit <- .Fortran("bakfit", x, npetc = as.integer(c(n, p, length(which), se, 0, maxit, 0)), y = y, w = w, which, spar = as.double(smooth.frame$spar), df = as.double(smooth.frame$df), as.integer(smooth.frame$o), as.integer(smooth.frame$nef), etal = double(n), s = s, eta = double(n), beta = double(p), var = s, tol, qr = x, qraux = double(p), qpivot = as.integer(1:p), effects = double(n), double((10 + 2 * 4 + 5) * (max(smooth.frame$nef) + 2) + 15 * n + 15 + length(which)), PACKAGE = "gam") ### 4. 後処理 #### 収束に関する警告 nit <- fit$npetc[5] qrank <- fit$npetc[7] if ((nit == maxit) & maxit > 1) warning(paste("s.wam convergence not obtained in ", maxit, " iterations")) #### ループの２回目以降は序盤の処理をスキップさせる if (first) { smooth.frame$spar <- fit$spar first <- FALSE } #### 最終的な返り値となる rl に情報を追加していくための準備 names(fit$df) <- dimnames(s)[[2]] names(fit$beta) <- labels(x)[[2]] qrx <- structure(list(qr = fit$qr, qraux = fit$qraux, rank = qrank, pivot = fit$qpivot, tol = 1e-07), class = "qr") effects <- fit$effects r1 <- seq(len = qrx$rank) dn <- colnames(x) if (is.null(dn)) dn <- paste("x", 1:p, sep = "") names(effects) <- c(dn[qrx$pivot[r1]], rep.int("", n - qrx$rank)) #### rl を生成 rl <- list(coefficients = fit$beta, residuals = fit$y - fit$eta, fitted.values = fit$eta, effects = effects, weights = w, rank = qrank, assign = attr(x, "assign"), qr = qrx, smooth = fit$s, nl.df = fit$df - 1) rl$df.residual <- n - qrank - sum(rl$nl.df) - sum(fit$w == 0) if (se) rl <- c(rl, list(var = fit$var)) c(list(smooth.frame = smooth.frame), rl) }
以下、一つずつ見ていきましょう。
1. 平滑化の対象変数について元のデータを順位に置き換え、smooth.frameを再定義する
s.wam では始めに smooth.frame に対する処理を行います。具体的には
- 平滑化の対象変数を抽出
- 指定のdigitsで丸める
- ソートした上でユニークなデータ数を得る
- 元データをラベル（数値の大小で並べかえたときの順位）に変換
- 行数・列数を付与
- 平滑化パラメータおよび自由度を付与
という加工を、下記の処理にて実施します。
### smooth.frame はループ初回は data.frame なのでここが評価される if (is.data.frame(smooth.frame)) { first <- TRUE data <- smooth.frame[, names(which), drop = FALSE] # 平滑化対象変数を抽出 smooth.frame <- gam.match(data) # 下で解説 dx <- as.integer(dim(x)) # 行列のサイズ smooth.frame$n <- dx[1] # number of row smooth.frame$p <- dx[2] # number of column oldClass(data) <- NULL ### 各列に対する平滑化パラメータおよび自由度 smooth.frame$spar <- unlist(lapply(data, attr, "spar")) smooth.frame$df <- unlist(lapply(data, attr, "df")) } else first <- FALSE
上のブロックのポイントは gam.match という関数で、 gam.match(data) とすると以下の結果が返ります：
> gam.match(data) $o s(year, 4) s(age, 5) [1,] 4 1 [2,] 2 7 [3,] 1 28 [4,] 1 26 [5,] 3 33 [6,] 6 37 ### 中略 $nef s(year, 4) s(age, 5) 7 61
ちなみに data は元データです：
> head(data) s(year, 4) s(age, 5) 231655 2006 18 86582 2004 24 161300 2003 45 155159 2003 43 11443 2005 50 376662 2008 54
これは gam.match という関数において下記の部分が評価された結果です。
### gam.match から xr <- signif(as.vector(x), 6) # signif 関数は指定の有効桁で丸める。デフォルトは6 sx <- unique(sort(xr)) # それをソートしてユニークにする nef <- as.integer(length(sx)) # ユニークなデータ数 if (nef <= 3) stop("A smoothing variable encountered with 3 or less unique values; at least 4 needed") o <- match(xr, sx, nef + 1) # 元のデータを順位に変更する。 +1 は欠損用。 o[is.na(o)] <- nef + 1
match の部分が少しわかりにくいかもしれませんが、以下のような結果が得られます：
> head(xr) # 元の値 [1] 2006 2004 2003 2003 2005 2008 > match(head(xr), sx, nef + 1) # 順位に変換 [1] 4 2 1 1 3 6
これは xr の値が sx の何番目に位置するかを表わしていますが、参照先の sx は元の値をソートしてユニークにしたものなので、順位に変換していることになります。
また gam.match の結果を smooth.frame に代入されているので、この時点で smooth.frame は数値の大小関係のみを持つことになり、また一連の処理で smooth.frame は data.frame ではなくなるため、ループの２回目以降ではこの処理はスキップされます。
2. 後でFortranに渡すために必要な指定を行う
次のブロックは後ろの工程でFortranの関数に渡すための準備となります。
### 後でFortranに渡すために必要な指定 storage.mode(tol) <- "double" storage.mode(maxit) <- "integer" which <- unlist(which) storage.mode(which) <- "integer" storage.mode(y) <- "double" storage.mode(w) <- "double" p <- smooth.frame$p n <- smooth.frame$n
また以下のコードでは、先程の smooth.frame 同様に有効桁を6桁に丸めます。
### 平滑化対象の変数を6桁に丸める for (ich in which) x[, ich] = signif(x[, ich], 6)
3. Fortran で書かれた bakfit を呼び出す（本体）
以上で準備が整いました。以下が s.wam の本体になります。 .Fortran でFortranで書かれた bakfit というモジュールを呼び出します。
### ここが本体。Fortran で書かれた bakfit を呼び出す fit <- .Fortran("bakfit", x, npetc = as.integer(c(n, p, length(which), se, 0, maxit, 0)), y = y, w = w, which, spar = as.double(smooth.frame$spar), df = as.double(smooth.frame$df), as.integer(smooth.frame$o), as.integer(smooth.frame$nef), etal = double(n), s = s, eta = double(n), beta = double(p), var = s, tol, qr = x, qraux = double(p), qpivot = as.integer(1:p), effects = double(n), double((10 + 2 * 4 + 5) * (max(smooth.frame$nef) + 2) + 15 * n + 15 + length(which)), PACKAGE = "gam")
さてこの bakfit 、ソースコードとしてはおそらくこちら
が正しいのかなぁと思うのですが、 gam をインストールする際に同時にダウンロードされる以下のファイル：
/Users/hogehoge/Library/R/3.6/library/gam/ratfor/backfit.r
の方が説明もあってわかりやすいので、こちらベースに進めたいと思います。なお前回記事でも述べたことですが、 s.wam は第二引数として y を受け取ります：
### s.wam function (x, y, w, s, which, smooth.frame, maxit = 30, tol = 1e-07, trace = FALSE, se = TRUE, ...)
ただしこの s.wam にはYそのものではなく、以下のように z として渡されるのでした：
### 以下は gam.fit から再掲 z <- eta - offset z[good] <- z[good] + (y - mu)[good]/mu.eta.val[good] ### 中略 # eval(bf.call) ← これは下の表現となる s.wam(x, z, wz, fit$smooth, which, fit$smooth.frame, bf.maxit, bf.epsilon, trace)
要するにリンク関数で変換したあとのデータを渡している、ということです。
では bakfit に移りたいのですが、その前にまずは残りの工程を見ておきましょう。
4. 後処理
s.wam の最後の部分です。
nit <- fit$npetc[5] # number of iteration qrank <- fit$npetc[7] # qrank ### nit が maxit と同じであった場合、収束していないことを警告 if ((nit == maxit) & maxit > 1) warning(paste("s.wam convergence not obtained in ", maxit, " iterations")) ### first == T の場合、spar を加えて FALSE にして返す ### ループ２回目以降は評価されない if (first) { smooth.frame$spar <- fit$spar first <- FALSE } names(fit$df) <- dimnames(s)[[2]] # 平滑化対象変数の列名 names(fit$beta) <- labels(x)[[2]] # 説明変数 qrx <- structure(list(qr = fit$qr, qraux = fit$qraux, rank = qrank, pivot = fit$qpivot, tol = 1e-07), class = "qr") # 属性を指定しながらオブジェクトを作成 effects <- fit$effects r1 <- seq(len = qrx$rank) dn <- colnames(x) if (is.null(dn)) dn <- paste("x", 1:p, sep = "") names(effects) <- c(dn[qrx$pivot[r1]], rep.int("", n - qrx$rank)) #### 最終的な返り値である rl を生成 rl <- list(coefficients = fit$beta, residuals = fit$y - fit$eta, fitted.values = fit$eta, effects = effects, weights = w, rank = qrank, assign = attr(x, "assign"), qr = qrx, smooth = fit$s, nl.df = fit$df - 1) rl$df.residual <- n - qrank - sum(rl$nl.df) - sum(fit$w == 0) if (se) rl <- c(rl, list(var = fit$var)) #### return c(list(smooth.frame = smooth.frame), rl)
基本的には s.wam という関数の返り値となる rl を生成するための工程と言えそうです。 bakfit はまた次回。
応用統計学フロンティアセミナー「データサイエンスと応用統計学」参加記録
10/19に開催された応用統計学フロンティアセミナーに参加してきましたので、そのメモを共有しておきます。話を聴きながらのメモなので単語しか書けず意味がわかりにくいところもありますが、ご参考まで。なおセミナーの様子は以下のtogetterでまとめられていますので、そちらも合わせてご覧ください。
- データサイエンスにおける応用統計人材の育成
- 楽天技術研究所におけるデータサイエンスおよび統計の様々な応用事例
- 自動車とデータサイエンス
- 企業ビッグデータから捉える企業活動と未来の活用可能性
- 早稲田大学におけるデータサイエンス人材育成への取り組み
- 学校における「統計」教育の課題 我々は木に縁りて魚を求めてはいないか
- 岩崎先生あいさつ
データサイエンスにおける応用統計人材の育成
横浜市立大学 岩崎学 氏
これまでとこれからの統計的データ解析
- 研究目的の設定
- データ収集法の立案（実験、観察、調査）
- データの収集（モニタリング）
- データの電子化
- データのチェック・マージ
- データの集計とグラフ化（記述統計）
- 統計的推測ないしは予測（推測統計）
- プレゼン・文章化、意思決定（終了もしくは最初に戻る） → 役に立たないと意味がない
いまのデータサイエンス、データありきになっているかもしれない
- 推測統計の部分にフォーカスが当たっている
- 全体の流れを見ることが大事
データサイエンスとは？
横浜市立大学での状況
- 線形代数の学部生の試験成績、文系と理系でほとんど差がない（120点満点で76 vs 72とか）
- 2020年、大学院
- ヘルスデータサイエンス専攻（医学部があるので）
- データサイエンス専攻（M20名、D3名）
楽天技術研究所におけるデータサイエンスおよび統計の様々な応用事例
楽天 森正弥 氏
1.3Bユーザー、48Kの店舗（ビジネスパートナー）
250Mのアイテムについて需要予測
- 価格、在庫の調整
- 非線形回帰。基盤を作ったロシア人が手作業でフィッティング。DLに勝つ。
金融市場予測
- Now-Cast
- Forecastに対してNowcast
- 今の状態を推測する
- Hal Varian
- Bayesian variable selection for nowcasting economic time series
- https://www.nber.org/chapters/c12995
- 独自の景気動向指数
- 会社の株価
- 売上と株価、２ヶ月のラグ
- 相関が0.96
- でも予測精度は芳しくない
- 化粧品会社の株価
- 多段のDoc2Vecでメーカー名のリンキング
- 楽天トラベルのデータ、インバウンド
- Now-Cast
自動車とデータサイエンス
日産自動車 上田哲郎 氏
- Connected Carのエキスパート、ITのエキスパート
- 自動車、AI/ITに関しては新卒の人たちの方が専門性高い
- 今はリアルタイムで車の車種まで特定できる
- BigData、 AI、IoT
- 箱根登って下るとバッテリー回復するので沼津まで行ける
- シミュレーションではなく過去のリーフのデータから、どこまで行けるか示せる
- ProPilot 2.0
- マジでやばい
- ハンドルから手を離して良いと言っている（高速に限る）
- ドライバーを監視している、前方を見ている限り
- Adaptive Cruise Control
- 単眼可視光カメラで前方の車との距離を測る
- ナンバープレートのサイズがわかっているのでカメラの画素数から車間距離がわかる
- AI
- オイルフィルターの真贋判定
- 空気抵抗を計算で求める
- スパコンでも数日かかる
- CADとCdの関係性
- Voxel dataとCdの関係性を学習
- 精度はまだまだ
- R値で0.7（R2のこと？）
- GANで車のデザインを生成
- 画像の生成を三次元に拡張する
- データアナリスト
- 問題のモデル化
- 人材育成
- 自動車業界、サプライヤーが持ってくる
- うまくいってない
- データを使っていて車に興味を持っている人に来てもらいたい
- 自動車がなくなっていく？
- MaaS、CASE
- 痛し痒しと思っている
- 体力のあるうちに全部やっておこう
企業ビッグデータから捉える企業活動と未来の活用可能性
帝国データバンク 中川みゆき 氏
- 帝国データバンクの保有するデータ
- 調査員が現地に足を運んで入手した、ネットでは得られないデータベース
- 100近い項目
- 定量、定性
- 中小企業が99%
- 年間60万件超の調査
- 活用事例
- 守りの活用
- 信用リスク
- 攻めの活用
- 取引関係からつながりを見る
- バリューチェーン全体のマッピング
- データはファクト、データサイエンスとは価値を創り出すこと
- 人間が本来もつ、マシンではできないことに注力
- 「未来の力」はマシンインテリジェンスに対してのリーダーシップ
- 人がやる部分
- 内訳
- 複雑な問題解決
- クリティカルシンキング
- 想像力
- マネジメント
- 人間関係調整
- 情緒的知性
- 決断
- サービスディレクション
- 交渉力
- 認識の柔軟性
- 守りの活用
早稲田大学におけるデータサイエンス人材育成への取り組み
早稲田大学 松嶋敏泰 氏
- 大隈重信と統計学
- 統計院を設立
- 統計センターの前身
- 統計院を設立
- DSとは
- Interdisciplinary Field
- concept to unify
- 日本だとSexiest Jobで広まったので、ビジネスの面が強い
- Fourth paradigm of science
- メタサイエンス
- 知覚と思考の拡大
- 早稲田のフォーカス領域
- 専門性（応用先）
- データサイエンス
- 教育
学校における「統計」教育の課題 我々は木に縁りて魚を求めてはいないか
代々木ゼミナール 西岡康夫 氏
- 初等・中等教育の現状
- サッカー、世界で活躍できるようになった理由は多くの子供がサッカーをやるようになったから
- 棟梁レベルのDSを多く育てるなら裾野を広げないといけない
- 国によってはエリート制を採用しているが、日本のカルチャーに合わないのでは
- 教育改革
- センター試験の廃止
- 混沌の極み
- 大学入学選抜は暗記、再生
- ドーリットル、甚兵衛、小ピピン
- 理系尚もて数学す、況や文系をや
- 統計が入ってきた
- 産業界の要請
- 教員研修や教職課程のカリキュラム改革
- まったく進んでない
- 統計のイロハのイもダメ
- 教員研修や教職課程のカリキュラム改革
- 産業界の要請
- 若者の強い自己否定感
- センター試験の廃止
岩崎先生あいさつ
- Technologyは変わっていくがPrincipleは変わらない
GAMをもう少し理解したい②
前回の続きです。 gam.fit() から。
前回記事はこちら。
GAMの実装
gam.fit()
それでは gam.fit の中身を覗いてみましょう。
### gam.fit は x, y に加えて smooth.frame を受け取る。これは gam で作った mf で、中身は平滑化に関する情報を持った data.frame function (x, y, smooth.frame, weights = rep(1, nobs), start = NULL, etastart = NULL, mustart = NULL, offset = rep(0, nobs), family = gaussian(), control = gam.control())
始めに、受け取ったオプションをもとにデータのサイズ、 control や family に関するパラメータの取り出しを行います。
### 列名 ynames <- if (is.matrix(y)) dimnames(y)[[1]] else names(y) xnames <- dimnames(x)[[2]] ### データの行数、列数 nobs <- NROW(y) nvars <- ncol(x) ## その他の gam.control オプション maxit <- control$maxit bf.maxit <- control$bf.maxit epsilon <- control$epsilon bf.epsilon <- control$bf.epsilon trace <- control$trace ### digits, weigths, offset digits <- -log10(epsilon) + 1 if (is.null(weights)) weights <- rep.int(1, nobs) if (is.null(offset)) offset <- rep.int(0, nobs) ### family に関するパラメータ variance <- family$variance dev.resids <- family$dev.resids aic <- family$aic linkinv <- family$linkinv mu.eta <- family$mu.eta if (!is.function(variance) || !is.function(linkinv)) stop("illegal `family' argument") valideta <- family$valideta if (is.null(valideta)) valideta <- function(eta) TRUE validmu <- family$validmu if (is.null(validmu)) validmu <- function(mu) TRUE eval(family$initialize) if (is.null(mustart)) { eval(family$initialize) } else { mukeep <- mustart eval(family$initialize) mustart <- mukeep } ### eta の初期値 eta <- if (!is.null(etastart)) etastart ### エラーチェック else if (!is.null(start)) if (length(start) != nvars) stop("Length of start should equal ", nvars, " and correspond to initial coefs for ", deparse(xnames)) else { coefold <- start offset + as.vector(if (NCOL(x) == 1) x * start else x %*% start) } else family$linkfun(mustart) ### mu の初期値 mu <- linkinv(eta) if (!(validmu(mu) && valideta(eta))) stop("Can't find valid starting values: please specify some") ### デビアンス（残差平方和） new.dev <- sum(dev.resids(y, mu, weights))
ここまでは関数に渡されたオプションを元に処理を進めていますが、以降のプロセスで gam 特有の処理が行われています。まずは smoothers の取り出しです。
### smoothers が指定されている場合に smoother を抽出する ### 今回のケースでは s が入る a <- attributes(attr(smooth.frame, "terms")) smoothers <- a$specials
オブジェクトの属性を取り出すにあたり、 attr は属性名を指定する必要がありますが、 attributes では全ての属性をリスト形式で取り出します。 a というオブジェクトには、具体的には以下のようなリストが格納されます。
> a $variables list(wage, s(year, 4), s(age, 5), education) $factors s(year, 4) s(age, 5) education wage 0 0 0 s(year, 4) 1 0 0 s(age, 5) 0 1 0 education 0 0 1 $term.labels [1] "s(year, 4)" "s(age, 5)" "education" $specials $specials$s [1] 2 3 $specials$lo NULL $specials$random NULL $order [1] 1 1 1 $intercept [1] 1 $response [1] 1 $class [1] "terms" "formula" $.Environment <environment: R_GlobalEnv> $predvars list(wage, s(year, 4), s(age, 5), education) $dataClasses wage s(year, 4) s(age, 5) education "numeric" "numeric" "numeric" "factor"
a$specials$s には 2 3 という数字が入っていますが、これは smooth.frame の2・3列目が平滑化の対象であるということだと思います。
今回は平滑化が含まれるので以下が処理されます：
if (length(smoothers) > 0) { ### NA ではない要素を取り出す smoothers <- smoothers[sapply(smoothers, length) > 0] ### 以下の処理はちょっとよくわからない for (i in seq(along = smoothers)) { tt <- smoothers[[i]] ff <- apply(a$factors[tt, , drop = FALSE], 2, any) smoothers[[i]] <- if (any(ff)) seq(along = ff)[a$order == 1 & ff] else NULL } }
smoothers に格納されている数値が 1 2 になりました。
> smoothers $s [1] 1 2
次に、以下の処理ではパラメータの推定に使われるエンジンを決めます。平滑化が含まれるケースでは general.wam または {s, lo}.wam が、そうでなければ lm.wfit が選ばれます。二種類以上の平滑化関数（s , lo, random）が指定されている、または s lo 以外の平滑化関数が指定されている場合には general.wam が選ばれるようです。
条件にしたがってエンジンを決めたあとは、後ろの工程で評価するための bf.call を作成しています。
### smoother が指定されている場合、ここが処理される if (length(smoothers) > 0) { gam.wlist = gam.smoothers()$wlist smooth.labels <- a$term.labels[unlist(smoothers)] assignx <- attr(x, "assign") assignx <- assign.list(assignx, a$term.labels) which <- assignx[smooth.labels] ### ２つ以上の smoothers が指定されている場合は general.wam が指定される ### wam は weighted additive model if (length(smoothers) > 1) bf <- "general.wam" ### 今回のケース（平滑化の関数として s を指定）はこちらで、 s.wam が指定される。 else { sbf <- match(names(smoothers), gam.wlist, FALSE) bf <- if (sbf) paste(gam.wlist[sbf], "wam", sep = ".") else "general.wam" } ### bf にオプション部分を文字列で結合 bf.call <- parse(text = paste(bf, "(x, z, wz, fit$smooth, which, fit$smooth.frame,bf.maxit,bf.epsilon, trace)", sep = ""))[[1]] s <- matrix(0, length(y), length(which)) dimnames(s) <- list(names(y), names(which)) fit <- list(smooth = s, smooth.frame = smooth.frame) } ### 平滑化しない場合、通常の lm.wfit に渡す。なお general.wam でも lm.wfit が使われる。 else { bf.call <- expression(lm.wfit(x, z, wz, method = "qr", singular.ok = TRUE)) bf <- "lm.wfit" }
なお上記コードブロックの最後に lm.wfit が使われていますが、LMやGLMの中身がどうなっているのかについては過去記事を参照してください。
bf.call が作られたので、反復しながら解を求めにいきます。ここからが本体となるブロックです。
### デビアンス（ループの打ち切り基準） old.dev <- 10 * new.dev + 10 ### ここから反復に入る for (iter in 1:maxit) { ### weight が 0 のデータは除外する good <- weights > 0 varmu <- variance(mu) if (any(is.na(varmu[good]))) stop("NAs in V(mu)") if (any(varmu[good] == 0)) stop("0s in V(mu)") mu.eta.val <- mu.eta(eta) if (any(is.na(mu.eta.val[good]))) stop("NAs in d(mu)/d(eta)") good <- (weights > 0) & (mu.eta.val != 0) ### z を生成。ただし bf.call で二番目の引数は y なので、この z は目的変数の意味 z <- eta - offset z[good] <- z[good] + (y - mu)[good]/mu.eta.val[good] ### wz を生成。重み。今回のケースでは無視して良い wz <- weights wz[!good] <- 0 wz[good] <- wz[good] * mu.eta.val[good]^2/varmu[good]
上のブロックでは weights を元に分析対象を定義し、 z および wz を生成しています。
ところで、今回のケースでは bf.call は以下のようになるのでした：
> bf.call s.wam(x, z, wz, fit$smooth, which, fit$smooth.frame, bf.maxit, bf.epsilon, trace)
z は第二引数となっていますが、後ほど紹介するようにこの s.wam という関数の中では二番目の引数は y なので、この z は目的変数となります（そのものではありませんが）。
続いて以下の処理が行われますが、これが gam.fit の本体となります：
### ここで bf.call が評価される。 s.wam の場合、bakfit が呼ばれる。 ### bf.call で指定されている smooth.frame はこの時点では単なる data.frame fit <- eval(bf.call)
eval によって bf.call が評価されるので、 s.wam が実行されます。 s.wam を見に行くまえに、このループ内での残りの処理を簡単に見ておきましょう。
### 予測値にオフセットを加算する eta <- fit$fitted.values + offset ### eta から mu に変換する mu <- linkinv(eta) ### デビアンスを更新 old.dev <- new.dev new.dev <- sum(dev.resids(y, mu, weights)) ### ここの trace は対角和ではなく、 gam のオプションでループごとのデビアンスをモニターできる if (trace) cat("GAM ", bf, " loop ", iter, ": deviance = ", format(round(new.dev, digits)), " \n", sep = "") ### ループの打ち切り判定 ### デビアンスが NA となった場合、警告を出して打ち切る if (is.na(new.dev)) { one.more <- FALSE warning("iterations terminated prematurely because of singularities") } ### 差が十分に小さければ終了 else one.more <- abs(old.dev - new.dev)/(old.dev + 0.1) > epsilon if (!one.more) break }
評価結果を受けとり、デビアンスをもとにループを継続するかを判定しています。 gam.fit の残りの工程は一旦無視して、 s.wam を見に行きましょう。また次回。
GAMをもう少し理解したい
とても久しぶりの更新です。
背景
業務でモデリングを行うとき、私は大抵の場合GLMから始めます。目的変数に合わせて柔軟に分布を選択することが可能で、回帰係数という極めて解釈性の高い結果を得ることができるというのが理由です。
一方でGLMを使っていて不満に感じることの一つが、（ の世界で）非線形な効果を表現できないという点です。もちろん２次・３次の項や交互作用項を追加することである程度それらの不満は解消できるのですが、もう少しデータからそれらの特徴を学習したいと思うことがあります。
今回取り上げる一般化加法モデル（Generalized Additive Model, GAM）は、そのような複雑な関連性を表現できるよう説明変数に非線形な変換を行うもので、GLMを拡張したものとなります。ちょっと古いですが、2015年にMicrosoft RearchがKDDに出した論文(PDF)では、GAMを指して「the gold standard for intelligibility when low-dimensional terms are considered」と言っており、解釈性を保ちつつ高い予測精度を得ることができるモデルとしています。なおこの論文ではGAMに２次までの交互作用を追加するGA2Mという手法を提案しています。
このGAMの実装について調べた内容を書き留めておきます。GAMがどういうものかとか、平滑化に関する説明は、他に良いページがありますのでそちらを参照してください。例えば：
- https://logics-of-blue.com/%E5%B9%B3%E6%BB%91%E5%8C%96%E3%82%B9%E3%83%97%E3%83%A9%E3%82%A4%E3%83%B3%E3%81%A8%E5%8A%A0%E6%B3%95%E3%83%A2%E3%83%87%E3%83%AB/
- http://testblog234wfhb.blogspot.com/2014/07/generalized-additive-model.html
- https://blog.datarobot.com/jp/2017/10/24/ga2m-and-rating-table
GAMの実行結果
始めに、GAMを使うとどのような結果を得ることができるのか確認しましょう。ちなみに検証した環境は以下の通りです。
> sessionInfo() R version 3.6.0 (2019-04-26) Platform: x86_64-apple-darwin15.6.0 (64-bit) Running under: macOS Mojave 10.14.6 Matrix products: default BLAS: /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib LAPACK: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRlapack.dylib locale: [1] ja_JP.UTF-8/ja_JP.UTF-8/ja_JP.UTF-8/C/ja_JP.UTF-8/ja_JP.UTF-8 attached base packages: [1] stats graphics grDevices utils datasets methods base loaded via a namespace (and not attached): [1] compiler_3.6.0 tools_3.6.0 knitr_1.23 xfun_0.7
GAMの実装としてRでは {mgcv}が使われることが多いようですが、今回は「Rによる統計的学習入門」を参考に{gam}を使用しました。ちなみに、この本は統計・機械学習の主要な手法を網羅的に押さえつつ章末にRでの実行方法が紹介されており、大変勉強になる良書です。
- 作者: Gareth James,Daniela Witten,Trevor Hastie,Robert Tibshirani,落海浩,首藤信通
- 出版社/メーカー: 朝倉書店
- 発売日: 2018/08/03
- メディア: 単行本（ソフトカバー）
- この商品を含むブログ (1件) を見る
同書で用いているサンプルデータ Wage を使用するため、{ISLR}を同時に読み込みます。この{ISLR}パッケージは上記の本の原著であるIntroduction of Statistical Learning with Applications in Rから来ているようです。またこのデータは、アメリカの大西洋岸中央部における男性3000人の賃金、および年齢や婚姻状況、人種や学歴などの属性が記録されています。
library(gam) library(ISLR)
データの中身を見てみましょう。
> head(Wage) year age maritl race education region jobclass 231655 2006 18 1. Never Married 1. White 1. < HS Grad 2. Middle Atlantic 1. Industrial 86582 2004 24 1. Never Married 1. White 4. College Grad 2. Middle Atlantic 2. Information 161300 2003 45 2. Married 1. White 3. Some College 2. Middle Atlantic 1. Industrial 155159 2003 43 2. Married 3. Asian 4. College Grad 2. Middle Atlantic 2. Information 11443 2005 50 4. Divorced 1. White 2. HS Grad 2. Middle Atlantic 2. Information 376662 2008 54 2. Married 1. White 4. College Grad 2. Middle Atlantic 2. Information health health_ins logwage wage 231655 1. <=Good 2. No 4.318063 75.04315 86582 2. >=Very Good 2. No 4.255273 70.47602 161300 1. <=Good 1. Yes 4.875061 130.98218 155159 2. >=Very Good 1. Yes 5.041393 154.68529 11443 1. <=Good 1. Yes 4.318063 75.04315 376662 2. >=Very Good 1. Yes 4.845098 127.11574
このデータを用いて早速フィッティングしてみましょう。 glm と同様に関数 gam でモデルを当てはめることができます。
res_gam <- gam(wage ~ s(year, 4) + s(age, 5) + education, data = Wage)
ここで s(age, 5) は、説明変数 year を平滑化した上でモデルに取り込むことを意味し、5は平滑化の自由度です。なお s() の時点ではまだ平滑化は行われておらず、平滑化に必要な情報を属性として付与しているだけのようです。
> head(s(Wage$age, 5)) [1] 18 24 45 43 50 54 attr(,"spar") [1] 1 attr(,"df") [1] 5 attr(,"call") gam.s(data[["s(Wage$age, 5)"]], z, w, spar = 1, df = 5) attr(,"class") [1] "smooth"
ではフィッティングした結果を見てみましょう。
> summary(res_gam) Call: gam(formula = wage ~ s(year, 4) + s(age, 5) + education, data = Wage) Deviance Residuals: Min 1Q Median 3Q Max -119.43 -19.70 -3.33 14.17 213.48 (Dispersion Parameter for gaussian family taken to be 1235.69) Null Deviance: 5222086 on 2999 degrees of freedom Residual Deviance: 3689770 on 2986 degrees of freedom AIC: 29887.75 Number of Local Scoring Iterations: 2 Anova for Parametric Effects Df Sum Sq Mean Sq F value Pr(>F) s(year, 4) 1 27162 27162 21.981 2.877e-06 *** s(age, 5) 1 195338 195338 158.081 < 2.2e-16 *** education 4 1069726 267432 216.423 < 2.2e-16 *** Residuals 2986 3689770 1236 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Anova for Nonparametric Effects Npar Df Npar F Pr(F) (Intercept) s(year, 4) 3 1.086 0.3537 s(age, 5) 4 32.380 <2e-16 *** education --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
パッと見ると、 glm を当てはめたときと同様の結果が得られるようです。実際、 （2019/10/17修正 今回のケースでは lm.wfit は使われないので正しくありませんでした。）gam の本体（と思われる） gam.fit の中では stats::lm.wfit が呼ばれます。
しかし、 glm の結果と大きく異なる点として、各説明変数の回帰係数が出ていません。これはどういうことでしょう。以下のように説明変数の効果をプロットしてみます。
par(mfrow = c(1, 3)) plot(res_gam, se = TRUE, col = "blue")
真ん中の age が分かりやすいですが、この変数は34まではy軸の値が0未満となっており、 age は（ year や education を固定した下で）35歳程度まで wage の平均を下回っていることがわかります。その後45歳をピークに減少傾向に入り、65歳を過ぎると再び平均を下回るようになり、以降は急激に減少していきます。このような現象は、年齢とともに役職が上がることで賃金が増加し、退職によって減少することを考えれば非常に納得感のあるものだと思います。
なお上のプロットに必要なデータは以下のように取れるので、説明変数の各点における目的変数に対する影響を数値でも確認できます（必要な関数がエクスポートされていないので gam::: で直接呼び出しています）。
tmp <- gam:::preplot.Gam(res_gam, terms = gam:::labels.Gam(res_gam)) age_sm <- cbind(tmp$`s(age, 5)`$x, tmp$`s(age, 5)`$y) age_sm_uniq <- unique(age_sm[order(age_sm[, 1]), ])
プロットしてみましょう。同じ線が描けます。
plot(age_sm_uniq, type = "l", xlab = "age", ylab = "s(age, 5)")
さて、上記の結果は例えば age の二次の項を含めることで lm でも再現できるかもしれません。例えば以下のようになります：
### lmでフィッティング res_lm <- lm(wage ~ year + poly(age, 2) + education, Wage) ## poly(age, 2)で二次の多項式とする ### 予測用のデータ作成。ageだけを変化させ、yearとeducationは固定する。 x_lm <- seq(min(Wage$age), max(Wage$age), length.out = 100) nd <- data.frame(year = rep(2003, 100), age = x_lm, education = rep("2. HS Grad")) prd_lm <- predict(res_lm, nd, se.fit = T) ### 予測値から平均を引いてからプロット prd_m <- 90 plot(x_lm, prd_lm$fit - prd_m, type = "l", col = "blue", ylim = c(-40, 10)) lines(x_lm, prd_lm$fit - prd_m + 2*prd_lm$se.fit, lty = "dashed") lines(x_lm, prd_lm$fit - prd_m - 2*prd_lm$se.fit, lty = "dashed")
大体同じようなプロットを作成する事ができました。しかしピークを過ぎてからの緩やかな減少は表現できていませんし、一つ一つの変数について何次までの多項式を含めるかを検討していくのは少し手間がかかります。GAMを使えばデータに存在する細やかな変化を自動的に捉えることができます（もちろん代償もあります）。
GAMの実装
gam()
それでは gam がどのようにフィッティングを行っているのかを見ていきましょう。本体は最後の方に出てくる gam.fit なのですが、途中も少し細かく追ってみます。コンソールで gam を実行すると、以下のように関数の中身を見ることができます。
まず以下のブロックでは、 gam にオプションとして指定した内容に沿った処理を実行しています。
function (formula, family = gaussian, data, weights, subset, na.action, start = NULL, etastart, mustart, control = gam.control(...), model = TRUE, method = "glm.fit", x = FALSE, y = TRUE, ...) { ### 関数の引数を名前付きで確定。gam(wage ~ s(year, 4) + education, Wage)として与えた場合、 ### formula = と data = がそれぞれ保持される。 ### match.call returns a call in which all of the specified arguments are specified by their full names. call <- match.call() ### familyの判定 if (is.character(family)) family <- get(family, mode = "function", envir = parent.frame()) if (is.function(family)) family <- family() if (is.null(family$family)) { print(family) stop("`family' not recognized") } ### データが指定されていない場合 if (missing(data)) data <- environment(formula) ### 指定されている引数の取り出し mf <- match.call(expand.dots = FALSE) m <- match(c("formula", "data", "subset", "weights", "etastart", "mustart", "offset"), names(mf), 0L) mf <- mf[c(1L, m)] ### 指定されていない引数を指定し、stats::model.frame()の形式に仕立てる mf$na.action = quote(na.pass) mf$drop.unused.levels <- TRUE mf[[1L]] <- quote(stats::model.frame)
次に、ここで一つ gam ならではの処理として平滑化に使う関数を取り出しています。
### 平滑化の関数を取ってくる（s, lo, random) gam.slist <- gam.smoothers()$slist
gam.smoothers として新しい平滑化関数を指定することも出来るようですが、デフォルトは s 、 lo および random で、それぞれ 平滑化スプライン、局所回帰、ランダム効果としての指定を意味しているようです。 3つめの random がわからなかったので調べてみたところ、これはカテゴリ変数に対する指定で、パラメータの推定においていわゆる縮小推定を行なうもののようでした。
https://www.rdocumentation.org/packages/gam/versions/1.16.1/topics/random
これらに接頭語として gam を加えた（e.g. gam.s）関数が平滑化のための関数として実行されます。 gam.s については後述します。
次のブロックですが、 call クラスであった mf を評価することで data.frame に変換しています。
### term を mf$formula に渡す mt <- if (missing(data)) terms(formula, gam.slist) else terms(formula, gam.slist, data = data) mf$formula <- mt ### ここで mf 、つまり model.frame が実行されて data.frame になる ### ただし平滑化は実行されず、平滑化のパラメータは attribute として持っている mf <- eval(mf, parent.frame()) if (missing(na.action)) { naa = getOption("na.action", "na.fail") na.action = get(naa) } mf = na.action(mf) mt = attributes(mf)[["terms"]]
ここまでは mf は call クラス、すなわち未評価の関数およびその引数を要素に持つオブジェクトでした。それが eval で評価されたため stats::model.frame が実行され、 formula にしたがい data.frame が生成された、という流れのようです（間違っていたらすみません）。
### method の指定によって処理を分ける。 glm.fit または glm.fit.null 以外の場合はエラー switch(method, model.frame = return(mf), glm.fit = 1, glm.fit.null = 1, stop("invalid `method': ", method))
ここでは method によって処理を分けていますが、現状は glm.fit または model.frame のみ受け付けているようです。 なおhelpを参照すると、 model.frame を指定した場合フィッティングは行われないようですね。以下のようになります。
> gam(wage ~ s(year, 4) + s(age, 5) + education, data = Wage, method = "model.frame") wage s(year, 4) s(age, 5) education 231655 75.04315 2006 18 1. < HS Grad 86582 70.47602 2004 24 4. College Grad 161300 130.98218 2003 45 3. Some College 155159 154.68529 2003 43 4. College Grad 11443 75.04315 2005 50 2. HS Grad
この通りデータが返ってきます。
以下では、YおよびXをそれぞれ抽出し、さらにフィッティングに必要なオプションを指定しています。
### Y を取り出す Y <- model.response(mf, "any") ### X を matrix で取り出す。 ### gam を実行したときのエラーメッセージ（ `non-list contrasts argument ignored` ）はここで出ている。 contrasts の指定が良くない様子。 X <- if (!is.empty.model(mt)) model.matrix(mt, mf, contrasts) else matrix(, NROW(Y), 0) ### その他パラメータ（weights, offset, mustart, etastart） weights <- model.weights(mf) offset <- model.offset(mf) if (!is.null(weights) && any(weights < 0)) stop("Negative wts not allowed") if (!is.null(offset) && length(offset) != NROW(Y)) stop("Number of offsets is ", length(offset), ", should equal ", NROW(Y), " (number of observations)") mustart <- model.extract(mf, "mustart") etastart <- model.extract(mf, "etastart")
以上で準備が完了し、 gam.fit で当てはめを行います。
### ここが本体。 gam.fit を呼び出している fit <- gam.fit(x = X, y = Y, smooth.frame = mf, weights = weights, start = start, etastart = etastart, mustart = mustart, offset = offset, family = family, control = control)
いったん後続の部分は無視して gam.fit に移りたいと思いますが、長くなったので今回はここまでにします。
「ガウス過程と機械学習」第3章のグラフ（一部）を作図する④
前回の記事の続きです。
今回は図3.23に挑戦します。データはこちらから該当する部分を取ってきました。
まずはプロットしてみます。一部のデータは除外しました。
dat <- read.table("./Data/World_Record.dat", sep = "\t", col.names = c("name", "time", "date")) dat <- dat[-c(6, 17:21), -1] # ドーピング！ dat$date <- as.integer(gsub(".*, ", "", dat$date)) colnames(dat) <- c("y", "x") plot(dat$x, dat$y, xlab = "year", ylab = "time", xlim = c(1960, 2020), ylim = c(9.5, 10.1), pch = 4, cex = 2)
本で使われているデータとは少し違っている様子ですが、このまま進めます。
本にしがたい、Xおよびyを平均0、分散1にスケーリングします。scaleを使いますが、datはdata.frameで持っておきたいのでスケーリングのパラメータは別に取っておきましょう。
dat <- scale(dat) scale_pars <- unlist(attributes(dat)[c("scaled:center", "scaled:scale")]) dat <- as.data.frame(dat)
それではこのデータを使ってガウス過程回帰を実行します。まずはRBFカーネルを使って推定してみましょう。
前回の記事で定義した関数を色々使い回すことにします。まずはL、これはハイパーパラメータを与えたときに対数尤度を返す関数でした。
L <- function(param, x, y) { # C <- -nrow(train)/2*log(2*pi) C <- 0 K <- get_cov_mat_exp(x, x, theta1 = param[1], theta2 = param[2]) diag(K) <- diag(K) + exp(param[3]) return(-log(det(K)) - t(as.matrix(y)) %*% solve(K) %*% as.matrix(y) + C) }
Lの中で使われているget_cov_mat_expも定義しましょう。これは共分散行列を得る関数です。
get_cov_mat_exp <- function(x1, x2, theta1 = 1, theta2 = 1) { ## matrixに変換 x1 <- as.matrix(x1) x2 <- as.matrix(x2) ## 行の組み合わせを作成する n <- nrow(x1) m <- nrow(x2) d <- ncol(x1) tmp <- cbind(kronecker(x1, matrix(1, m)), kronecker(matrix(1, n), x2)) ret <- apply(tmp, 1, function(x, d, theta1, theta2) { rbf_knl_exp(x[(1:d)], x[(d+1):ncol(tmp)], theta1, theta2) # rbf_knl_expを使う }, d, theta1, theta2) return(matrix(ret, n, m, byrow = T)) }
さらにrbf_knl_expも使います。
rbf_knl_exp <- function(x1, x2, theta1 = 1, theta2 = 1) { exp(theta1) * exp(-norm(x1 - x2, "2")^2/exp(theta2)) }
上記の関数は、パラメータ探索における効率化のためにパラメータのlogを取ってから与えることを前提としていましたが、可視化用に元の関数も定義しておきましょう。
get_cov_mat <- function(x1, x2, theta1 = 1, theta2 = 1) { ## matrixに変換 x1 <- as.matrix(x1) x2 <- as.matrix(x2) ## 行の組み合わせを作成する n <- nrow(x1) m <- nrow(x2) d <- ncol(x1) tmp <- cbind(kronecker(x1, matrix(1, m)), kronecker(matrix(1, n), x2)) ret <- apply(tmp, 1, function(x, d, theta1, theta2) { rbf_knl(x[(1:d)], x[(d+1):ncol(tmp)], theta1, theta2) }, d, theta1, theta2) return(matrix(ret, n, m, byrow = T)) }
rbf_knl <- function(x1, x2, theta1 = 1, theta2 = 1) { theta1 * exp(-norm(x1 - x2, "2")^2/theta2) }
ではこれらの関数を使い、今回のデータでハイパーパラメータの最適化を実行してみます。
param <- c(1, 1, 1) res_par <- optim(par = optim(par = param, fn = L, x = dat$x, y = dat$y, control = list(fnscale = -1))$par, fn = L, x = dat$x, y = dat$y, control = list(fnscale = -1))
> print(exp(res_par$par), digits = 3) [1] 2.877 6.887 0.101 > L(res_par$par, dat$x, dat$y) [,1] [1,] 6.582661
教科書P94に示された値とは結構異なるようです（それぞれ1.62、0.44、0.06）。特にtheta2が大きく推定されていますね。
このパラメータを使って可視化します。やっつけですが、以下のように関数を定義しました。
gen_gpreg_plot <- function(param, x, y) { min_y <- 9.4 max_y <- 10.1 min_x <- -2 max_x <- 1.7 test <- seq(min_x, max_x, 0.05) theta1 <- exp(param[1]) theta2 <- exp(param[2]) theta3 <- exp(param[3]) K <- get_cov_mat(x, x, theta1 = theta1, theta2 = theta2) diag(K) <- diag(K) + theta3 k <- get_cov_mat(x, test, theta1 = theta1, theta2 = theta2) s <- get_cov_mat(test, test, theta1 = theta1, theta2 = theta2) diag(s) <- diag(s) + theta3 Mu <- t(k) %*% solve(K) %*% y * scale_pars["scaled:scale.y"] + scale_pars["scaled:center.y"] Var <- (s - t(k) %*% solve(K) %*% k) * (scale_pars["scaled:scale.y"])^2 CI_u <- Mu + 2 * sqrt(diag(Var)) CI_l <- Mu - 2 * sqrt(diag(Var)) x <- x * scale_pars["scaled:scale.x"] + scale_pars["scaled:center.x"] y <- y * scale_pars["scaled:scale.y"] + scale_pars["scaled:center.y"] test <- test * scale_pars["scaled:scale.x"] + scale_pars["scaled:center.x"] min_x <- min_x * scale_pars["scaled:scale.x"] + scale_pars["scaled:center.x"] max_x <- max_x * scale_pars["scaled:scale.x"] + scale_pars["scaled:center.x"] xlim <- c(min_x, max_x) ylim <- c(min_y, max_y) plot(x, y, xlim = xlim, ylim = ylim, type = "n") polygon(c(test, rev(test)), c(CI_u, rev(CI_l)), col = 'gray', border = NA) points(x, y, xlim = xlim, ylim = ylim, xlab = "x", pch = 4, cex = 2) lines(test, Mu, xlim = xlim, ylim = ylim, type = "l", ylab = "") }
上の関数の中で、ガウス過程によりK、kおよびsを求めるところまでは標準化後の値を使っていますが、Muを求めるところからは元のスケールに戻しています。
プロットしてみましょう。
gen_gpreg_plot(res_par$par, dat$x, dat$y)
おや、随分とスムーズな線になってしまいました。試しに教科書のパラメータを使ってみましょう。
book_par <- c(log(1.62), log(0.44), log(0.06)) gen_gpreg_plot(book_par, dat$x, dat$y)
教科書と同じようなプロットが得られました。
ガウス過程回帰の予測値（Mu）が滑らかであるということは、隣接する値が似ていることを意味しています。つまりそれらの共分散が（相対的に）大きい状態です。
RBFカーネルではtheta1とtheta2がパラメータとして与えられますが、theta2はexpの中で使われるため、特に共分散に対する影響が大きいのではないでしょうか。今回の例で言えば（expの中で分母として使われる）theta2が教科書の値よりも10倍以上大きな値となっているので、expの項が0に近づきやすかったのではないかと思います。
ちょっと確認してみましょう。theta2を教科書の値と差し替えた場合の共分散行列の変化をプロットします。
library(gplots) a <- get_cov_mat(dat$x, dat$x, res_par$par[1], res_par$par[2]) b <- get_cov_mat(dat$x, dat$x, res_par$par[1], 0.44)
heatmap.2(a, Rowv = NA, Colv = NA, dendrogram = "none", trace = "none") heatmap.2(b, Rowv = NA, Colv = NA, dendrogram = "none", trace = "none")
これを見ると、データから推定した値による共分散行列（a）はdat$xの1~5番目と6~11番目の値の共分散を大きく評価していますが、教科書の値を使った場合（b）では「似ていない」と判断しています。その結果、データから推定した方では相対的にグネグネした線が描かれたのだと思います。
続いてRBFカーネルに線形カーネルを追加して当てはめてみたいと思います。対数尤度関数を以下のように修正します：
L2 <- function(param, x, y) { # C <- -nrow(train)/2*log(2*pi) C <- 0 K <- get_cov_mat_exp(x, x, theta1 = param[1], theta2 = param[2]) diag(K) <- diag(K) + exp(param[3]) K <- K + exp(param[4]) * outer(x, x) # 線形カーネルを追加 return(-log(det(K)) - t(as.matrix(y)) %*% solve(K) %*% as.matrix(y) + C) }
optimによる最適化を実行しましょう。
param <- c(1, 1, 1, 1) res_par_2 <- optim(par = optim(par = param, fn = L2, x = dat$x, y = dat$y, control = list(fnscale = -1))$par, fn = L2, x = dat$x, y = dat$y, control = list(fnscale = -1))
> print(exp(res_par_2$par), digits = 3) [1] 0.147 1.234 0.104 0.772 > L2(res_par_2$par, dat$x, dat$y) [,1] [1,] 9.436611
相変わらずtheta2が教科書の値と外れますが、先ほどよりは小さな値となりました。しかし対数尤度は一致しないですね。。。
可視化用の関数を少し修正します。
gen_gpreg_plot_2 <- function(param, x, y) { min_y <- 9.4 max_y <- 10.1 min_x <- -2 max_x <- 1.7 test <- seq(min_x, max_x, 0.05) theta1 <- exp(param[1]) theta2 <- exp(param[2]) theta3 <- exp(param[3]) theta4 <- exp(param[4]) K <- get_cov_mat(x, x, theta1 = theta1, theta2 = theta2) diag(K) <- diag(K) + theta3 K <- K + theta4 * outer(x, x) k <- get_cov_mat(x, test, theta1 = theta1, theta2 = theta2) k <- k + theta4 * outer(x, test) s <- get_cov_mat(test, test, theta1 = theta1, theta2 = theta2) diag(s) <- diag(s) + theta3 s <- s + theta4 * outer(test, test) Mu <- t(k) %*% solve(K) %*% y * scale_pars["scaled:scale.y"] + scale_pars["scaled:center.y"] Var <- (s - t(k) %*% solve(K) %*% k) * (scale_pars["scaled:scale.y"])^2 CI_u <- Mu + 2 * sqrt(diag(Var)) CI_l <- Mu - 2 * sqrt(diag(Var)) x <- x * scale_pars["scaled:scale.x"] + scale_pars["scaled:center.x"] y <- y * scale_pars["scaled:scale.y"] + scale_pars["scaled:center.y"] test <- test * scale_pars["scaled:scale.x"] + scale_pars["scaled:center.x"] min_x <- min_x * scale_pars["scaled:scale.x"] + scale_pars["scaled:center.x"] max_x <- max_x * scale_pars["scaled:scale.x"] + scale_pars["scaled:center.x"] xlim <- c(min_x, max_x) ylim <- c(min_y, max_y) plot(x, y, xlim = xlim, ylim = ylim, type = "n") polygon(c(test, rev(test)), c(CI_u, rev(CI_l)), col = 'gray', border = NA) points(x, y, xlim = xlim, ylim = ylim, xlab = "x", pch = 4, cex = 2) lines(test, Mu, xlim = xlim, ylim = ylim, type = "l", ylab = "") }
プロットしてみます。
gen_gpreg_plot_2(res_par_2$par, dat$x, dat$y)
依然としてかなりスムーズな線が引かれています。教科書の値を使ってみましょう。
book_par <- c(log(0.07), log(0.02), log(0.06), log(0.92)) gen_gpreg_plot_2(book_par, dat$x, dat$y)
同じようなプロットが得られました！
以上です。 今回まで「ガウス過程と機械学習」第3章のグラフをいくつか作図してきましたが、これ以降はMCMCやGPyを使うなどしており、ちょっと理解に時間がかかりそうなので一旦ここで終わりにしておきます。
