統計コンサルの議事メモ

統計や機械学習の話題を中心に、思うがままに

glmnetをもう少し理解したい④

それでは前回の記事に続いてelnet1の紹介です。過去の記事はこちらです。

ushi-goroshi.hatenablog.com ushi-goroshi.hatenablog.com ushi-goroshi.hatenablog.com

elnet1の実装

前回の記事で最後に触れた通り、elnet1自体は 180 行程度とそれほど大きくはないサブルーチンなのですが、多数のループが込み入っています。 具体的には以下の通り 9 つのループ処理(fortran なので do 文)がネストした構造となっており、 しかもgotoによって行き来しています(わかりやすいように R で書いてありますが、添字は統一してあります)。

# 1番目
for (m in 1:nlam) {
  # 2番目のループ
  for (j in 1:ni) {
  }
  # 3番目のループ
  for (k in 1:ni) {
    # 4番目のループ
    for (j in 1:ni) {
    }
    # 5番目のループ
    for (j in 1:ni) {
    }
  }
  # 6番目のループ
  for (l in 1:nin) {
    # 7番目のループ
    for (j in 1:nin) {
    }
  }
  # 8番目のループ
  for (j in 1:ni) {
  }
  # 9番目のループ
  for (j in 1:nin) {
  }
}

前処理

まずはいつもの通り変数の定義ですが、それに加えて初期パラメータを取得するという処理が入ります。

subroutine elnet1(beta,ni,ju,vp,cl,g,no,ne,nx,x,nlam,flmin,ulam,th
                  *r,maxit,xv,  lmu,ao,ia,kin,rsqo,almo,nlp,jerr)
implicit double precision(a-h,o-z)
double precision vp(ni),g(ni),x(no,ni),ulam(nlam),ao(nx,nlam)
double precision rsqo(nlam),almo(nlam),xv(ni)
double precision cl(2,ni)
integer ju(ni),ia(nx),kin(nlam)
double precision, dimension (:), allocatable :: a,da
integer, dimension (:), allocatable :: mm
double precision, dimension (:,:), allocatable :: c
allocate(c(1:ni,1:nx),stat=jerr)
if(jerr.ne.0) return;

! 初期パラメータを取得
call get_int_parms(sml,eps,big,mnlam,rsqmax,pmin,exmx,itrace)

! a, mm, da を allocate
allocate(a(1:ni),stat=jerr)  ! a は説明変数の数の次元をもつベクトル

if(jerr.ne.0) return
allocate(mm(1:ni),stat=jerr) ! mm は説明変数の数の次元をもつベクトル

if(jerr.ne.0) return
allocate(da(1:ni),stat=jerr)
if(jerr.ne.0) return

ここでget_int_parmsはそれほど大きくないので全体を見てみましょう。 以下のようなサブルーチンです:

subroutine get_int_parms(sml,eps,big,mnlam,rsqmax,pmin,exmx,itrace)
implicit double precision(a-h,o-z)                                
data sml0,eps0,big0,mnlam0,rsqmax0,pmin0,exmx0,itrace0  /1.0d-5,1.0d-6,9.9d35,5,0.999,1.0d-9,250.0,0/  
  sml=sml0                                                          
eps=eps0                                                          
big=big0                                                          
mnlam=mnlam0                                                      
rsqmax=rsqmax0                                                    
pmin=pmin0                                                        
exmx=exmx0                                                        
itrace=itrace0                                                    
return                                                            
entry chg_fract_dev(arg)                                          
sml0=arg                                                          
return                                                            
entry chg_dev_max(arg)                                            
rsqmax0=arg                                                       
return                                                            
entry chg_min_flmin(arg)                                          
eps0=arg                                                          
return                                                            
entry chg_big(arg)                                                
big0=arg                                                          
return                                                            
entry chg_min_lambdas(irg)                                        
mnlam0=irg                                                        
return                                                            
entry chg_min_null_prob(arg)                                      
pmin0=arg                                                         
return                                                            
entry chg_max_exp(arg)                                            
exmx0=arg                                                         
return                                                            
entry chg_itrace(irg)                                             
itrace0=irg                                                       
return                                                            
end

上から3行目のdata文は変数に初期値を与える fortran の記法のようで、dataに続いて宣言した変数に対して/で挟んだ値を初期値として与えるようです。 そのためsml0には 1.0d-5 が、eps0には 1.0d-6 が入力されます。 ここで d は倍精度の指数表記を表します。13行目のentry以降は各変数について特定の値を指定するためのもののようです(entryの使い方がよくわからない…)。

続けていくつか変数に値を代入します。 まずはbtaですが、代入しているbetaは元々parmとして渡されたもので、これはelnet.rparm = alphaとして渡していたものでした。さらにこのalphaglmnet.rで定義されたもので、L1 と L2 それぞれに対する罰則の配分を決めるパラメータです:

(なぜか Tex が表示されないのでひとまず)

(1 − \alpha)/2 ||\beta||^2_2 + \alpha||\beta||_1

      bta=beta

このbtaを 1 から減じたものをombとしますが、このombはすぐ下で定義されるalmとの乗算でdemを定義する(つまりdem = alm * obm)ためだけに使われています。 さらにalmはループの中で更新されながら最終的にはbtaとの乗算によってabとなり、回帰係数の縮小に使われることになります。 またその次のalfalmの更新に使われますので、これらの変数がループの中で更新されつつ回帰係数の縮小に利用されるということになります(他にもあります)。

      omb=1.0-bta
      alm=0.0
      alf=1.0

以下のブロックではeqsalfを定義しますが、flminが 1.0以上であればスキップされるようです。 このflminというのはglmnet.rにおいて罰則lambdaが指定されていれば 1 が、されていない時にはlambda.min.ratioが入力される変数でした。 lambda.min.ratioはデフォルトではlambda.min.ratio = ifelse(nobs < nvars, 0.01, 1e-04)となっていますので 1 よりは小さい値が入りそうです。 したがって以下のブロックは「lambdaが指定されていないときはalfを定義しよう」という処理になっています(eqsはここしか出てきません)。

その場合、epsflmin(=1)の大きい方を新たにeqsと定義しますが、このepsget_int_parmseps0(1.0d-6 という小さい数)を受け取っていました。 一方lambda.min.ratioは先ほど述べたようにデフォルトではlambda.min.ratio = ifelse(nobs < nvars, 0.01, 1e-04)となっていますので、もう少し大きい値となりそうです。 したがってeqsは 0.01 or 1e-04 、alfはその1/(nlam-1)乗となるようです。

      if(flmin .ge. 1.0)goto 10271
      eqs=max(eps,flmin)
      alf=eqs**(1.0/(nlam-1))  ! alf を eqs の (1/(nlam-1)) で定義する

flminが 1 以上である(lambdaが指定されている)場合は上記をスキップしてこちらにきます。rsqはそのまま残差平方和ですね。

続くaelnet1の中で重要な役割を担っているのでじっくりと見ていきましょう。 実はこのaは(縮小された)回帰係数を格納する変数です。

10271 continue
      ! パラメータの初期化
      rsq=0.0 ! 残差平方和
      a=0.0

このaがどうなるのか、フライングして先の処理を見てみましょう。 elnet1の 70 行目前後に以下の処理があります:

      ak=a(k)                                                           
      u=g(k)+ak*xv(k)                                                   
      v=abs(u)-vp(k)*ab                                                 
      a(k)=0.0                                                          
      if(v.gt.0.0) a(k)=max(cl(1,k),min(cl(2,k),sign(v,u)/(xv(k)+vp(k)*dem)))
      if(a(k).eq.ak)goto 10371                                  

akという変数にaの k 番目の値を渡しておき、uvを定義し、aの k 番目の値を 0 に更新した上で色んな値を参照しながら再度更新しています(このuvは後で確認します)。 最終的にaは以下のようにaoという変数に代入されます(154 行目):

      if(nin.gt.0) ao(1:nin,m)=a(ia(1:nin)) 

このaoですが、elnetuの中でelnet1を呼び出すときにはcaという引数として渡されています。

! elnet1 で受け取る変数
! lmu の次に ao がある
subroutine elnet1(beta,ni,ju,vp,cl,g,no,ne,nx,x,nlam,flmin,ulam,th
     *r,maxit,xv,  lmu,ao,ia,kin,rsqo,almo,nlp,jerr)

! elnetu で elnet1 を call するときの引数
! こちらは lmu の次に ca がある
call elnet1(parm,ni,ju,vp,cl,g,no,ne,nx,x,nlam,flmin,vlam,thr,maxi,xv,  lmu,ca,ia,nin,rsq,alm,nlp,jerr)

このcaelnet.rの中で.Fortran("elnet", ...)と call される際に定義される変数でした:

else .Fortran("elnet", ka, parm = alpha, nobs, nvars, as.double(x), 
              y, weights, jd, vp, cl, ne, nx, nlam, flmin, 
              ulam, thresh, isd, intr, maxit, lmu = integer(1), 
              a0 = double(nlam), 
              # ここで ca が定義されている
              ca = double(nx * nlam), 
              ia = integer(nx), 
              nin = integer(nlam),  rsq = double(nlam), alm = double(nlam), 
              nlp = integer(1), jerr = integer(1), PACKAGE = "glmnet")

ここでnxは説明変数の数、nlamは罰則lambdaの数なので、説明変数の数 × lambda の数のベクトルを定義しています(そしてそれがelnet1の中でaoとして評価・格納される)。

このcaelnet.rの後続の処理において以下の箇所で抽出されます:

outlist = getcoef(fit, nvars, nx, vnames)

ここでglmnet:::getcoefは以下の通りで、fitとして返ってきたオブジェクトのcaそのものをbetaに格納しています(ninmaxが 0 の場合は 0 のベクトルが返る)。

# glmnet:::getcoef
function (fit, nvars, nx, vnames) {
  # ここまで省略
  nin = fit$nin[seq(lmu)]
  ninmax = max(nin)
  # ここまで省略
  if (ninmax > 0) {
    # ここで ca を抽出している 
    ca = matrix(fit$ca[seq(nx * lmu)], nx, lmu)[seq(ninmax), 
                                                , drop = FALSE]
    df = apply(abs(ca) > 0, 2, sum)
    ja = fit$ia[seq(ninmax)]
    oja = order(ja)
    ja = rep(ja[oja], lmu)
    ia = cumsum(c(1, rep(ninmax, lmu)))
    # beta に格納する
    beta = drop0(new("dgCMatrix", Dim = dd, Dimnames = list(vnames, 
                                                            stepnames), x = as.vector(ca[oja, ]), p = as.integer(ia - 
                                                                                                                   1), i = as.integer(ja - 1)))
  }
  else {
    beta = zeromat(nvars, lmu, vnames, stepnames)
    df = rep(0, lmu)
  }
  # ここも省略
  list(a0 = a0, beta = beta, df = df, dim = dd, lambda = lam)
}

これにいくつかの情報を追加したものがglmnetの返り値です。elnet1において評価されたaaoに格納され、elnetcaとして渡され、elnet.rbetaに抽出・格納される流れが伝わりましたでしょうか。

重要な変数を説明したところなので、以下ブロックで初期化している変数の詳細は出てきたときに説明するとして、さっさと次に進んでしまいましょう。

      mm=0
      nlp=0
      nin=nlp
      
      iz=0
      mnl=min(mnlam,nlam)

ループ①(almの更新)

上記までで必要な変数の初期化が完了したので、以下よりループに入ります。 一番外側のループはlambdaの個数(nlam)に対して実行されますが、nlamのデフォルトは 100 となっています(glmnet.r)。

以下ではおおよそalmを更新する処理を行うのですが、lambdaの指定の有無や、ループの回数によってalmに入力する値を変えています。

まずはlambdaの指定の有無で処理を分けます。以下のまとまりはflminが 1.0 より小さい場合にスキップされますが、先ほど述べたように、flminglmnet.rにおいてlambdaの指定がない場合に相当します。 lambdaの指定がある場合にはalm = ulam(m)としてalmを更新した上で、10291 までスキップするのですが、この 10291 は 2 番目のループの中にありますので、少し大きめのスキップとなるようです。 なおulamlambdaが指定されている場合、lambdaの降順になっているため、ループの 1 回目であればlambdaの最大値が入ります。

      do 10281 m=1,nlam ! nlambda なので lambda の個数だけループ

      if(itrace.ne.0) call setpb(m-1)  ! プログレスバー
      if(flmin .lt. 1.0)goto 10301
      alm=ulam(m) ! flmin が 1.0 以上の場合は alm = ulam(m) とする
      goto 10291

lambdaの指定がなければ以下の処理に入るのですが、ここではループの回数によってalmに入力する値を変えています。 具体的には、ループの 1 回目にはbig(9.9d35)という極端に大きな値を入力し、 2 回目には 0.0 を、3 回目以降は 元の値にalfを乗じたものを入力します。

10301 if(m .le. 2)goto 10311 ! ループの1回目と2回目はここをスキップ
      alm=alm*alf ! ループの3回目からは alm を alf を乗じる
      goto 10291
10311 if(m .ne. 1)goto 10321 ! ループの2回目はここをスキップ
      alm=big     ! ループの1回目は alm = big(9.9d35) にする 
      goto 10331
10321 continue
      alm=0.0     ! ループの2回目は alm を いったん 0 にする

このalfは先ほど説明した通りeqs^(1.0/(nlam-1))として定義されますが、eqsが 0.01 or 1e-4 とすると、nlambdaを 10 とした場合には以下のような数値になります:

0.01^(1/(10-1))
# [1] 0.5994843
1e-4^(1.0/(10-1))
# [1] 0.3593814

つまりalmはだんだん絶対値が小さくなるわけですね。

ループ②(罰則の定義)

続いて 2 番目のループに入ります…と言いつつ 2 番目のループは一瞬で終わります。 先ほど更新したalmについて変数ごとの内積と比較し、大きい方を採用します。 したがってここでは各変数に対するループとなります。

