glmnetをもう少し理解したい②
前回の記事では glmnet
の中身を確認し、引数の family
によって呼び出す関数を変えていることがわかりました。
今回はそのなかでも gaussian
が指定された場合の関数である elnet
を見ていきましょう。
なお前回の記事はこちらです。
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
はさらに先の処理で elnetu
と elnetn
という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
が指定されているか否かで spelnet
と elnet
のどちらが呼ばれるかが決まりますが、引数の違いとしては spelnet
において x
が as.double
とされており、 ix
と jx
(いずれも疎行列において非ゼロの要素の座標を特定するための数値)が追加されています。
# 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")
として呼ばれていました。これまで glm
や GAM
で見てきたときと同じように、 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 } }
として取り出せます。
関数をみてみると、 errlist
は switch(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
の呼びだしを確認すると、elnetu
と elnetn
のいずれを呼ぶかは ka
で決まっています。
先ほど少し触れた通り、 ka
は ka = 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
を見てみましょう。