glmnet 使用什么偏差来比较λλ?

机器算法验证 r 网络
2022-04-05 20:07:55

的最佳值的一个标准是检查偏差对范围的图，并在偏差最小时选择（或在最低限度）。λλλλ

但是，我很难理解用 显示的确切内容，glmnet因为plot.cv.glmnet显示的图与绘制偏差对的结果完全不同。λ

set.seed(4567)
N       <- 500
P       <- 100
coefs   <- NULL
for(p in 1:P){
    coefs[p]    <- (-1)^p*100*2^(-p)
}
inv.logit <- function(x) exp(x)/(1+exp(x))
X   <- matrix(rnorm(N*P), ncol=P, nrow=N)
Y   <- rbinom(N, size=1, p=inv.logit(cbind(1, X)%*%c(-4, coefs)))
plot(test   <- cv.glmnet(x=X, y=Y, family="binomial", nfolds=10, alpha=0.8))
plot(log(test$lambda), deviance(test$glmnet.fit))

在此处输入图像描述 在此处输入图像描述

似乎第二个图没有包含弹性净惩罚，并且垂直缩放也不正确。我的主张是基于较大的值的曲线形状类似于输出的形状。然而，当我尝试自己计算惩罚时，我的尝试同样显得非常不准确。λglmnet

penalized.dev.fn    <- function(lambda, alpha=0.2, data, cv.model.obj){
    dev <- deviance(cv.model.obj$glmnet.fit)[seq_along(cv.model.obj$lambda)[cv.model.obj$lambda==lambda]]
    beta <- coef(cv.model.obj, s=lambda)[rownames(coef(cv.model.obj))!="(Intercept)"]
    penalty <- lambda * ( (1-alpha)/2*(beta%*%beta) + alpha*sum(abs(beta)) )
    penalized.dev <- penalty+dev
    return(penalized.dev)
}

out <- sapply(test$lambda, alpha=0.2, cv.model.obj=test, FUN=penalized.dev.fn)
    plot(log(test$lambda), out)

我的问题是：如何手动计算默认plot.cv.glmnet图表中报告的偏差？它的公式是什么，我在尝试计算它时做错了什么？

2个回答

我只是想添加到输入中，但目前没有简明的答案，评论太长了。希望这能提供更多的见解。

似乎感兴趣的函数在解压后的 glmnet 库中，名为 cv.lognet.R ,' 作者使用的，似乎与 cv.glmnet 计算二项式偏差的方式相匹配。

虽然我在论文中的任何地方都没有看到它，从跟踪 glmnet 代码到 cv.lognet，我收集到的是它使用了这里描述的称为 capped binomial deviance 的东西。

−[Ylog10⁡(E)+(1−Y)log10⁡(1−E)]

predmat 是每个 lambda 的上限概率值 (E, 1-E) 输出的矩阵，与 y 和 y 的补码值进行比较，得到 lp。然后将它们放入 2*(ly-lp)偏差形式中，并在交叉验证的保留折叠上进行平均，以获得您在第一张图像中绘制的 cvm - 平均交叉验证误差 - 和 cv 范围。

我认为手动偏差函数（第二个图）的计算方式与内部偏差函数（第一个图）的计算方式不同。

    # from cv.lognet.R

    cvraw=switch(type.measure,
    "mse"=(y[,1]-(1-predmat))^2 +(y[,2]-predmat)^2,
    "mae"=abs(y[,1]-(1-predmat)) +abs(y[,2]-predmat),
    "deviance"= {
      predmat=pmin(pmax(predmat,prob_min),prob_max)
      lp=y[,1]*log(1-predmat)+y[,2]*log(predmat)
      ly=log(y)
      ly[y==0]=0
      ly=drop((y*ly)%*%c(1,1))
      2*(ly-lp)

   # cvm output
   cvm=apply(cvraw,2,weighted.mean,w=weights,na.rm=TRUE)

所以我访问了 CRAN 站点并下载了我认为是glmnet 包的来源。在 ./glmnet/R/plot.cv.glmnet.R 中，您似乎会找到您想要的源代码。它非常简短，因此我将在此处粘贴，但最好自己检查一下以确保它确实是正在运行的代码。

plot.cv.glmnet=function(x,sign.lambda=1,...){
  cvobj=x
  xlab="log(Lambda)"
  if(sign.lambda<0)xlab=paste("-",xlab,sep="")
  plot.args=list(x=sign.lambda*log(cvobj$lambda),y=cvobj$cvm,ylim=range(cvobj$cvup,cvobj$cvlo),xlab=xlab,ylab=cvobj$name,type="n")
      new.args=list(...)
      if(length(new.args))plot.args[names(new.args)]=new.args
    do.call("plot",plot.args)
    error.bars(sign.lambda*log(cvobj$lambda),cvobj$cvup,cvobj$cvlo,width=0.01,col="darkgrey")
  points(sign.lambda*log(cvobj$lambda),cvobj$cvm,pch=20,col="red")
axis(side=3,at=sign.lambda*log(cvobj$lambda),labels=paste(cvobj$nz),tick=FALSE,line=0)
abline(v=sign.lambda*log(cvobj$lambda.min),lty=3)
    abline(v=sign.lambda*log(cvobj$lambda.1se),lty=3)
  invisible()
}
其它你可能感兴趣的问题
上一篇时间序列数据中的顺序模式匹配 下一篇多组的多重比较