まずjuvpですが、juは前回記事で確認した通り、chkvarsによって各変数列の内容が全く同じでないかを確認したものでした。 ある変数列の中身が全く同じであれば 0 であったため、ここで次の変数にスキップされます。 次にvpですが、これは 1 回目の記事で確認した通りglmnet.rにおいて各変数に対する罰則の重み(デフォルトは 1) が入ったベクトルとして定義されたものです(vp = as.double(penalty.factor))。 罰則をかけない場合は 0 となり、スキップされるようです。

変数にバラつきがあり、罰則を検討する場合にはここで再度almを更新します。 ここで出てくるgstandardの中でyx内積(共分散)を格納したものとして定義されたものでした。 それを罰則の大きさで除しているため、penalty.factorを小さくする(分母が小さくなる)と共分散が大きくなり変数として残りやすい、というロジックになっているようですね。

ちなみにループ①の 1 回目のループはalmに 9.9d35 という数値が入るので必ずこの数値が採用されると思います。またループ 2 回目は今度はalmが 0.0 になるため、今度は必ず変数の共分散側の数値がalmになると思われます。

      ! 2番目のループ
      ! alm の更新
      do 10341 j=1,ni  ! ni は変数の数
      if(ju(j).eq.0) goto 10341
      if(vp(j).le.0.0) goto 10341
      alm=max(alm,abs(g(j))/vp(j))
10341 continue  ! 2番目のループここまで

上記の処理で変数を横断してalmを更新したのち、以下でさらにalmを更新します。 ここではbtaalpha; L1 と L2 への重みの配分パラメータ)と 0.001 の max で alm を除し、alfを乗じています。 一応ここで式を確認しておくと以下のようになります:

(なぜか Tex が表示されないのでひとまず)

alm = alf * alm/max(bta, 1.0d-3) = eqs^{(1.0/(nlam-1))} * alm/max((1-alpha), 0.001)

一体これは何をやっているんでしょうか…。

      continue
      alm=alf*alm/max(bta,1.0d-3)  

続いていくつかの変数を更新します。 demalm * ombとして定義されますが、ここでombは (1-bta)でした。 またabalmbtaを乗じたものですので、これらはそれぞれ「lambda×(1-alpha)」および「lambda×alpha」ということになり、demabが実質的な罰則の大きさを表すことになりそうですね。

10331 continue
10291 continue
      dem=alm*omb ! dem = alm * (1-bta)
      ab=alm*bta  ! ab = alm * bta

これらがどのように使われているか少し先を見てみましょう。

! ab
u=g(k)+ak*xv(k)   ! L69(ループ③の中)、L119(ループ⑥の中)
v=abs(u)-vp(k)*ab ! L70、L120(ともに上に同じ)

! dem
a(k)=0.0 ! L71、L121
if(v.gt.0.0) a(k)=max(cl(1,k),min(cl(2,k),sign(v,u)/(xv(k)+vp(k)*dem))) ! L72、L122

両方ともvpに乗じており、ababs(u)からの減算、demxv(k)との加算の後にsign(v,u)と除算し、clとの max/min を取っています。 vpは罰則の重みを定義したものでしたので、alphalambdaで決まる罰則の大きさをそのまま使うか弱くするかを決めています。 demの方は演算の結果をaに格納していますが、前述の通りaは回帰係数を保存する変数でしたので、sign(v,u)/(xv(k)+vp(k)*dem)cl(1,k)よりも大きければa、すなわち回帰係数が更新されるということになりますね。

またこの演算が実行されるかの基準としてvが使われており、このvを計算するためにabが使われている、ということのようです。 じゃあこのuとかvって何なの?という話なのですが、これは次のループの話なので少しお待ちください。

残る変数のうちrsq0は残差平方和ですね。またjzizと組み合わせて使われていますが、この条件分岐がちょっと理解出来なかったのでスキップします。 一応、izはループ①の途中(ループ③が終了した時点)で 1 になるため、iz * jzが 0 になるのはほぼjzが 0 の時に限ると言えそうです。 nlpは iteration のカウンターとして使われており、dlxは回帰係数の更新前後の差分を見ています。 どちらもループを抜けるための基準として使われています。

      rsq0=rsq 
      jz=1
      continue
10351 continue
      if(iz*jz.ne.0) goto 10360   ! iz = 0, jz = 1
      nlp=nlp+1 
      dlx=0.0

ちょっと長くなってしまったので一度きります。 次回はループ③から始めます。

glmnetをもう少し理解したい③

前回の記事では R の関数 elnet の中で elnet という Fortran のサブルーチンが呼ばれ(やっぱりややこしいですね)、さらに type.gaussian の値( covariancenaive )によって elnetuelnetn のいずれかが呼ばれるところまで確認しました。 今回は elnetu の中身を見ていきます。 過去の記事はこちらです。

ushi-goroshi.hatenablog.com

ushi-goroshi.hatenablog.com

elnetu の実装

それでは早速 elnetu を見ていきましょう。 elnetuelnet と同様にそれほど大きくないのでいきなり内容の確認に入りますが、処理としては以下の手順になっているようです:

  1. 前処理
  2. 標準化
  3. フィッティング
  4. 後処理

まずは前処理ですが、メモリの割り付けのあとに chkvars というサブルーチンを呼び出しています。

      subroutine elnetu(parm,no,ni,x,y,w,jd,vp,cl,ne,nx,nlam,  flmin,ulam,thr,isd,intr,maxit,  lmu,a0,ca,ia,nin,rsq,alm,nlp,jerr)
      implicit double precision(a-h,o-z)                                
      double precision x(no,ni),y(no),w(no),vp(ni),ulam(nlam),cl(2,ni)  
      double precision ca(nx,nlam),a0(nlam),rsq(nlam),alm(nlam)         
      integer jd(*),ia(nx),nin(nlam)                                    
      double precision, dimension (:), allocatable :: xm,xs,g,xv,vlam   
      integer, dimension (:), allocatable :: ju                         
      allocate(g(1:ni),stat=jerr)              
      if(jerr.ne.0) return                                              
      allocate(xm(1:ni),stat=jerr)                                      
      if(jerr.ne.0) return                                              
      allocate(xs(1:ni),stat=jerr)                                      
      if(jerr.ne.0) return                                              
      allocate(ju(1:ni),stat=jerr)                                     
      if(jerr.ne.0) return                                              
      allocate(xv(1:ni),stat=jerr)                                    
      if(jerr.ne.0) return                                              
      allocate(vlam(1:nlam),stat=jerr)                                  
      if(jerr.ne.0) return

      ! 1. 前処理
      call chkvars(no,ni,x,ju)

      if(jd(1).gt.0) ju(jd(2:(jd(1)+1)))=0
      if(maxval(ju) .gt. 0)goto 10071                                   
      jerr=7777                                                         
      return                                                            
10071 continue

この chkvars では x の各変数について一行目の値と異なる値が二行目以降にあるかを確認し、 ju に格納しています。

      subroutine chkvars(no,ni,x,ju)
      implicit double precision(a-h,o-z)
      double precision x(no,ni)
      integer ju(ni)
      
      ! ここから各変数のチェックを開始
      do 11061 j=1,ni
      ju(j)=0
      t=x(1,j) ! 1行目の値を取得

      ! ここから2行目の値を確認する
      do 11071 i=2,no
      ! t は x(1, j) なので、各変数 j について 1 行目の値と等しいかを確認している
      if(x(i,j).eq.t) goto 11071 ! 等しければ次の行へ
      ju(j)=1 ! 等しくない数値があれば ju を 1 にして次の変数へ
      goto 11072
11071 continue
11072 continue
11061 continue
      continue
      return
      end

異なる値がなければ全ての値は同じということになりますので、例えば回帰係数を推定する意味はありません。 後続の処理ではこの ju を参照してスキップするかを決めている箇所が多々出てきます。

続いて standard というサブルーチンを呼び出して標準化を行います。

      ! 2. 標準化
      call standard(no,ni,x,y,w,isd,intr,ju,g,xm,xs,ym,ys,xv,jerr)

この standard とうサブルーチンは結構大きく見えますが、切片の有無で処理を分けているため重複部分があります。 処理の内容としては:

  1. 重みの変換
  2. y と x の更新
  3. y と x の内積(共分散)を計算

となっています。

まずは重みの変換を確認してみると、重み w を「重みの総和あたりの重み」に変換し、 さらにその平方根をとったものを v として定義しています。 またその次から、先に述べたように切片の有無によって処理を分けています。

      subroutine standard(no,ni,x,y,w,isd,intr,ju,g,xm,xs,ym,ys,xv,jerr)
      implicit double precision(a-h,o-z)                                
      double precision x(no,ni),y(no),w(no),g(ni),xm(ni),xs(ni),xv(ni)  
      integer ju(ni)                                                    
      double precision, dimension (:), allocatable :: v                 
      allocate(v(1:no),stat=jerr)                                       
      if(jerr.ne.0) return
      
      ! 1. 重みの変換
      w=w/sum(w)
      v=sqrt(w) 

      ! intr は intercept なので切片が 0 であるかで判定
      ! 切片が 0 でない場合は 10141 に飛ばされる
      if(intr .ne. 0) goto 10141                                         

以降の処理ではこの vyx に対して掛け合わせるのですが、全ての観測値の重みが等しい単純なパターンを想定すると w には $1/n$ 、v にはその平方根が入ります。 例えば観測値の数が 100 であれば $w = 1/100 = 0.01$ 、$v = sqrt(1/100) = 0.1$ となります。

ではこのような wv を使って何をやっているかというと、 y に対しては:

  1. yv を乗じたものを新たに y とする
  2. その y内積(二乗和)から vy内積の二乗を減じ、平方根をとる(ys
  3. yys で割る

ということをしています。

      ! 2. y と x の更新
      ! 以下のセクションでは y と x それぞれについて観測値の重みを使って色々と調整する
      ! まずは y      
      ym = 0.0
      y  = v*y 
      ys = sqrt(dot_product(y,y)-dot_product(v,y)**2)
      y  = y/ys 

ただこの説明だけでは意味が分からないと思いますので少し式を整理してみましょう。 もとの y および w を $y0$ 、 $w0$ とおくと、


y_{i} = v_{i} * y0_{i} = \sqrt{w_{i}} * y0 = \sqrt{\frac{w0_{i}}{\sum{w0_{i}}}}y0_{i}

となります。 また ys の二乗(平方根を取る前) $ (ys)^{2} $ は:


 (ys)^{2} = \sum{y_{i}^{2}} - (\sum{v_{i}y_{i}})^2 = \sum{\frac{w_{i}}{\sum{w_{i}}}y0_{i}^{2}} - (\sum{\frac{w_{i}}{\sum{w_{i}}}}y0_{i})^{2}

と書けます。 ここで $w$ は観測値に対する重み $w0$ をその総和で除した形(単純なパターンでは $\frac{1}{n}$ )となっていることを思い出すと、これを乗じたものの総和は重み付き平均となります。 そうすると右辺の第一項はもともとの y ($y0$)の二乗の重み付き平均、第二項は重み付き平均の二乗が得られていることがわかります。 二乗の平均から平均の二乗を引いたものと言えば分散ですので、その平方根をとった ys は $y0$ の重み付き標準偏差を得ているようです (ところで $w$、$w0$ は添字 $i$ を付けるべきですが、はてなブログLaTeX がなぜか崩れるので省略しています)。

実際にサンプルデータで計算してみましょう。 まずは以下のような簡単なデータで二乗の平均から平均の二乗を引いたものが分散になることを確認します。

# 適当なデータ
a <- c(5, 5, 6, 7, 9)

# 一般的な分散の計算
mean((a - mean(a))^2)
# 二乗の平均から平均の二乗を引いてみる
mean(a^2) - mean(a)^2
# R の var を使う
var(a) * 4 / 5
[1] 2.24
[1] 2.24
[1] 2.24

上の例ではいずれも同じ値を返していることがわかります。 なお var を使った計算では不偏分散ではなく標本分散に修正しました。

続いて先の計算にしたがった場合に、やはり同じように分散・標準偏差が得られるかを見てみます。

set.seed(123)
n <- 10
y0 <- rnorm(n)
w0 <- rep(1, n)

w <- w0/sum(w0)
v <- sqrt(w)

y <- v*y0
ys <- sqrt(y %*% y - (v %*% y)^2)
y_new <- y/ys[1]
> mean((y0 - mean(y0))^2) # 一般的な分散の計算
[1] 0.8187336
> mean(y0^2) - mean(y0)^2 # 二乗の平均から平均の二乗を引く
[1] 0.8187336
> var(y0) * (n-1) / n # R の var を使って
[1] 0.8187336
> (ys^2)[1] # 計算した値
[1] 0.8187336

$(ys)^{2}$ が $y0$ の分散になっていることが確認できますね。 ということで、先ほどの処理では wv を使ってもともとの y の重み付き標準偏差を計算し、その値で重み付きの y を除していることがわかりました。 このサブルーチンの名前が standard なので当然ですが、標準化をしているようです。

x についても基本的に同様の処理を行っており、v を使って重み付き標準偏差を計算・標準化をしています。 ただし最後に重み付き平均の二乗 / 分散 に 1 を加算したものを xv に格納しており、これを x の分散としているようなのですが、これが良くわかりませんでした。

ちなみに ju は先ほど説明したように各変数に異なる数値・バラつきがあるかを示すもので、バラつきがなければさっさとループを抜けて次の変数に移っていることがわかります。

      ! x
      do 10151 j=1,ni ! ni は nvars
      if(ju(j).eq.0)goto 10151
      xm(j) = 0.0 
      x(:,j) = v*x(:,j) ! x にも重みを乗じる
      xv(j) = dot_product(x(:,j),x(:,j)) ! x の二乗の重み付き平均

      ! isd は標準化するかの指定で、標準化する場合は 1 が入っており 10171 に飛ばされない
      if(isd .eq. 0) goto 10171 
      xbq = dot_product(v, x(:,j))**2 ! x の重み付き平均の二乗
      vc = xv(j)-xbq ! 重み付き分散
      xs(j) = sqrt(vc) ! 重み付き標準偏差。 ys と対応している。
      x(:,j) = x(:,j)/xs(j) ! 標準偏差で割って標準化。 y/ys と対応している。
      
      ! これはよくわからない
      xv(j) = 1.0 + xbq/vc ! 重み付き平均の二乗 / 分散 に 1 を加算
      goto 10181 
10171 continue
      xs(j)=1.0
10181 continue
      continue
10151 continue
      continue
      goto 10191

切片が 0 でない場合はこちらにきます(基本はこっち)が、処理は上記と大体同じです。 yx ともに値を更新する前に重み付き平均を引いているところが違う点ですね。

      ! 切片が 0 でない場合ここに来る
      ! 基本はこっち
10141 continue
      ! x
      do 10201 j=1,ni
      if(ju(j).eq.0)goto 10201 
      xm(j) = dot_product(w,x(:,j)) ! x の重み付き平均
      x(:,j) = v*(x(:,j)-xm(j))  ! 重み付き平均を引いてから重みを乗じる
      xv(j) = dot_product(x(:,j),x(:,j)) ! 二乗の重み付き平均
      if(isd.gt.0) xs(j) = sqrt(xv(j)) ! 重み付き標準偏差
10201 continue
      continue
      if(isd .ne. 0)goto 10221
      xs = 1.0
      goto 10231
10221 continue
      do 10241 j=1,ni
      if(ju(j).eq.0)goto 10241
      x(:,j) = x(:,j)/xs(j) ! 標準化はここで実行
10241 continue
      continue
      xv=1.0
10231 continue
      continue
      ym = dot_product(w,y) ! y の重み付き平均 
      y  = v*(y-ym)          ! y から重み付き平均を引いたものに重みを乗じる
      ys = sqrt(dot_product(y,y)) ! 二乗和(分散)の平方根(SD)
      y  = y/ys ! 標準化

次の処理は共通のもので、y と x の内積を計算し、 g に格納します。 単純に yx内積を計算しているように見えますが、ここでの y


\frac{\sqrt{\frac{w_{i}}{\sum{w_{i}}}}y0_{i}}{SD(y0)}

x


\frac{\sqrt{\frac{w_{i}}{\sum{w_{i}}}}x0_{i}}{SD(x0)}

となっているので、その内積は重み付き共分散をそれぞれの標準偏差の積で除したもの、つまり重み付きの相関係数となっているはずです。

      ! 3. 内積(重み付き相関係数)を格納
10191 continue                                                          
      continue                                                          
      g = 0.0                                                             
      do 10251 j=1,ni 
      ! j 番目の変数にバラツキがあるなら g に y と x の内積(共分散)を格納する
      ! ただしこの時点での y と x はそれぞれ標準偏差で除したものとなっている
      if(ju(j).ne.0) g(j) = dot_product(y, x(:,j))                          
10251 continue 
      continue 
      deallocate(v) 
      return
      end 

先のサンプルデータで確かめてみましょう。 重みが全て等しいという単純なパターンでは、更新された yx内積相関係数になっていることが確認できます。

set.seed(123)
n <- 10
y0 <- rnorm(n)
x0 <- rnorm(n)
w0 <- rep(1, n)

w <- w0/sum(w0)
v <- sqrt(w)

y <- v*(y0 - (w %*% y0)[1])
ys <- sqrt(y %*% y)
y_new <- y/ys[1]

x <- v*(x0 - (w %*% x0)[1])
xs <- sqrt(x %*% x)
x_new <- x/xs[1]
> (y_new %*% x_new)[1] # 内積
[1] 0.5776151
> cor(y_new, x_new) # 更新後の y と xの相関係数
[1] 0.5776151
> cor(y0, x0) # 元の値の相関係数
[1] 0.5776151

一方重みが観測値によって異なる場合はというと、これは近い値になるものの完全に一致はしませんでした(でもこれなんでだろう、一致するような気がするんだけど)。

set.seed(123)
n <- 10
y0 <- rnorm(n)
x0 <- rnorm(n)
w0 <- rep(1, n) - 0.5 * ifelse(runif(n) > 0.8, 1, 0) # 一部のデータに対して重みを小さくしている

w <- w0/sum(w0)
v <- sqrt(w)

y <- v*(y0 - (w %*% y0)[1])
ys <- sqrt(y %*% y)
y_new <- y/ys[1]

x <- v*(x0 - (w %*% x0)[1])
xs <- sqrt(x %*% x)
x_new <- x/xs[1]
> (y_new %*% x_new)[1]
[1] 0.5687947
> cor(y_new, x_new)
[1] 0.5687133
> cor(y0, x0)
[1] 0.5776151

ところで重み調整後の yx内積相関係数と近似(一致?)するなら、個別のデータのペアが相関に対してどのような影響を持っているかを評価できるのではないでしょうか。

内積ではなく各ペアの掛け算語の値を見てみると、6 番目と 8 番目の値が高い値を示していることがわかります。 このデータの重み付き相関係数0.568 ぐらいだったので、この 2 つの観測値の影響が大きそうです。

> cbind(1:n, y_new * x_new)
      [,1]          [,2]
 [1,]    1 -0.0744142551
 [2,]    2 -0.0043928887
 [3,]    3  0.0261049036
 [4,]    4  0.0004833048
 [5,]    5 -0.0025033504
 [6,]    6  0.2852904868
 [7,]    7  0.0104270429
 [8,]    8  0.3466035906
 [9,]    9 -0.0413263433
[10,]   10  0.0225221645

実際にデータを見てみると、 6 番と 8 番のデータは他の観測値と比べて関連性が強そうに見えます。

> cbind(y_new, x_new)
             y_new       x_new
 [1,] -0.233569767  0.31859541
 [2,] -0.117117800  0.03750829
 [3,]  0.513582873  0.05082900
 [4,] -0.011106121 -0.04351698
 [5,]  0.009617489 -0.26029149
 [6,]  0.568708951  0.50164585
 [7,]  0.126538479  0.08240215
 [8,] -0.481982832 -0.71912020
 [9,] -0.278126098  0.14858851
[10,] -0.136535493 -0.16495465

6 と 8 番目のデータを塗り分けてみるとわかりやすいですね。

cols <- c(1, 1, 1, 1, 1, 3, 1, 3, 1, 1) + 1
plot(y ~ x, col = cols, pch = 16)

f:id:ushi-goroshi:20200824113309p:plain

以上で yx について標準化が終わったのでstandard から elnet に帰ってくると今度は回帰係数の上限・下限についても標準化を行います。 また flmin が 1 以上の場合は vlam を更新するのですが、 flminlambda が指定された場合に 1 が入り、そうでなければ $[0, 1)$ の実数が期待されるパラメータでした。 なので lambba が指定された場合(= flmin が 1 のとき)に vlamy の重み付き標準偏差で調整される事になります。 この vlam は後続の処理(フィッティング)では ulam として渡されるものですが、ulamlambda の指定がなければ 1 、指定があればその降順となるものでした。 要するに lambda の大きさについても標準化するよ、という事のようですね。

      ! jerr に 0 でない値が入っていると return
      if(jerr.ne.0) return

      ! cl は glmnet で cl = rbind(lower.limits, upper.limits) と定義される
      ! 回帰係数の上限・加減
      cl=cl/ys

      ! 標準化の指定が 0 であれば以下はスキップ                     
      if(isd .le. 0) goto 10091
      
      ! 説明変数ごとに標準偏差を乗じる
      do 10101 j=1,ni
      cl(:,j)=cl(:,j)*xs(j)
10101 continue                                                          
      continue                                                          
10091 continue                                                          
      
      ! flmin は glmnet のなかで flmin = as.double(lambda.min.ratio) で定義される
      ! ここで lambda.min.ratio = ifelse(nobs < nvars, 0.01, 1e-04)
      if(flmin.ge.1.0) vlam=ulam/ys 

ではフィッティングです。 ここで呼ばれる elnet1 こそが {glmnet} の本体となり、回帰係数の計算はここで行われます。 この中ではもうサブルーチンはほとんど呼ばれず、初期パラメータを取ってくるものとプログレスバーを表示するためのものだけです。 ようやくたどり着きました、今回も長かったですね。

      ! 3. フィッティング
      ! 本体である elnet1 の呼び出し
      call elnet1(parm,ni,ju,vp,cl,g,no,ne,nx,x,nlam,flmin,vlam,thr,maxi,xv,  lmu,ca,ia,nin,rsq,alm,nlp,jerr)

このサブルーチンは量はそこそこ(180行程度)なのですが、ループが込み入っていて紹介が長くなるので今回はここまでです。 また次回。

glmnetをもう少し理解したい②

前回の記事では glmnet の中身を確認し、引数の family によって呼び出す関数を変えていることがわかりました。 今回はそのなかでも gaussian が指定された場合の関数である elnet を見ていきましょう。 なお前回の記事はこちらです。

ushi-goroshi.hatenablog.com

elnet の実装

それでは早速 elnet という関数を見ていきましょう。 ちなみにここでの elnet はコンソールで elnet と打っても表示されませんが、C や Fortran で書かれたものではなくて単に glmnet からエクスポートされていない関数なので glmnet:::elnet で中身を見ることができます。

この関数はそれほど長くないのでいきなり内容の確認に入りますが、他の多くの関数同様に elnet でも最初はパラメータの受け取り・確認を行います。 下のブロックでは反復回数( maxit )、観測値の重み( weights )を受け取った後、 type.gaussian の指定内容によって ka というパラメータに格納する値を変えています。

function (x, is.sparse, ix, jx, y, weights, offset, type.gaussian = c("covariance", 
    "naive"), alpha, nobs, nvars, jd, vp, cl, ne, nx, nlam, flmin, 
    ulam, thresh, isd, intr, vnames, maxit) 
{
    # 1. パラメータの受け取り
    ### maxit
    maxit = as.integer(maxit)
    ### weights
    weights = as.double(weights)
    ### type.gaussian
    type.gaussian = match.arg(type.gaussian)
    ka = as.integer(switch(type.gaussian, covariance = 1, naive = 2, 
        ))

ka はさらに先の処理で elnetuelnetn という2つのサブルーチンのどちらを呼ぶかを決めていますので、 type.gaussian の指定に合わせてサブルーチンを変更しているということですね。

以下では y および offset (存在する場合)を double に変換しています。 また y の重み付き平均を使って Null Deviance (残差逸脱度)を計算しています。

### y の storage.mode
storage.mode(y) = "double"
### offset
if (is.null(offset)) {
  is.offset = FALSE
}
else {
  storage.mode(offset) = "double"
  is.offset = TRUE
  y = y - offset
}
### 重み付き平均
ybar = weighted.mean(y, weights)
### Null Deviance(帰無モデルの残差逸脱度)
nulldev = sum(weights * (y - ybar)^2)
if (nulldev == 0) 
  stop("y is constant; gaussian glmnet fails at standardization step")

次のブロックで早速フィッティングに入ります。 is.sparse が指定されているか否かで spelnetelnet のどちらが呼ばれるかが決まりますが、引数の違いとしては spelnet において xas.double とされており、 ixjx (いずれも疎行列において非ゼロの要素の座標を特定するための数値)が追加されています。

# 2. フィッティング
## 疎行列であるかで関数を変える
fit = if (is.sparse) 
  .Fortran("spelnet", ka, parm = alpha, nobs, nvars, x, 
           # 疎行列である場合、以下の ix, jx が引数として追加される
           # ix, jx は疎行列における非ゼロの要素の累積個数と行番号
           ix, jx, 
           y, weights, jd, vp, cl, ne, nx, nlam, flmin, 
           ulam, thresh, isd, intr, maxit, lmu = integer(1), 
           a0 = double(nlam), ca = double(nx * nlam), ia = integer(nx), 
           nin = integer(nlam), rsq = double(nlam), alm = double(nlam), 
           nlp = integer(1), jerr = integer(1), PACKAGE = "glmnet")
else .Fortran("elnet", ka, parm = alpha, nobs, nvars, as.double(x), 
              y, weights, jd, vp, cl, ne, nx, nlam, flmin, 
              ulam, thresh, isd, intr, maxit, lmu = integer(1), 
              a0 = double(nlam), ca = double(nx * nlam), ia = integer(nx), 
              nin = integer(nlam),  rsq = double(nlam), alm = double(nlam), 
              nlp = integer(1), jerr = integer(1), PACKAGE = "glmnet")
# nx は 非ゼロの変数の個数
# nlam は検証する lambda の個数
# なので ca は変数の数 * lambda の数

処理を抜けたあとは、エラーをチェックした上で必要なパラメータを取得します。

# 3. 後処理
## エラーチェック
if (fit$jerr != 0) {
  errmsg = jerr(fit$jerr, maxit, pmax = nx, family = "gaussian")
  if (errmsg$fatal) 
    stop(errmsg$msg, call. = FALSE)
  else warning(errmsg$msg, call. = FALSE)
}
## パラメータ(切片、回帰係数、自由度、次元、lambda)を取ってくる
outlist = getcoef(fit, nvars, nx, vnames)

## パラメータ(xxxxxxxxxx)を取ってきて outlist に結合する
dev = fit$rsq[seq(fit$lmu)]

outlist = c(outlist, list(dev.ratio = dev, nulldev = nulldev, 
                          npasses = fit$nlp, jerr = fit$jerr, offset = is.offset))
## elnet クラスを付与する
class(outlist) = "elnet"
outlist
}

それでは次に elnet の本体である elnet(ややこしいですね) の中身を見ていきましょう。

elnet(二度目)の実装

上記のフィッティングのセクションで elnet.Fortran("elnet") として呼ばれていました。これまで glmGAM で見てきたときと同じように、 glmnet でもやはり fortran に行き着くようですね。

と言ってもここではまだ関数自体は大きくなく、下のように(コメント抜きで)30行程度で書かれています。

subroutine elnet(ka,parm,no,ni,x,y,w,jd,vp,cl,ne,nx,nlam,  flmin,u
                 *lam,thr,isd,intr,maxit,  lmu,a0,ca,ia,nin,rsq,alm,nlp,jerr)
implicit double precision(a-h,o-z)                                
double precision x(no,ni),y(no),w(no),vp(ni),ca(nx,nlam),cl(2,ni) 
double precision ulam(nlam),a0(nlam),rsq(nlam),alm(nlam)          
integer jd(*),ia(nx),nin(nlam)                                    
double precision, dimension (:), allocatable :: vq;               

! vp が 0.0 だった場合には jerr = 100000 として return してしまう
if(maxval(vp) .gt. 0.0)goto 10021                                 
jerr=10000                                                        
return                                                            
10021 continue                                                          
allocate(vq(1:ni),stat=jerr)

! ここでも jerr に 0 以外の数値が入っていたら return してしまう
if(jerr.ne.0) return                                              

! vp の値によって vq を生成
! デフォルトは 1
! ni は nvars で変数の数なので、 vq にはデフォルトでは変数の数が入る
! でもなんで sum(vq) なんだろ
vq=max(0d0,vp)                                                    
vq=vq*ni/sum(vq)

! elnetu か elnetn のどちらを呼ぶかは ka .ne. 1 であるかで判断している
! 1 でなければ elnetn 、 1 なら elnetu
if(ka .ne. 1)goto 10041                                           
call elnetu  (parm,no,ni,x,y,w,jd,vq,cl,ne,nx,nlam,flmin,ulam,thr,
              *isd,intr,maxit,  lmu,a0,ca,ia,nin,rsq,alm,nlp,jerr)
goto 10051                                                        
10041 continue                                                          
call elnetn (parm,no,ni,x,y,w,jd,vq,cl,ne,nx,nlam,flmin,ulam,thr,i
             *sd,intr,maxit,  lmu,a0,ca,ia,nin,rsq,alm,nlp,jerr)
10051 continue                                                          
continue                                                          
deallocate(vq)                                                    
return                                                            
end                                                               

goto を多用していますね。。。 変数宣言以下で気になるところとしては、 vp が 0 だったときの挙動と、 elnetu を呼ぶところでしょうか。

vp は前回の記事で確認した通り、 glmnet のなかで vp = as.double(penalty.factor) として定義されています。 この penalty.factor はデフォルトでは 1 が入りますので基本的には goto 10021 で飛ばされてしまいます。 このセクションで引っかかるのは明示的に penalty.factor に 0 を指定した場合ですね。

! vp は各変数に対する罰則の重み(デフォルトは 1) が入ったベクトル
! vp = as.double(penalty.factor)
! jerr の数値で後続の処理で出力するエラーメッセージが決まる
if(maxval(vp) .gt. 0.0)goto 10021                                 
jerr=10000                                                        
return                                                            
10021 continue                                                          
allocate(vq(1:ni),stat=jerr)

では penalty.factor に 0 を指定した場合はどうなるかと言うと、 jerr に 10000 が入力されて return されます。 この jerr は先ほど確認した後処理において errmsg = jerr(fit$jerr, maxit, pmax = nx, family = "gaussian") としてエラーメッセージに変換されるのでした。 またこの jerr という関数は glmnet で定義されていますので、

> glmnet:::jerr
function (n, maxit, pmax, family) 
{
  if (n == 0) 
    list(n = 0, fatal = FALSE, msg = "")
  else {
    errlist = switch(family, gaussian = jerr.elnet(n, maxit, 
                                                   pmax), binomial = jerr.lognet(n, maxit, pmax), multinomial = jerr.lognet(n, 
                                                                                                                            maxit, pmax), poisson = jerr.fishnet(n, maxit, pmax), 
                     cox = jerr.coxnet(n, maxit, pmax), mrelnet = jerr.mrelnet(n, 
                                                                               maxit, pmax))
    names(errlist) = c("n", "fatal", "msg")
    errlist$msg = paste("from glmnet Fortran code (error code ", 
                        n, "); ", errlist$msg, sep = "")
    errlist
  }
}

として取り出せます。 関数をみてみると、 errlistswitch(family, ~) で更に異なる関数を呼び出し、その結果を格納しているようです。 そのため更に jerr.elnet を確認すると

> glmnet:::jerr.elnet
function (n, maxit, pmax) 
{
  if (n > 0) {
    if (n < 7777) 
      msg = "Memory allocation error; contact package maintainer"
    else if (n == 7777) 
      msg = "All used predictors have zero variance"
    else if (n == 10000) 
      msg = "All penalty factors are <= 0"
    else msg = "Unknown error"
    list(n = n, fatal = TRUE, msg = msg)
  }
  else if (n < 0) {
    if (n > -10000) 
      msg = paste("Convergence for ", -n, "th lambda value not reached after maxit=", 
                  maxit, " iterations; solutions for larger lambdas returned", 
                  sep = "")
    if (n < -10000) 
      msg = paste("Number of nonzero coefficients along the path exceeds pmax=", 
                  pmax, " at ", -n - 10000, "th lambda value; solutions for larger lambdas returned", 
                  sep = "")
    list(n = n, fatal = FALSE, msg = msg)
  }
}

else if (n == 10000) msg = "All penalty factors are <= 0" と、罰則項が 0 であることを教えてくれていますね。

さて続いて elnetu の呼びだしを確認すると、elnetuelnetn のいずれを呼ぶかは ka で決まっています。 先ほど少し触れた通り、 kaka = as.integer(switch(type.gaussian, covariance = 1, naive = 2, )) で定義されています。 また type.gaussian は glmnet の引数であり、type.gaussian = ifelse(nvars < 500, "covariance", "naive") と定義されています。 変数の数が 500 未満であれば covarinace となり、 ka には 1 が引き渡されるので if(ka .ne. 1) には該当せず、したがって elnetu が呼ばれることになるようですね。

! elnetu か elnetn のどちらを呼ぶかは ka .ne. 1 であるかで判断している
! 1 でなければ elnetn 、 1 なら elnetu
! ka は elnet の第一引数
! ka = as.integer(switch(type.gaussian, covariance = 1, naive = 2, ))
! この covariance / naive は変数の数で決まる
! type.gaussian = ifelse(nvars < 500, "covariance", "naive")
if(ka .ne. 1)goto 10041                                           
call elnetu  (parm,no,ni,x,y,w,jd,vq,cl,ne,nx,nlam,flmin,ulam,thr,
              *isd,intr,maxit,  lmu,a0,ca,ia,nin,rsq,alm,nlp,jerr)
goto 10051                                                        
10041 continue                                                          
call elnetn (parm,no,ni,x,y,w,jd,vq,cl,ne,nx,nlam,flmin,ulam,thr,i
             *sd,intr,maxit,  lmu,a0,ca,ia,nin,rsq,alm,nlp,jerr)

次回はこの elnetu を見てみましょう。

glmnetをもう少し理解したい①

久しぶりの更新です(いつも言っています)。

背景

データサイエンス入門シリーズの「スパース回帰分析とパターン認識」を読んでいたら大変面白かったので、いつものように glmnet の中身を見てみることにしました。 なお私は業務でLasso/Ridgeを使った経験があまりないため理解が間違っているかもしれませんが、その点あらかじめご了承ください。

こちらの本です。良い本です。

glmnet の実行結果

前回の GAM の時と同様に、まずは glmnet でどのような結果を得ることができるのか確認してみましょう。「スパース回帰分析とパターン認識」(以下、教科書)P12 コード1.2を(少し改変して)実行してみます。 なおこれらのコードはこちらからダウンロードすることもできます。 環境は以下のような感じです。

> 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     grid_3.6.0      lattice_0.20-38
library(glmnet)
library(plotmo)

x <- scale(LifeCycleSavings[, 2:5])
y <- LifeCycleSavings[, 1] - mean(LifeCycleSavings[, 1])

lasso <- glmnet(x, y, family = "gaussian", alpha = 1) # alpha = 1 で lasso
ridge <- glmnet(x, y, family = "gaussian", alpha = 0) # alpha = 0 で ridge

## directoryは適当に指定
png("./Image/glmnet_dive_01_01.png", width = 600, height = 400)
plot_glmnet(lasso, xvar = "lambda", label = TRUE)
dev.off()
png("./Image/glmnet_dive_01_02.png", width = 600, height = 400)
plot_glmnet(ridge, xvar = "lambda", label = TRUE)
dev.off()

f:id:ushi-goroshi:20200402204310p:plain

f:id:ushi-goroshi:20200402204326p:plain

結果の解釈などについて詳しくは教科書を見て頂くとして、 glmnet は目的関数に回帰係数の規模に応じた罰則を設けることで、回帰係数を0に向かって縮小させながらフィッティングを行います。 またグラフのように罰則の大きさを色々と動かすことで各変数への回帰係数がどのように変化するかを評価することができます。 このグラフでは左から右に向かって罰則が強くなりますが、それにつれてLasso/Ridgeの両方とも回帰係数が0に向かって小さくなっている(縮小している)ことがわかります。

なお Lasso では回帰係数が0に収束している一方、 Ridge では微小ながら最後まで係数が0とならずに残っていることがわかりますが(グラフ上部の Degrees of Freedom が 4 のままとなっています)、 Lasso のように一部の回帰係数を正確に 0 と推定することが可能な手法をスパース推定と呼びます。

glmnet の実装

それでは glmnet という関数がどのように実装されているのか見ていきましょう。 まずはいつものように全体を眺め、見通しをよくします。

function (x, y, family = c("gaussian", "binomial", "poisson", 
                           "multinomial", "cox", "mgaussian"), weights, offset = NULL, 
          alpha = 1, nlambda = 100, lambda.min.ratio = ifelse(nobs < 
                                                                nvars, 0.01, 1e-04), lambda = NULL, standardize = TRUE, 
          intercept = TRUE, thresh = 1e-07, dfmax = nvars + 1, pmax = min(dfmax * 
                                                                            2 + 20, nvars), exclude, penalty.factor = rep(1, nvars), 
          lower.limits = -Inf, upper.limits = Inf, maxit = 1e+05, type.gaussian = ifelse(nvars < 
                                                                                           500, "covariance", "naive"), type.logistic = c("Newton", 
                                                                                                                                          "modified.Newton"), standardize.response = FALSE, type.multinomial = c("ungrouped", 
                                                                                                                                                                                                                 "grouped"), relax = FALSE, trace.it = 0, ...) 
{
  
  ### 1. パラメータの設定、前処理、エラーチェック
  family = match.arg(family)
  if (alpha > 1) {
    warning("alpha >1; set to 1")
    alpha = 1
  }
  if (alpha < 0) {
    warning("alpha<0; set to 0")
    alpha = 0
  }
  alpha = as.double(alpha)
  this.call = match.call()
  nlam = as.integer(nlambda)
  y = drop(y)
  np = dim(x)
  if (is.null(np) | (np[2] <= 1)) 
    stop("x should be a matrix with 2 or more columns")
  nobs = as.integer(np[1])
  if (missing(weights)) 
    weights = rep(1, nobs)
  else if (length(weights) != nobs) 
    stop(paste("number of elements in weights (", length(weights), 
               ") not equal to the number of rows of x (", nobs, 
               ")", sep = ""))
  nvars = as.integer(np[2])
  dimy = dim(y)
  nrowy = ifelse(is.null(dimy), length(y), dimy[1])
  if (nrowy != nobs) 
    stop(paste("number of observations in y (", nrowy, ") not equal to the number of rows of x (", 
               nobs, ")", sep = ""))
  vnames = colnames(x)
  if (is.null(vnames)) 
    vnames = paste("V", seq(nvars), sep = "")
  ne = as.integer(dfmax)
  nx = as.integer(pmax)
  if (missing(exclude)) 
    exclude = integer(0)
  if (any(penalty.factor == Inf)) {
    exclude = c(exclude, seq(nvars)[penalty.factor == Inf])
    exclude = sort(unique(exclude))
  }
  if (length(exclude) > 0) {
    jd = match(exclude, seq(nvars), 0)
    if (!all(jd > 0)) 
      stop("Some excluded variables out of range")
    penalty.factor[jd] = 1
    jd = as.integer(c(length(jd), jd))
  }
  else jd = as.integer(0)
  vp = as.double(penalty.factor)
  internal.parms = glmnet.control()
  if (internal.parms$itrace) 
    trace.it = 1
  else {
    if (trace.it) {
      glmnet.control(itrace = 1)
      on.exit(glmnet.control(itrace = 0))
    }
  }
  if (any(lower.limits > 0)) {
    stop("Lower limits should be non-positive")
  }
  if (any(upper.limits < 0)) {
    stop("Upper limits should be non-negative")
  }
  lower.limits[lower.limits == -Inf] = -internal.parms$big
  upper.limits[upper.limits == Inf] = internal.parms$big
  if (length(lower.limits) < nvars) {
    if (length(lower.limits) == 1) 
      lower.limits = rep(lower.limits, nvars)
    else stop("Require length 1 or nvars lower.limits")
  }
  else lower.limits = lower.limits[seq(nvars)]
  if (length(upper.limits) < nvars) {
    if (length(upper.limits) == 1) 
      upper.limits = rep(upper.limits, nvars)
    else stop("Require length 1 or nvars upper.limits")
  }
  else upper.limits = upper.limits[seq(nvars)]
  cl = rbind(lower.limits, upper.limits)
  if (any(cl == 0)) {
    fdev = glmnet.control()$fdev
    if (fdev != 0) {
      glmnet.control(fdev = 0)
      on.exit(glmnet.control(fdev = fdev))
    }
  }
  storage.mode(cl) = "double"
  isd = as.integer(standardize)
  intr = as.integer(intercept)
  if (!missing(intercept) && family == "cox") 
    warning("Cox model has no intercept")
  jsd = as.integer(standardize.response)
  thresh = as.double(thresh)
  if (is.null(lambda)) {
    if (lambda.min.ratio >= 1) 
      stop("lambda.min.ratio should be less than 1")
    flmin = as.double(lambda.min.ratio)
    ulam = double(1)
  }
  else {
    flmin = as.double(1)
    if (any(lambda < 0)) 
      stop("lambdas should be non-negative")
    ulam = as.double(rev(sort(lambda)))
    nlam = as.integer(length(lambda))
  }
  is.sparse = FALSE
  ix = jx = NULL
  if (inherits(x, "sparseMatrix")) {
    is.sparse = TRUE
    x = as(x, "CsparseMatrix")
    x = as(x, "dgCMatrix")
    ix = as.integer(x@p + 1)
    jx = as.integer(x@i + 1)
    x = as.double(x@x)
  }
  if (trace.it) {
    if (relax) 
      cat("Training Fit\n")
    pb <- createPB(min = 0, max = nlam, initial = 0, style = 3)
  }
  kopt = switch(match.arg(type.logistic), Newton = 0, modified.Newton = 1)
  if (family == "multinomial") {
    type.multinomial = match.arg(type.multinomial)
    if (type.multinomial == "grouped") 
      kopt = 2
  }
  kopt = as.integer(kopt)
  
  ### 2. フィッティング
  fit = switch(family, gaussian = elnet(x, is.sparse, ix, jx, 
                                        y, weights, offset, type.gaussian, alpha, nobs, nvars, 
                                        jd, vp, cl, ne, nx, nlam, flmin, ulam, thresh, isd, intr, 
                                        vnames, maxit), poisson = fishnet(x, is.sparse, ix, jx, 
                                                                          y, weights, offset, alpha, nobs, nvars, jd, vp, cl, ne, 
                                                                          nx, nlam, flmin, ulam, thresh, isd, intr, vnames, maxit), 
               binomial = lognet(x, is.sparse, ix, jx, y, weights, offset, 
                                 alpha, nobs, nvars, jd, vp, cl, ne, nx, nlam, flmin, 
                                 ulam, thresh, isd, intr, vnames, maxit, kopt, family), 
               multinomial = lognet(x, is.sparse, ix, jx, y, weights, 
                                    offset, alpha, nobs, nvars, jd, vp, cl, ne, nx, nlam, 
                                    flmin, ulam, thresh, isd, intr, vnames, maxit, kopt, 
                                    family), cox = coxnet(x, is.sparse, ix, jx, y, weights, 
                                                          offset, alpha, nobs, nvars, jd, vp, cl, ne, nx, nlam, 
                                                          flmin, ulam, thresh, isd, vnames, maxit), mgaussian = mrelnet(x, 
                                                                                                                        is.sparse, ix, jx, y, weights, offset, alpha, nobs, 
                                                                                                                        nvars, jd, vp, cl, ne, nx, nlam, flmin, ulam, thresh, 
                                                                                                                        isd, jsd, intr, vnames, maxit))
  if (trace.it) {
    utils::setTxtProgressBar(pb, nlam)
    close(pb)
  }
  
  ### 3. 後処理
  if (is.null(lambda)) 
    fit$lambda = fix.lam(fit$lambda)
  fit$call = this.call
  fit$nobs = nobs
  class(fit) = c(class(fit), "glmnet")
  if (relax) 
    relax.glmnet(fit, x = x, y = y, weights = weights, offset = offset, 
                 lower.limits = lower.limits, upper.limits = upper.limits, 
                 check.args = FALSE, ...)
  else fit
}

glmnet では以上のように、

  1. パラメータの設定、前処理、エラーチェック
  2. フィッティング
  3. 後処理

といったステップで処理が進んでおり、これは過去にみてきた glmgam と同様ですね。 それでは各ステップを細かく見ていきましょう。

1. パラメータの設定、前処理、エラーチェック

まずはパラメータの設定や前処理に関わる部分ですが、はじめに family の指定が問題ないかをチェックします。

## 指定したfamilyが引数としてOKかチェック
family = match.arg(family)

glmnet で使用可能な familyglm とは異なっており、Gamma / inverse.gaussian / quasi- が使えない代わりに、 multinomial / cox / mgaussian が使えるようになっています。 ここで multinomial は多項分布、mgaussian は多変量正規分布を意味するようです。

family のチェックには match.arg 関数が使われています。 この関数の挙動を理解するのは少し難しいのですが、こちらのブログが参考になります。

続いて alpha をチェックします:

## alpha
### LassoとRidgeそれぞれに対するペナルティの配分を決めるパラメータ
### glmnetにおける罰則項は以下で定義
### alphaは0~1で、1ならLasso、0ならRidgeに対応
if (alpha > 1) {
  warning("alpha >1; set to 1")
  alpha = 1
}
if (alpha < 0) {
  warning("alpha<0; set to 0")
  alpha = 0
}
alpha = as.double(alpha)

glmnet においてこの alpha は、回帰係数のL1およびL2ノルムそれぞれに対する罰則の割合をコントロールします。 より具体的には、 glmnet では罰則項は以下によって定義されます(https://cran.r-project.org/web/packages/glmnet/glmnet.pdf の P19より):


(1 − \alpha)/2||\beta||^{2}_{2} + \alpha||\beta||_{1}

冒頭のコードでは alpha = 1 または alpha = 0 としましたが、上の式から alpha = 1 のときにL2ノルムに対する罰則が消えてL1ノルムのみが残り(Lasso)、逆に alpha = 0 とするとL1ノルムに対する罰則が消えてL2ノルムが残る(Ridge)ことがわかります。 また alpha を (0, 1) とすると両者がそれぞれの割合でブレンドされます。

なお、ここでL2ノルムに対する罰則が1/2になっている理由はよくわかりませんでした。 glmnet の help で引用されているこちらの論文ではすでに $(1-\alpha)1/2||\beta||^2_2$ として定義されています。 またscikit-learnでも同様にL2ノルムに対しては0.5を乗じているようです(https://scikit-learn.org/stable/modules/generated/sklearn.linear_model.ElasticNet.html)。 誰か理由を教えてください。

続いて match.call() を用いて引数の指定を正式なものに直します:

## match.call
this.call = match.call()

これだけだと何を言っているかちょっとわからないと思いますので、以下の例で確認してみましょう:

myfun <- function(abc, def, ghi) { 
  return(abc + 2*def + 3*ghi)
}

上のように引数として abcdefghi を取る関数を定義します。 このとき R では、引数の指定がない場合には順番通りに入力されます:

> myfun(1, 2, 3)
[1] 14

一部の引数のみ指定がある場合では指定された引数だけがその通りに入力され、残りは順番通りに割り当てられるようです。

> myfun(def = 3, 4, 5)
[1] 25

ところでこの引数の指定は、一意に決まれば指定は省略することができます:

> myfun(d = 3, 4, 5)
[1] 25

一方、例えば以下のような呼び出しでは g から始まる引数が2つあるため一意に決まらず、エラーとなってしまいます。

> myfun2 <- function(abc, def, ghi, gjk) {
+   return(abc + 2*def + 3*ghi + 4*gjk)
+ }
> myfun2(g = 3, 4, 5, 6)
 myfun2(g = 3, 4, 5, 6) でエラー:  引数 1 が複数の仮引数に一致します 

では match.call を使って関数を呼び出すとどうなるかと言うと:

> match.call(myfun, call("myfun", 1, def = 3, ghi = 5))
myfun(abc = 1, def = 3, ghi = 5)

この通り、各引数に対して何を割り当てたかを得ることができます。 便利ですね。

さらに続いて、 nlambda の指定です。 ここでは $\lambda$ (罰則の大きさ)そのものではなく、検証する $\lambda$ の数(nubmer of lambda)を指定します(デフォルトは100)。

## nlambda
nlam = as.integer(nlambda)

ここからは yx および weight のチェックです:

## drop
y = drop(y)

## x
### x は2列以上持つ必要があるので、単回帰はできない様子
np = dim(x)
if (is.null(np) | (np[2] <= 1)) 
  stop("x should be a matrix with 2 or more columns")
### x のレコード数
nobs = as.integer(np[1])

### weights
### 未入力のときは 1 を与え、weights と nobs が一致しないときはエラー
if (missing(weights)) 
  weights = rep(1, nobs)
else if (length(weights) != nobs) 
  stop(paste("number of elements in weights (", length(weights), 
             ") not equal to the number of rows of x (", nobs, 
             ")", sep = ""))

### 変数の数
nvars = as.integer(np[2])

## y
dimy = dim(y)
### y のレコード数 
nrowy = ifelse(is.null(dimy), length(y), dimy[1])
### y と x でレコード数が合わないときはエラー
if (nrowy != nobs) 
  stop(paste("number of observations in y (", nrowy, ") not equal to the number of rows of x (", 
             nobs, ")", sep = ""))
## 変数名
vnames = colnames(x)
if (is.null(vnames)) 
  vnames = paste("V", seq(nvars), sep = "")

y に対する drop ですが、これは length が 1 であるような冗長な次元を落とす関数です。 続いて x の行数が weighty と合わない場合にエラーを返しています。

以下ではモデルに含める変数や非ゼロとする変数などを指定します ( nx(=pmax) の方はちょっと理解がアヤシイので help の説明を書いておきます):

## 自由度
### モデルに含まれる変数の上限を指定
### dfmax = nvars + 1
ne = as.integer(dfmax)

### 非ゼロとする変数の数の上限(?)
### Limit the maximum number of variables ever to be nonzero
### pmax = min(dfmax * 2 + 20, nvars)
nx = as.integer(pmax)

### 除外対象となる変数の指定
if (missing(exclude)) 
  exclude = integer(0)

次に変数ごとに異なるペナルティを与えるために penalty.factor を指定します。 この数値が lambda に乗じられるため、例えば特定の変数に対して penalty.factor = 0 としておけば罰則を与えないようにすることが可能となります(結果として常にモデルに採用されるようになる):

## 変数ごとに異なるペナルティを与える
### デフォルトは 1 が入る
### Inf が指定されている変数は exclude として扱われる
if (any(penalty.factor == Inf)) {
  exclude = c(exclude, seq(nvars)[penalty.factor == Inf])
  exclude = sort(unique(exclude))
}
if (length(exclude) > 0) {
  jd = match(exclude, seq(nvars), 0)
  if (!all(jd > 0)) 
    stop("Some excluded variables out of range")
  penalty.factor[jd] = 1
  jd = as.integer(c(length(jd), jd))
}
else jd = as.integer(0)
vp = as.double(penalty.factor)

これはせっかくなので実際にやってみましょう。 冒頭のコードを持ってきて、以下のように lambda を適当に設定してみます。

x <- scale(LifeCycleSavings[, 2:5])
y <- LifeCycleSavings[, 1] - mean(LifeCycleSavings[, 1])
> coef(glmnet(x, y, family = "gaussian", alpha = 1, lambda = 0.3))
5 x 1 sparse Matrix of class "dgCMatrix"
                       s0
(Intercept)  1.182354e-15
pop15       -1.691002e+00
pop75        .           
dpi          .           
ddpi         9.816514e-01

このとき、2・3番目の変数である pop75dpi は 0 と推定されてしまいました。 ここでこれらの変数の penalty.factor を 0 とすると

> coef(glmnet(x, y, family = "gaussian", alpha = 1, lambda = 0.3,
+             penalty.factor = c(1, 0, 0, 1)))
5 x 1 sparse Matrix of class "dgCMatrix"
                       s0
(Intercept)  9.523943e-16
pop15       -7.827680e-01
pop75        8.127991e-01
dpi         -1.560908e-01
ddpi         6.812498e-01

ちゃんと推定されるようになっています。 逆に pop15penalty.factor を大きくすると

> coef(glmnet(x, y, family = "gaussian", alpha = 1, lambda = 0.3,
+             penalty.factor = c(2, 0, 0, 1)))
5 x 1 sparse Matrix of class "dgCMatrix"
                      s0
(Intercept) 7.266786e-16
pop15       .           
pop75       1.374655e+00
dpi         2.586151e-02
ddpi        9.300500e-01

このようにモデルから除外されることになります。 さらに penalty.factor = Inf とすると、その変数は exclude として扱われるようになります。

続いて glmnet.control で持っているパラメータを渡します。

## 内部でデフォルトで持っているパラメータ 
internal.parms = glmnet.control()
### プログレスバーを表示する!
if (internal.parms$itrace) 
  trace.it = 1
else {
  if (trace.it) {
    glmnet.control(itrace = 1)
    on.exit(glmnet.control(itrace = 0))
  }
}

次に、回帰係数に対する上限・下限を設定します。 なお下限は non-positive 、上限は non-negative しか設定できないようですね。

## 上限・下限
### lower.limit としては非正の値のみ指定できる
if (any(lower.limits > 0)) {
  stop("Lower limits should be non-positive")
}
### upper.limtit は逆に非負の値のみ指定できる
if (any(upper.limits < 0)) {
  stop("Upper limits should be non-negative")
}
### Inf (デフォルト)になっているものについては特定の値(9.9e35)に差し替え  
lower.limits[lower.limits == -Inf] = -internal.parms$big
upper.limits[upper.limits == Inf] = internal.parms$big

### nvars との整合性チェック
if (length(lower.limits) < nvars) {
  ### lower.limits としてスカラが指定されている場合は nvars 全てに適用
  if (length(lower.limits) == 1) 
    lower.limits = rep(lower.limits, nvars)
  else stop("Require length 1 or nvars lower.limits")
}
### lower.limits が nvars よりも長い場合は前から利用する
else lower.limits = lower.limits[seq(nvars)]
### nvars との整合性チェック(lower.limits と同様)
if (length(upper.limits) < nvars) {
  if (length(upper.limits) == 1) 
    upper.limits = rep(upper.limits, nvars)
  else stop("Require length 1 or nvars upper.limits")
}
else upper.limits = upper.limits[seq(nvars)]
### 上限・下限
### coefficient limit?
cl = rbind(lower.limits, upper.limits)

### lower または upper に 0 を含む場合
### 0除算が発生するときのエラー対策?
if (any(cl == 0)) {
  ### fdev は最小となるデビアンスの変化量(割合)
  ### minimum fractional change in deviance for stopping path; factory default = 1.0e5
  fdev = glmnet.control()$fdev
  if (fdev != 0) {
    glmnet.control(fdev = 0)
    on.exit(glmnet.control(fdev = fdev)) # 関数終了時に実行される処理
  }
}
storage.mode(cl) = "double"

標準化と切片に対する指定です。 標準化の処理そのものは以降の関数の内部で実行されるため、ここでは指定のみを行います。

## 標準化
### standardize と intercept はデフォルトは TRUE なので 1 になる
isd = as.integer(standardize)
intr = as.integer(intercept)
### Cox回帰における警告
if (!missing(intercept) && family == "cox") 
  warning("Cox model has no intercept")
### standardize.response は family="mgaussian" のときに目的変数を標準化するかの指定
jsd = as.integer(standardize.response)

収束を判定する閾値を指定します。

## 収束判定
### coordinate descent における収束の閾値
thresh = as.double(thresh)

次に、 lambda に関する指定となりますが、 flmin および ulam の使われ方がよく理解できなかったため、これらの説明は省略します。 なお help にもありますが、通常は lambda には単一の値ではなく、候補となる値のベクトルを与えます。

Avoid supplying a single value for lambda (for predictions after CV use predict() instead).

## lambda
### ペナルティの大きさ
### 指定がない場合、flmin と ulam は lambda.min.ratio および 1 に指定される
### lambda.min.ratio = ifelse(nobs < nvars, 0.01, 1e-04)
if (is.null(lambda)) {
  if (lambda.min.ratio >= 1) 
    stop("lambda.min.ratio should be less than 1")
  flmin = as.double(lambda.min.ratio)
  ulam = double(1)
}
### 指定がある場合、flmin(下限?)とulam(上限?)は 1 および lambdaの降順 に指定される
else {
  flmin = as.double(1)
  if (any(lambda < 0)) 
    stop("lambdas should be non-negative")
  ulam = as.double(rev(sort(lambda)))
  nlam = as.integer(length(lambda))
}

次に疎行列の指定です。 入力 X が疎行列である場合、dgCMatrix 形式に変換されます。 ここで dgCMatrix とは列方向の志向性を持つ疎行列の形式です。

## sparse matrix 
### x が Matrix::sparseMatrix の場合は Matrix::dgCMatrix に変換する
### dgCMatrix: csc順に並び替えて(csc形式)の疎行列圧縮保管
is.sparse = FALSE
ix = jx = NULL
if (inherits(x, "sparseMatrix")) {
  is.sparse = TRUE
  x = as(x, "CsparseMatrix")
  x = as(x, "dgCMatrix")
  ### x@p は各列の非ゼロの値の個数を積み上げたものが格納されている(列数 + 1)
  ### diff(x@p + 1) すれば各列の非ゼロの値の個数がわかる
  ix = as.integer(x@p + 1)
  ### x@i は各列の非ゼロの値の行番号が格納されている(なので length(x@i) が非ゼロの値の個数と一致する)
  ### 0-index なので R のスタイルと合わせるために +1 しているのでしょう
  jx = as.integer(x@i + 1)
  ### x@x は非ゼロである値そのもののベクトル
  x = as.double(x@x)
}

ここも、せっかくなので疎行列における数値の格納方法についても見ておきましょう。 以下のように疎行列を作成します:

set.seed(1234)
i <- c(1, 5, 18)
j <- c(4, 13, 19)
n <- rnorm(3)

m <- matrix(0, 20, 20)
for (k in 1:length(n)) {
  m[i[k], j[k]] <- n[k]
}

s_m <- as(m, "dgCMatrix")

ここで s_m は行列 m を疎行列として扱ったものです。 str() で確認すると、 s_m には

  • @ i :非ゼロの要素の入っていた行番号( 0-index であることに注意)
  • @ p :各列における非ゼロの要素の個数を積み上げたもの
  • @ Dim :行列の次元
  • @ Dimnames :行列の各次元の名前
  • @ x :非ゼロの要素の数値
  • @ factors :(これはちょっとわかりませんでした)

が格納されています。

> str(s_m)
Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
  ..@ i       : int [1:3] 0 4 17
  ..@ p       : int [1:21] 0 0 0 0 1 1 1 1 1 1 ...
  ..@ Dim     : int [1:2] 20 20
  ..@ Dimnames:List of 2
  .. ..$ : NULL
  .. ..$ : NULL
  ..@ x       : num [1:3] -1.207 0.277 1.084
  ..@ factors : list()

ここで @ i には非ゼロである各要素の行番号が入るため行列 m を作ったときの行番号の指定 i に対応しますが、0-index であるため数字は1つずつ小さくなっています。

> print(i- 1)
[1]  0  4 17
> print(s_m@i)
[1]  0  4 17

ちょっとわかりにくいのが @ p で、ここには各列における非ゼロの要素の個数の累積が格納され、列数に対応します(ただし最初に 0 が追加されるため、列数 + 1 の長さになります)。 今回の例では行列の列数が 20 なので、length が 21 となります。

> length(s_m@p)
[1] 21

このベクトルには非ゼロの要素の個数の累積が入っているため、差分を取ると元の行列で非ゼロの要素が入っていた列を得ることができます。

> diff(s_m@p)
 [1] 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 1 0

列番号を指定した j と比較してみましょう:

> which(diff(s_m@p) == 1)
[1]  4 13 19
> j
[1]  4 13 19

合っていますね。 続く処理では、 ix には各列における非ゼロの要素の累積個数(+1)を 、 jx には行番号を代入しています。 また x には元の疎行列における非ゼロの要素の値そのものをベクトルとして入力しており、説明変数の行列が疎行列であった場合、この時点で行列ではなくベクトルとして扱われることになります。

次に、プログレスバーの指定です(出せるんですね)。

## プログレスバー
if (trace.it) {
  if (relax) 
    cat("Training Fit\n")
  pb <- createPB(min = 0, max = nlam, initial = 0, style = 3)
}

そして最後に最適化の手法についての指定です。 family`binomial または multinomial の場合、 glmnet の引数である type.logistic および type.multinomial が評価され、(後の工程で)それに応じて呼ばれる関数が変わります。 具体的には lognet2mlognetn および multlognetn のどれが選ばれるかが決まります。 これは別の機会に解説します(予定です)。

## 最適化の手法(ロジスティックおよび多項ロジスティックの時)
### type.logistic = c("Newton", "modified.Newton")
### Newton を指定なら 0、modified.Newton を指定なら 1 を返す
### If "Newton" then the exact hessian is used (default), while "modified.Newton" uses an upper-bound on the hessian, and can be faster.
kopt = switch(match.arg(type.logistic), Newton = 0, modified.Newton = 1)
### type.multinomial = c("ungrouped", "grouped")
### 多項ロジスティックで更にgroupedの場合は kopt は 2 となる
### If "grouped" then a grouped lasso penalty is used on the multinomial coefficients for a variable. This ensures they are all in our out together. 
### The default is "ungrouped"
if (family == "multinomial") {
  type.multinomial = match.arg(type.multinomial)
  if (type.multinomial == "grouped") 
    kopt = 2
}
kopt = as.integer(kopt)

最初の方で family のチェックに使われ、ここでも使われている match.arg ですが、せっかくなので挙動を確認しておきましょう:

### 引数に type.logistic を持つ関数を定義
myfun <- function(a = "aaa", type.logistic = c("Newton", "modified.Newton")) {
  ### 呼び出し元の関数の引数をチェックし、 Newton なら 0、modified.Newton なら 1を割り当てる
  kopt <- switch(match.arg(type.logistic), Newton = 0, modified.Newton = 1)
  kopt
}

上のような関数を定義し、以下のように呼び出すと、結果はそれぞれ 0, 0, 1 となります。

> myfun()
[1] 0
> myfun(type.logistic = "Newton")
[1] 0
> myfun(type.logistic = "modified.Newton")
[1] 1

2. フィッティング

以上でパラメータの設定や前処理が終わりましたので次はフィッティングです。 といってもここでは family に応じて呼び出す関数を変えているだけなので、詳細は一旦スキップしましょう。

# フィッティング
## family に応じてその後に呼び出す関数を変える
fit = switch(family,
             ### gaussian のときは elnet 
             gaussian = elnet(x, is.sparse, ix, jx, 
                              y, weights, offset, type.gaussian, alpha, nobs, nvars, 
                              jd, vp, cl, ne, nx, nlam, flmin, ulam, thresh, isd, intr, 
                              vnames, maxit), 
             ### poisson のときは fishnet
             poisson = fishnet(x, is.sparse, ix, jx, 
                               y, weights, offset, alpha, nobs, nvars, jd, vp, cl, ne, 
                               nx, nlam, flmin, ulam, thresh, isd, intr, vnames, maxit),
             ### binomial のときは lognet
             binomial = lognet(x, is.sparse, ix, jx, y, weights, offset, 
                               alpha, nobs, nvars, jd, vp, cl, ne, nx, nlam, flmin, 
                               ulam, thresh, isd, intr, vnames, maxit, kopt, family), 
             ### multinomial のときも lognet
             multinomial = lognet(x, is.sparse, ix, jx, y, weights, 
                                  offset, alpha, nobs, nvars, jd, vp, cl, ne, nx, nlam, 
                                  flmin, ulam, thresh, isd, intr, vnames, maxit, kopt, 
                                  family), 
             ### cox のときは coxnet
             cox = coxnet(x, is.sparse, ix, jx, y, weights, 
                          offset, alpha, nobs, nvars, jd, vp, cl, ne, nx, nlam, 
                          flmin, ulam, thresh, isd, vnames, maxit), 
             ### mgaussian のときは mrelnet
             mgaussian = mrelnet(x, 
                                 is.sparse, ix, jx, y, weights, offset, alpha, nobs, 
                                 nvars, jd, vp, cl, ne, nx, nlam, flmin, ulam, thresh, 
                                 isd, jsd, intr, vnames, maxit))
## プログレスバー
if (trace.it) {
  utils::setTxtProgressBar(pb, nlam)
  close(pb)
}

なおここでそれぞれの関数に渡されている引数を比較すると以下のようになります(一部はよくわかりませんでした):

引数 説明 elnet fishnet lognet coxnet mrelnet
x 説明変数の行列
is.sparse 疎行列であるかの指定
ix 疎行列における非ゼロの要素の累積個数
jx 疎行列における非ゼロの要素の行番号
y 目的変数の行列
weights 観測値に対する重み
offset オフセット
type.gaussian 1:covariance, 2:naïve - - - -
alpha L1とL2に対する重みの調整パラメータ
nobs レコード数
nvars 説明変数の数
jd ?
vp 各変数に対する罰則の重み(penalty.factor)
cl ?
ne モデルに含まれる変数の上限。ne = dfmax = nvars + 1
nx 非ゼロとする変数の個数の上限?
nlam lambdaの数
flmin ?
ulam ?
thresh 収束判定の閾値
isd standardizeするかの指定
jsd ? - - - -
intr 切片(Intercept)を含めるかの指定 -
vnames 変数名
maxit 反復回数の上限
kopt 最適化の手法 - - - -
family family - - - -

3. 後処理

最後に後処理です。

# 後処理
## lambda が指定されておらず fit$lambda が 3 パターン以上検証されている場合、先頭を差し替える
## glmnet::fix.lam
## function (lam) {
## if (length(lam) > 2) {
##     llam = log(lam)
##     lam[1] = exp(2 * llam[2] - llam[3])
## }
## lam
## }
if (is.null(lambda)) 
  fit$lambda = fix.lam(fit$lambda)
## call
fit$call = this.call
## レコード数
fit$nobs = nobs
## class に glmnet を追加
class(fit) = c(class(fit), "glmnet")

# リターン
## relax が TRUE の場合、解パスの各セットについて罰則なしでモデルをフィッティングする   
## If TRUE then for each active set in the path of solutions, the model is refit without any regularization. See details for more information. 
## This argument is new, and users may experience convergence issues with small datasets, especially with non-gaussian families. 
## Limiting the value of ’maxp’ can alleviate these issues in some cases.
if (relax) 
  relax.glmnet(fit, x = x, y = y, weights = weights, offset = offset, 
               lower.limits = lower.limits, upper.limits = upper.limits, 
               check.args = FALSE, ...)
else fit

この後処理で目立つ工程としては relax の部分でしょう。 ここで relax は help によると、

If relax=TRUE a duplicate sequence of models is produced, where each active set in the elastic-net path is refit without regularization. The result of this is a matching "glmnet" object which is stored on the original object in a component named "relaxed", and is part of the glmnet output.

ということで、glmnet によって変数選択された結果を用いて、罰則なしで再度フィッティングを行うオプションのようです。 これも実際にやってみるのが早いと思いますので、以下のように実行してみます:

lasso_02 <- glmnet(x, y, family = "gaussian", relax = T)

すると、先程の結果( lasso )に、 lasso_02$relaxed という結果が追加されていることがわかりますが、内容は lasso とほとんど同じです。

> str(lasso)
List of 12
 $ a0       : Named num [1:68] 6.11e-16 6.71e-16 7.26e-16 7.76e-16 8.22e-16 ...
  ..- attr(*, "names")= chr [1:68] "s0" "s1" "s2" "s3" ...
 $ beta     :Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
  .. ..@ i       : int [1:216] 0 0 0 0 0 3 0 3 0 3 ...
  .. ..@ p       : int [1:69] 0 0 1 2 3 4 6 8 10 12 ...
  .. ..@ Dim     : int [1:2] 4 68
  .. ..@ Dimnames:List of 2
  .. .. ..$ : chr [1:4] "pop15" "pop75" "dpi" "ddpi"
  .. .. ..$ : chr [1:68] "s0" "s1" "s2" "s3" ...
  .. ..@ x       : num [1:216] -0.181 -0.347 -0.497 -0.634 -0.757 ...
  .. ..@ factors : list()
 $ df       : int [1:68] 0 1 1 1 1 2 2 2 2 2 ...
 $ dim      : int [1:2] 4 68
 $ lambda   : num [1:68] 2.02 1.84 1.68 1.53 1.39 ...
 $ dev.ratio: num [1:68] 0 0.0352 0.0645 0.0888 0.1089 ...
 $ nulldev  : num 984
 $ npasses  : int 562
 $ jerr     : int 0
 $ offset   : logi FALSE
 $ call     : language glmnet(x = x, y = y, family = "gaussian", alpha = 1)
 $ nobs     : int 50
 - attr(*, "class")= chr [1:2] "elnet" "glmnet"

> str(lasso_02$relaxed)
List of 12
 $ a0       : Named num [1:68] 6.11e-16 1.29e-15 1.29e-15 1.29e-15 1.29e-15 ...
  ..- attr(*, "names")= chr [1:68] "s0" "s1" "s2" "s3" ...
 $ beta     :Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
  .. ..@ i       : int [1:216] 0 0 0 0 0 3 0 3 0 3 ...
  .. ..@ p       : int [1:69] 0 0 1 2 3 4 6 8 10 12 ...
  .. ..@ Dim     : int [1:2] 4 68
  .. ..@ Dimnames:List of 2
  .. .. ..$ : chr [1:4] "pop15" "pop75" "dpi" "ddpi"
  .. .. ..$ : chr [1:68] "s0" "s1" "s2" "s3" ...
  .. ..@ x       : num [1:216] -2.04 -2.04 -2.04 -2.04 -1.98 ...
  .. ..@ factors : list()
 $ df       : int [1:68] 0 1 1 1 1 2 2 2 2 2 ...
 $ dim      : int [1:2] 4 68
 $ lambda   : num [1:68] 2.02 1.84 1.68 1.53 1.39 ...
 $ dev.ratio: num [1:68] 0 0.208 0.208 0.208 0.208 ...
 $ nulldev  : num 984
 $ npasses  : int 562
 $ jerr     : int 0
 $ offset   : logi FALSE
 $ call     : language glmnet(x = x, y = y, family = "gaussian", relax = T)
 $ nobs     : int 50
 - attr(*, "class")= chr [1:2] "elnet" "glmnet"

ここで lasso_02$relaxed の中身を少し見てみると、例えば beta には以下のような数値が入っています。

> lasso_02$relaxed$beta[, 1:6]
4 x 6 sparse Matrix of class "dgCMatrix"
      s0        s1        s2        s3        s4        s5
pop15  . -2.040996 -2.040996 -2.040996 -2.040996 -1.980216
pop75  .  .         .         .         .         .       
dpi    .  .         .         .         .         .       
ddpi   .  .         .         .         .         1.270865

これは何かと言うと、少しずつ罰則の重みを変えたことで変数が選択された状態で通常の線形回帰を当てはめたときの回帰係数となっています。 例えば lasso_02$relaxed$beta[, 6] には、変数として選択された pop15ddpi それぞれの回帰係数が入っています。 実際に lm の結果と一致するか見てみましょう:

> coef(lm(y ~ x[, c(1, 4)]))
      (Intercept) x[, c(1, 4)]pop15  x[, c(1, 4)]ddpi 
     1.364331e-15     -1.980216e+00      1.270865e+00 

合っていますね。 ところで切片の推定値が入っている lasso_02$relaxed$a0 の値は少し異なるようです:

> lasso_02$relaxed$a0[6]
         s5 
1.28119e-15 

なんでやろか。。

もしかしたら標準化の違いかとも思いましたがそれでもないようで、この理由はわかりませんでした。

lasso_03 <- glmnet(x, y, family = "gaussian", relax = T, standardize = F)
> lasso_03$relaxed$a0[6]
         s5 
1.28119e-15 

glmnet() の実装は以上となります。 次回はフィッティングの部分で呼ばれている elnet を詳しく見ていきましょう。 なお gam のときとは違い、 glmnet では library をインストールしてもソースコードは付いてきませんでしたので、こちらを参考に fortranソースコードを取得しました。

ではまた次回。

統計数理研究所公開講座「統計の哲学を理解するために」参加記録

1/31に開催された統計数理研究所公開講座「統計の哲学を理解するために」に参加してきましたのでそのメモを共有しておきます。全体的にはエリオット・ソーバーという尤度主義者から見た頻度主義・ベイズ主義に対する批判的観点の紹介という構成で、それぞれの立場が答えようとしている問いが浮き彫りになるような内容でした。

なおこの講義は「科学と証拠」という書籍に基づいています。とても面白い本です。

科学と証拠―統計の哲学 入門―

科学と証拠―統計の哲学 入門―

また最後に、あまり時間がなくて駆け足での紹介となっていましたが、Deborah Mayoによる「誤り統計学」について簡単に紹介がありました(松王先生曰く世界初の資料、とのこと)。

Statistical Inference as Severe Testing: How to Get Beyond the Statistics Wars

Statistical Inference as Severe Testing: How to Get Beyond the Statistics Wars


日時

2020/1/31 10:00@立川

講師

北海道大学 松王政浩

エリオット・ソーバー

  • アメリカの代表的な科学哲学者
  • 尤度主義 + AICを支持
  • 著書
    • 科学と証拠
    • OCKHAM'S RAZORSの翻訳がでる予定

統計学論争は終わっていない

  • Statistical Inference as Severe Testing
    • Deborah Mayo
      • 頻度主義
  • 有意検定論争、レシピ的な統計学の食い止め

統計の哲学

  • 統計の基礎をめぐる議論の総体
  • 哲学は常に論争で成り立ち、結論には到達しない
    • ベイズ主義
      • 信念度合い
        • 主観
        • 客観
    • 頻度主義
      • ネイマン-ピアソン
        • 無限回の施行を前提とした考えた
      • フィッシャー
        • 相対頻度、推測確率
  • 統一的見解を示したものはないので各主義を学ぶしかない
    • ベイズ
      • 実態が捉えにくい
      • 意思決定中心、科学的仮説の確証中心(科学哲学)、ベイズ統計の流儀
      • バーガー&ウォルパート(1988)、ゲルマン(一連)
    • 頻度主義
      • Mayoのものは包括的だがわかりにくい
    • 尤度主義
      • コンパクトでわかりやすい
      • ソーバーやロイヤルは話が短く論点が明快だが排他主義

3つの主義

  • ロイヤルの3つの問い
    • 証拠をもとに何がわかるか → 尤度主義
      • 静的現在主義、外在主義
      • データが直接示している「仮説に関する情報(証拠)」をキャッチせよ
      • いったんデータが得られたら他のデータの可能性は捨象せよ
    • 証拠をもとに何を信じるか → ベイズ主義
      • 動的現在主義、内在主義
      • 他のデータの可能性は一切捨象せよ
      • データが「(可能な)仮説に関する情報」をどう変化させるかキャッチせよ
    • 証拠をもとに何をするべきか → 頻度主義
      • 反事実主義、規約主義
      • この場で生じていることだけでなく、生じた可能性がある事柄もすべて併せて仮説について判断せよ
      • 低い確率を棄却のサインとみなそう、同じルールを何度も適用すれば誤りの可能性が低いのだから今回もそのルールを適用する

尤度「原理」について

  • ソーバーの議論の核となる2つの原理
    • 穏当な原理
      • もしEが真であると知ることによって命題Pを棄却することが正当化され、かつ、この情報を得てはじめてPの棄却が正当化されたのであれば、EはPに反する証拠とされねばならない
    • 全証拠の原理
      • 実験によって得られたデータはすべて証拠の判断に用いられなければならない
  • 尤度の法則とは異なる
    • 尤度主義、ベイズ主義の両方に共通する源泉
  • 尤度原理(LP)
    • xが観測されたあと、θについて推論(決定)する際に、実験について関係のあるすべての情報は、観測されたxに対する尤度関数に含まれている。さらに2つの尤度関数がθの関数として互いに比例の関係にあるなら2つの関数はθについて同じ情報を含んでいる
    • フィッシャーに由来
    • これを満たすか満たさないかで、陣営が分かれる
      • 満たす
      • 満たさない
        • 頻度主義
    • いまだ論争の火種
    • もしLPがパラメータ推定の根本原理であるなら、
      • P(x' | θi) = kP(x | θi) が成り立つときEv(E, x) = Ev(E', x')となるはず。Evは証拠、Eは実験
      • これが成り立たなければそのような結果をもたらす推論は不適切
      • バーンバウムが証明
        • 証明の是非は未決着
    • さらに遡ると
      • 十分性の原理
      • 条件付けの原理
      • これら2つの原理は統計学者なら誰でも受け入れられるはずの原理
      • これら2つの原理と尤度原理の等価性を「証明」(バーンバウム)、2つの原理を受け入れるなら尤度原理も受け入れないといけない
    • バーンバウムの問題意識
      • 統計学的に導ける証拠」が「実験における証拠」になっている
        • (E, x) と Ev(E, x) は区別される
          • (E, x)
            • パラメータ空間Ωについての記述
            • Eの可能な結果xのサンプル空間についての記述
          • Ev(E, x):実験的証拠
            • どう評価するか?
      • 解明のポイント
        • 2つの統計的証拠(E, x) と (E', y) が関係するあらゆる点で等しいと言える条件?
          • 統計的証拠(E, x) と (E', y) が関係するあらゆる重要な点で等しい時、 Ev(E, x) = Ev(E', y)
  • 尤度原理からの重要な帰結
    • サンプルスペース(可能だが実際には得られなかった確率変数の値)の無関係性
      • 重要な争点
      • 反事実の重要性を否定

ベイズ主義

  • デフィネッティ
    • 確率、賭け、信念の関係
      • 「合理的な賭けが成立するための条件」「確率の3つの規則」同じ条件
    • D.ギリース「確率の哲学理論」
  • サヴェッジ
    • Inductive inference と Inductive behavior
    • 後者がより重要
      • Inference
        • 意見を変えること
      • Behavior
        • 分布と最終的な行為の経済的事実を用いて最も期待効用の高いものを選ぶ
  • 主観的 vs 客観的
    • 信念変化の合理性
    • 基礎付け主義
    • 客観的事前確率
  • 哲学的ベイズ主義
  • 実用への転換
    • Lindley & Smithによる階層ベイズ
    • Gelfand & SmithによるMCMCの導入

尤度主義とは

  • グレムリン仮説
    • 哲学ではしばしば「説明可能性」と「確率」が結びつく
  • 尤度主義の限界
    • 区間を持った仮説(複合仮説)は考察の対象とならない
    • 尤度の平均を求めることは可能だが、事前分布が必要となる
    • ソーバーの立場
      • 頻度主義の統計テストは用いたくない
      • ベイズ主義はできるだけ控えたい
      • 結果
  • 尤度もモデルに依存するため恣意性が混じるが、なぜ尤度主義者はベイズ主義を批判できるのか?
    • ミスリーディング確率というものを用いることで客観性を保証できる、というのがソーバーの立場
  • ソーバーのAICの法則、どれぐらい使われているのか?
    • あまり

対頻度主義

  • フィッシャーの有意性検定
    • 単独の仮説に対する検定
    • 帰無仮説は正しいと証明されることはない
    • 確率論的モーダストレンスへの批判
  • ネイマン - ピアソン
    • 行為選択の規則

GAMをもう少し理解したい④

前回の続きです。過去記事はこちらから。

ushi-goroshi.hatenablog.com

ushi-goroshi.hatenablog.com

ushi-goroshi.hatenablog.com

GAMの実装

bakfit

さて bakfit ですが、このサブルーチンはバックフィッティングというアルゴリズムを定義しています。このアルゴリズムは、簡単に言えば「関心のある平滑化対象変数を順番に残差に対してフィッティングしていく」方法です。Wikipediaバックフィッティングのページからアルゴリズムの部分について抜粋しましょう:

f:id:ushi-goroshi:20191125142457p:plain

ポイントは4行目の (a) のところで、


{\displaystyle \mathrm{Smooth} \left[\left\{ y_{i} - \hat{\alpha} - \sum_{k \neq j}\hat{f_{k}}(x_{ik})\right\}_{1}^{N}\right]}

とあるように、目的変数から定数部(線形部分)および関心のある平滑化対象変数"以外"(非線形部分)を引いています。すなわち「平滑化対象変数で説明される部分」と「確率的な誤差」で構成される部分のみを残し、平滑化を行なっていることがわかります。これを非線形としたい変数全てについて収束するまで順次繰り返す、というのがバックフィッティングの流れです。

ちなみにこのバックフィッティング、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.

で、これは以下のように解くことができます:


 {\displaystyle x_{i}^{(k+1)}={\frac {1}{a_{ii}}}\left(b_{i}-\sum _{j=1}^{i-1}a_{ij}x_{j}^{(k+1)}-\sum _{j=i+1}^{n}a_{ij}x_{j}^{(k)}\right),\quad i=1,2,\dots ,n.}

ここで括弧の中に着目すると、$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 ではざっくりと

  1. 前処理
  2. バックフィッティング
  3. 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

3つ目の工程で eta の再定義とあるのですが、これまでの処理(〜 s.wam)において eta は「 y をリンク関数で変換したもの」として扱われてきました。一方 backf1 の中では 「平滑化対象変数によって構成される部分」(非線形/平滑化部分)として扱われており、混乱しやすいので注意が必要です。つまり backf1 の中では y を「線形部分」と「非線形(平滑化)部分」に分割し、前者を etal 、後者を eta として計算が行われます。

では一つずつ見ていきましょう。

1. 前処理

まずは前処理です。ここでは

  • データに重みをつける
  • QR分解する
  • eta非線形部分)を定義する

という処理が行われます。

### 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 の integer1~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分解によって最小二乗解を求める、という流れになっているようです。

ushi-goroshi.hatenablog.com

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 に格納。ループ1回目の時点では etal は Nulldo 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 という関数のコアとなる部分では:

  1. バックフィッティングによる非線形項のフィッティング( bakfit および backf1
  2. バックフィッティングに渡すためにデータを順位に置き換える( s.wam
  3. {general, s, lo}.wam のいずれのエンジンに渡すかを特定し、 bf.call を作成( gam.fit
  4. gam.fit に渡すための諸々の設定( gam

といった処理が行われているようでした。またバックフィッティングというアルゴリズムは、端的には「関心のある変数による影響のみを残してフィッティングすることを繰り返す」というものであることがわかりました。

一方で、今回の一連の記事では平滑化の仕組みを追いかけることがまだ出来ていません。また平滑化の関数にしても s.wam しか見ておらず、 lo.wamgeneral.wam の挙動を確認していません。その辺りもいずれ調べたいですね。

樋口先生退任記念シンポジウム参加記録

11/5に開催された元統計数理研究所の所長である樋口先生の退任記念シンポジウムに参加してきましたので、そのメモを共有しておきます。ご参考まで。やっぱり長らく研究をされてきた方のお話というのは面白いですね。

ご挨拶

統数研 椿広計 氏

バイオインフォマティクスベイズモデリング

東京大学 井本清哉 氏

  • 樋口先生との共同研究
    • ベイジアンネットワークによる遺伝子ネットワークの推定
      • 遺伝子間の因果関係(制御関係)
      • グラフ構造の推定
      • ヒトゲノム、21000個の遺伝子
        • それらが複雑に関係し合う
        • たかだか有限のデータでは推定できない
      • ネットワーク構造の事後確率
      • 生物学的知識、他のデータで事前確率を構成
        • 結合配列
        • タンパク質の相互作用
        • 進化情報
        • データベース
        • 遺伝子ノックダウン
      • バイオインフォマティクスベイズモデリングの格好の応用対象であり、問題の宝庫
      • マイクロバイオームを含んで「私」
      • 健康や病気をデータ化できる時代が来た
    • データ同化による細胞内分子シミュレーションモデルの構築
    • インフルエンザ感染拡大エージェントベースシミュレーション
      • インフルエンザの拡大を抑えるためには会社員のワクチン接種が有効
  • 予測だけを目的としたら、AI(機械学習技術)は人を超えたかも
    • Kaggle、AutoMLがコンペで2位になった
  • ゲノム医療、理由が説明できないと医療に入っていくことは難しい現状
    • 医療分野で AI は Augumented Intelligence
    • 人の能力を向上させるもの
  • 医科学研究所、2900万件の医学論文を学習
  • 思考のリミッターを外しなさい by 樋口先生

消費者行動理解のためのベイジアンモデリング

筑波大学 佐藤忠彦 氏

  • マーケティングに関わる領域も最近は樋口先生の表芸になってきている
  • 能書き
    • マーケティング、現実の市場で発生する現象
    • 問題を解決するために以下の問いに答える必要がある
      • マーケティングを取り囲む環境は?
      • 市場はどこにある?
      • 敵は誰?
      • 顧客はどこに?
      • 顧客の嗜好は?
      • 消費者の価格反応は?
      • 広告、販促への反応は?
      • ものやサービスを選ぶ基準は?
    • 上記は経験や感覚で評価できるか?
      • 科学的アプローチが必要
    • 消費者を理解する = 仮説を発見する
    • 消費者を理解する ≠ 真の姿を解明する
    • 消費者理解のプロセス
      • 複数のモデルを立て
      • 全てのモデルを推定・比較し
      • この時点でもっとも妥当と考えられる構造の推定
      • 情報の変換
      • 結果の活用
    • 現在持ちうる知識で記述能力・汎化能力が高いモデルを選ぶ
    • 消費者理解のキーワード
      • 消費者異質性
      • 時間的異質性
      • 潜在変数
    • パラメータがいっぱいで普通のモデルでは対応できない
    • 事後分布が関心の対象
    • 樋口流ベイジアンモデリング
      • 理論のみならずドメイン知識に基づく
      • こんな構造はないか、あんな構造はないか、という構成的なモデリング
  • マーケティングの事例①:佐藤、樋口(2009)
    • スーパーで牛乳を購買する確率が、潜在変数である家庭内在庫量/消費量の推移に基づき変化する構造
    • 一般状態空間
    • 店に来て牛乳を買う状況
    • n時点での在庫量 = 前日開始時点での在庫量 + 購買量 - 消費量
    • 消費量 = 残っている量 * 消費率
    • 店に来る効用、牛乳を買う効用
    • 動的に変化する
    • 在庫量が変化すると購入確率がどのように変化するか
  • マーケティングの事例②:宮津、佐藤(2015)
    • スーパーにおける購買個数が、心理的財布の状況により変化する構造
    • 給料日直前、給料日
    • 25、10、17日を起点とした3つの累積購買金額の線形結合で柔軟に構成する
    • 閾値ポアソン回帰
      • 心理的財布が閾値に到達しているか否かで行動が変わる
  • まとめ
    • Supply Side Thinking ではなく Demand Side Thinking
    • ぐにゃぐにゃモデルのすすめ
    • ベイズモデリングが有効

固体地球科学とデータ同化

東京大学 長尾大道 氏

  • 北川源四郎 時系列解析プログラミング
  • データ同化
    • 数値シミュレーションと観測データの融合
    • 大規模な状態空間モデルと大規模データ
  • 四次元変分法
  • データ同化入門

樋口先生とブリジストンと私

ブリジストン 花塚泰史 氏

  • 樋口先生と10年以上の関係がある
  • タイヤ事業が8割、多角化事業が2割
  • 売り切り型からソリューションプロバイダーとして
  • 状態空間モデルとカーネルマシンを併用した・・・
  • ブリジストンでの事例
    • タイヤ、ほぼ唯一の電子化されていない部品
    • CAISによる安全・安心なモビリティ社会への貢献
      • インテリジェント化したタイヤ
      • 空気圧、温度を拾うのに加え、タイヤ内のセンサーから
        • 磨耗状態推定
        • 路面状態推定 → 自動車の制御へフィードバック
          • 加速度
          • 路面状態に応じて波形が特徴的に変化
            • ウェット、凍結
            • 周波数領域において特徴的な波形
          • Windowをシフトしながら
          • カーネルSVM
            • (速度によって)系列長が異なるものを比較
            • グローバルアライメントカーネル
              • 計算量がでかい
              • リアルタイムでやりたい
            • 窓同士の類似度
            • ドライとウェットを分類する判別問題
          • GAカーネルの高速化
            • 類似度を計算する窓を制限する
            • タイヤの波形は特徴的な2つのピークが出るので・・・
          • 元のGA、93.9%で判別できるがタイヤ100回転分の時間が必要
          • 工夫により96.9%、タイヤ一回転分ぐらいで判断できるようになった
            • 7種類に判別
          • 冬季道路管理の最適化への適用
            • 路面のメンテナンス
              • 巡回車に乗っている人の目で判断
            • 北海道で実用化されている
              • ISCOS
    • 次世代のデータサイエンティスト育成
      • 社内研修
  • 将来的な発展
    • センシング技術
    • 解析、予測
    • デジタルサービス

ベイズモデリングと歩んだ30年

中央大学 樋口知之 氏

  • 平成元年に入所、平成31年に退任
    • 200MBのデータ解析に大型計算機が必要だった時代
  • 人間の未来予測は当たらない
  • モデルとモデリング
    • 何が違うのか
      • 常にモデルを改良する姿勢
      • ものの見方、捉え方を柔軟に変化
      • 思考のプロセス自体を科学する
    • モデルを比較するための羅針盤としての情報量基準
  • 赤池先生との出会い
    • 多変量自己回帰を応用したセメント工場、火力制御
    • 船のオートパイロット(北川先生)
  • 当時
    • 支配方程式がない、計算コストが大
  • 生成モデルと識別モデル
    • 生成モデルはデータが生成される過程
    • 識別モデルはデータやパラメータを与えたときにどのように識別されるか
    • 後者が簡単(機械学習コミュニティで盛ん)
    • 統計の人たちはデータがどんな機序で生まれたのかに強い関心がある
  • 頻度主義 vs ベイズ主義
  • ベイズの定理がなぜ役立つか?
    • コンピュータの性能向上
      • 周辺尤度の計算
    • 高精度センサーのコモディティ化
      • 尤度関数
    • ストレージの廉価化
      • 細かい情報を保存
  • 1980年代の統数研
  • 統数研
    • 統数研の助手は世界でもっとも恵まれている
    • 人材採用方針
      • 数理統計、確率論
      • 物理、計算科学、計算機科学
      • データの質、量の観点で大きく変化のある分野から
    • 赤池先生
  • アカデミアから産業界への人材流出は欧米でも同じ
    • クロスアポイント
  • メゾスコピックモデリング
    • 数理表現された法則
    • 第一原理法則
  • 演繹と帰納
  • 逐次データ同化技術
    • KFとPFの間にある
    • 状態ベクトルの次元数、アンサンブルメンバー数
  • 統計的推論法のパラダイムシフトを促す外的要因
    • ビッグデータ
    • データ駆動サービス社会
      • 因果から相関へ
    • 目的特化型計算機
  • 分析者の技術から計算パワーへ
  • 力学的ダウンスケーリング
  • エミュレーション、仮想計測
  • やっていること
    • GANを使ってエミュレータを作る
    • 深層学習と状態空間モデルの融合
    • 高等教育
  • 機械学習と統計
    • Data driven vs Science driven
    • Predictive model vs Descriptive model
    • Correlation vs Causality
    • 認識科学から設計科学
    • 対象理解から帰納の最適化
    • 「真のモデルなど存在しない」という見方は、当面の問題のイメージ能力の欠如

AI社会での数理モデリングと現場主義の意義と価値

東京大学 合原一幸 氏 東北大学 早瀬敏幸 氏 大阪大学 鷲尾 隆 氏 ブリジストン 花塚泰史 氏 中央大学 樋口知之 氏

  • 分野横断的なアプローチで豊富な共同研究実績の持ち主
  • 合原氏
    • 数理工学
      • 統計とは兄弟みたいな
      • 暮らしを変える驚きの数理工学
      • ニューラルコンピュータ
      • アルファ碁解体新書
    • AIと人間の共同作業、協働
      • Pair Go
      • ファッション
        • GANでデザイナーっぽいデザイン
  • 早瀬氏
    • 流体力学、制御工学、生体工学
    • データ同化
  • 鷲尾氏
    • 演繹的なモデリングの自動化
      • 要素還元主義・普遍主義
      • 外挿に利用可能
      • 深い理解が可能
    • 帰納的なモデリング
      • 高精度な近似はできるが外挿できない
      • XAIの観点から、浅い理解しかできない
    • 演繹的モデルと帰納的モデルのより広範な融合
  • 三者は生成モデルにこだわりがある
  • 3つのお題(6つの事前質問から)
    • 現場主義
      • 1.現場主義とKKDの違いは? 現場主義に重要性について
      • 2.学生や研究者に現場主義をどう教育していくべきか?
        • 樋口
          • 現場に出向く研究者魂
          • 現場の人と厳しく議論できるコミュニケーション力、分野跳躍力
          • 本質的な課題を数理の問題にしたてるセンス
          • しつこさ
        • 鷲尾
          • 向こうの分野の研究者に統計や機械学習を教えた方が良いという経験
          • 現場を知っている人の方がよく使える
        • 合原
          • 中小企業において、社長が現場を見に行く
            • そうすると階層構造にショートカットが生まれる
        • 早瀬
          • 研究時間が減っていて、研究に割ける時間がない
        • 花塚
          • ブリジストンでの取り組み
            • 工場のオペレーター、タイヤ館の店員
            • 何に困っているか、ソリューションが役にたつかを話を聞きに行く
            • 障害
              • 相手が自分のナレッジを理解していない
              • KKDが正しいかの見極め
            • ドメインの知識を作り手側が持っていないといけない
      • コミュニケーションが取りづらい相手は?
        • 無茶振りする(合原)
          • 専念しても一年かかるようなもの
          • データ化されていないものを大量に持ってくる
          • 割ける時間の中でうまく成果が出るテーマ選び
          • 他の人に無茶振りしてみる
            • その人だったら無茶ではないかもしれない
            • 人脈づくり、人的ネットワーク
              • 有名な人を知っているだけではダメ
              • あの人はいま時間がある
        • 分析・テーマの振り返りが繰り返されるので、じっくり取り組む必要がある(鷲尾)
    • 人とAI
      • 3.複雑な現象を数理モデルで表現する意義は?
      • 4.今後の研究スタイルはどう変わっていくか?
        • 合原
          • Data drivenは実用上は問題ない
          • そこから法則をスッキリ抜き出すことができれば面白い
          • シンプルな非線形性で抜き出したい
          • ニューラルネットは完全にTransparent
            • 言葉にできないだけ
        • 鷲尾
          • 機械学習の研究者が精度を競っているので精一杯
          • 最近XAIが成長してきた
          • 中身を理解しやすいモデルはどうやって作るか
          • 現在の機械学習は、現実との引き合わせができていない
        • 花塚
          • 製造業においてはケースバイケース
        • 早瀬
          • 順モデルなのか逆モデルなのか
          • どのような過程で生まれたのか説明される必要がある
          • 現実のシステムがもつ因果
    • 教育システム
      • 5.分野横断的な学術体系をどのようにして教育システムに組み込むか
      • 6.10年後に国立大学はどのような姿になっているか
        • 東北大学では全員にデータサイエンス系の講義を受けさせることになる(早瀬)
        • オーストラリアの某大学は数学科を三等分している(合原)
        • 若い人を早めに独立させる(合原)
          • 教授、准教授、助教みたいな階層構造をやめる
        • 数学が苦手な人たちにいかに脱落させないか(鷲尾)
          • OJTとして授業より前に使ってもらう
        • 古典的製造業は縦割りがピシッとしている(花塚)
          • 物売りからこと売り、ソリューションやるなら横串が必要
          • サプライチェーンの初めから終わりまでを担う
            • それがデータサイエンス組織
          • コミュニケーションを密にして話をしにいく
          • 現場ごとにデータサイエンスを理解している人が必要
          • CoEとしての機能
          • 育成
            • CoE
              • 高度な手法を理解できる
            • 現場
              • 話を理解できる人を多く育てる
              • 裾野を広げる
          • 知っている人をいかにつなげていくか
      • 統数研に期待すること、文科省に頑張ってほしいこと
        • 予算(合原)
        • 日本のトップの研究所として世界をリードする研究(早瀬)
        • 地道に学問を究めるところが魅力(鷲尾)
          • 外にアピールして外部資金を獲得
          • 基礎研究を地道にコツコツと続けられる
        • リカレント教育(花塚)
          • DSを外から取ってくるのは至難
          • 社内だけでの教育も困難
            • アカデミアとの連携