变量选择之SCAD算法
本文提出了一种用于同时达到选择变量和预测模型系数的目的的方法——SCAD。这种方法的罚函数是对称且非凹的,并且可处理奇异阵以产生稀疏解。此外,本文提出了一种算法用于优化对应的带惩罚项的似然函数。这种方法具有广泛的适用性,可以应用于广义线性模型,强健的回归模型。借助于波和样条,还可用于非参数模型。更进一步地,本文证明该方法具有Oracle性质。模拟的结果显示该方法相比主流的变量选择模型具有优势。并且,模型的预测误差公式显示,该方法实用性较强。
SCAD的理论理解
在总结了现有模型的一些缺点之后,本文提出构造罚函数的一些目标:
罚函数是奇异的(singular)
连续地压缩系数
对较大的系数产生无偏的估计
SCAD模型的Oracle性质,使得它的预测效果跟真实模型别无二致。
并且,这种方法可以应用于高维非参数建模。
SCAD的目标函数如下:
SCAD的罚函数与$\theta$的(近似)关系如下图所示。
可见,罚函数可以用二阶泰勒展开逼近。
Hard Penality,lasso,SCAD的系数压缩情况VS系数真实值的情况如下图所示。
可以看到,lasso压缩系数是始终有偏的,Hard penality是无偏的,但压缩系数不连续。而SCAD既能连续的压缩系数,也能在较大的系数取得渐近无偏的估计。
这使得SCAD具有Oracle性质。
SCAD的缺点
模型形式过于复杂
迭代算法运行速度较慢
在low noise level的情况下表现较优,但在high noise level的情况下表现较差。
SCAD的实现
SCAD迭代公式
SCAD的目标函数如下:
时,罚函数可以用二阶泰勒展开逼近。
从而,有如下迭代公式:
根据以上公式,代入迭代步骤,即可实现算法。
SCAD的R实现
##------数据模拟--------
library(MASS)
##mvrnorm()
##定义一个产生多元正态分布的随机向量协方差矩阵
Simu_Multi_Norm<-function(x_len, sd = 1, pho = 0.5){
#初始化协方差矩阵
V <- matrix(data = NA, nrow = x_len, ncol = x_len)
#mean及sd分别为随机向量x的均值和方差
#对协方差矩阵进行赋值pho(i,j) = pho^|i-j|
for(i in 1:x_len){ ##遍历每一行
for(j in 1:x_len){ ##遍历每一列
V[i,j] <- pho^abs(i-j)
}
}
V<-(sd^2) * V
return(V)
}
##产生模拟数值自变量X
set.seed(123)
X<-mvrnorm(n = 200, mu = rep(0,10), Simu_Multi_Norm(x_len = 10,sd = 1, pho = 0.5))
##产生模拟数值:响应变量y
beta<-c(1,2,0,0,3,0,0,0,-2,0)
#alpha<-0
#prob<-exp(alpha + X %*% beta)/(1+exp(alpha + X %*% beta))
prob<-exp( X %*% beta)/(1+exp( X %*% beta))
y<-rbinom(n = 200, size = 1,p = prob)
##产生model matrix
mydata<-data.frame(X = X, y = y)
#X<-model.matrix(y~., data = mydata)
##包含截矩项的系数
#b_real<-c(alpha,beta)
b_real<-beta
########----定义惩罚项相关的函数-----------------
##定义惩罚项
####运行发现,若lambda设置为2,则系数全被压缩为0.
####本程序根据rcvreg用CV选出来的lambda设置一个较为合理的lambda。
p_lambda<-function(theta,lambda = 0.025){
p_lambda<-sapply(theta, function(x){
if(abs(x)< lambda){
return(lambda^2 - (abs(x) - lambda)^2)
}else{
return(lambda^2)
}
}
)
return(p_lambda)
}
##定义惩罚项导数
p_lambda_d<-function(theta,a = 3.7,lambda = 0.025){
if(abs(theta) > lambda){
if(a * lambda > theta){
return((a * lambda - theta)/(a - 1))
}else{
return(0)
}
}else{
return(lambda)
}
}
# ##当beta_j0不等于0,定义惩罚项导数近似
# p_lambda_d_apro<-function(beta_j0,beta_j,a = 3.7, lambda = 2){
# return(beta_j * p_lambda_d(beta = beta_j0,a = a, lambda = lambda)/abs(beta_j0))
# }
#
#
# ##当beta_j0 不等于0,指定近似惩罚项,使用泰勒展开逼近
# p_lambda_apro<-function(beta_j0,beta_j,a = 3.7, lambda = 2){
# if(abs(beta_j0)< 1e-16){
# return(0)
# }else{
# p_lambda<-p_lambda(theta = beta_j0, lambda = lambda) +
# 0.5 * (beta_j^2 - beta_j0^2) * p_lambda_d(theta = beta_j0, a = a, lambda = lambda)/abs(beta_j0)
# }
# }
#define the log-likelihood function
loglikelihood_SCAD<-function(X, y, b){
linear_comb<-as.vector(X %*% b)
ll<-sum(y*linear_comb) + sum(log(1/(1+exp(linear_comb)))) - nrow(X)*sum(p_lambda(theta = b))
return (ll)
}
##初始化系数
#b0<-rep(0,length(b_real))
#b0<- b_real+rnorm(length(b_real), mean = 0, sd = 0.1)
##将无惩罚时的优化结果作为初始值
b.best_GS<-b.best
b0<-b.best_GS
##b1用于记录更新系数
b1<-b0
##b.best用于存放历史最大似然值对应系数
b.best_SCAD<-b0
# the initial value of loglikelihood
ll.old<-loglikelihood_SCAD(X = X,y = y, b = b0)
# initialize the difference between the two steps of theta
diff<-1
#record the number of iterations
iter<-0
#set the threshold to stop iterations
epsi<-1e-10
#the maximum iterations
max_iter<-100000
#初始化一个列表用于存放每一次迭代的系数结果
b_history<-list(data.frame(b0))
#初始化列表用于存放似然值
ll_list<-list(ll.old)
#######-------SCAD迭代---------
while(diff > epsi & iter < max_iter){
for(j in 1:length(b_real)){
if(abs(b0[j]) < 1e-06){
next()
}else{
#线性部分
linear_comb<-as.vector(X %*% b0)
#分子
nominator<-sum(y*X[,j] - X[,j] * exp(linear_comb)/(1+exp(linear_comb))) +
nrow(X)*b0[j]*p_lambda_d(theta = b0[j])/abs(b0[j])
#分母,即二阶导部分
denominator<- -sum(X[,j]^2 * exp(linear_comb)/(1+exp(linear_comb))^2) +
nrow(X)*p_lambda_d(theta = b0[j])/abs(b0[j])
#2-(3) :更新b0[j]
b0[j]<-b0[j] - nominator/denominator
#2-(4)
if(abs(b0[j]) < 1e-06){
b0[j] <- 0
}
# #更新似然值
# ll.new<- loglikelihood_SCAD(X = X, y = y, b = b0)
#
#
#
# #若似然值有所增加,则将当前系数保存
# if(ll.new > ll.old){
# #更新系数
# b.best_SCAD[j]<-b0[j]
# }
#
# #求差异
# diff<- abs((ll.new - ll.old)/ll.old)
# ll.old <- ll.new
# iter<- iter+1
# b_history[[iter]]<-data.frame(b0)
# ll_list[[iter]]<-ll.old
# ##当达到停止条件时,跳出循环
# if(diff < epsi){
# break
# }
#
}
}
#更新似然值
ll.new<- loglikelihood_SCAD(X = X, y = y, b = b0)
#若似然值有所增加,则将当前系数保存
if(ll.new > ll.old){
#更新系数
b.best_SCAD<-b0
}
#求差异
diff<- abs((ll.new - ll.old)/ll.old)
ll.old <- ll.new
iter<- iter+1
b_history[[iter]]<-data.frame(b0)
ll_list[[iter]]<-ll.old
}
b_hist<-do.call(rbind,b_history)
#b_hist
ll_hist<-do.call(rbind,ll_list)
#ll_hist
#
iter
##
ll.best<-max(ll_hist)
ll.best
##
b.best_SCAD
##对比
cbind(coeff_glm,b.best,b.best_SCAD,b_real)
##----------ncvreg验证-----------
library(ncvreg)
my_ncvreg<-ncvreg(X,y,family = c("binomial"),penalty = c("SCAD"),lambda = 2)
my_ncvreg$beta
my_ncvreg<-ncvreg(X,y,family = c("binomial"),penalty = c("SCAD"))
summary(my_ncvreg)
my_ncvreg$beta
###用cv找最优的lambda
scad_cv<-cv.ncvreg(X,y,family = c("binomial"),penalty='SCAD')
scad_cv$lambda.min
mySCAD=ncvreg(X,y,family = c("binomial"),penalty='SCAD',lambda=scad_cv$lambda.min)
summary(mySCAD)
ncv_SCAD<-mySCAD$beta[-1]
##对比
myFinalResults<-cbind(无惩罚项回归=coeff_glm, GS迭代 = b.best,
GS_SCAD迭代 = b.best_SCAD, ncvreg = ncv_SCAD,真实值 = b_real)
save(myFinalResults,file = "myFinalResults.rda")
数据分析咨询请扫描二维码
若不方便扫码,搜微信号:CDAshujufenxi
在职场中,许多言辞并非表面意思那么简单,有时需要听懂背后的“潜台词”。尤其在数据分析的领域里,掌握常用术语就像掌握一门新 ...
2025-01-04在当今信息化社会,数据分析已成为各行各业的核心驱动力。它不仅仅是对数字进行整理与计算,而是在数据的海洋中探寻规律,从而指 ...
2025-01-03又到一年年终时,各位打工人也迎来了展示成果的关键时刻 —— 年终述职。一份出色的年终述职报告,不仅能全面呈现你的工作价值, ...
2025-01-03在竞争激烈的商业世界中,竞品分析对于企业的发展至关重要。今天,我们就来详细聊聊数据分析师写竞品分析的那些事儿。 一、明确 ...
2025-01-03在数据分析的江湖里,有两个阵营总是争论不休。一派信奉“大即是美”,认为数据越多越好;另一派坚守“小而精”,力挺质量胜于规 ...
2025-01-02数据分析是一个复杂且多维度的过程,从数据收集到分析结果应用,每一步都是对信息的提炼与升华。可视化分析结果,以图表的形式展 ...
2025-01-02在当今的数字化时代,数据分析师扮演着一个至关重要的角色。他们如同现代企业的“解密专家”,通过解析数据为企业提供决策支持。 ...
2025-01-02数据分析报告至关重要 一份高质量的数据分析报告不仅能够揭示数据背后的真相,还能为企业决策者提供有价值的洞察和建议。 年薪 ...
2024-12-31数据分析,听起来好像是技术大咖的专属技能,但其实是一项人人都能学会的职场硬核能力!今天,我们来聊聊数据分析的核心流程,拆 ...
2024-12-31提到数据分析,你脑海里可能会浮现出一群“数字控”抱着电脑,在海量数据里疯狂敲代码的画面。但事实是,数据分析并没有你想象的 ...
2024-12-31关于数据分析师是否会成为失业高危职业,近年来的讨论层出不穷。在这个快速变化的时代,技术进步让人既兴奋又不安。今天,我们从 ...
2024-12-30数据分析师在现代企业中扮演着关键角色,他们的工作内容不仅丰富多样,还对企业的决策和发展起着重要的作用。正如一个经验丰富的 ...
2024-12-29数据分析师的能力要求 在当今的数据主导时代,数据分析师的角色变得尤为重要。他们不仅需要具备深厚的技术背景,还需要拥有业务 ...
2024-12-29随着技术的飞速发展与行业的持续变革,不少人心中都存有疑问:到了 2025 年,数据分析师还有前途吗?给你分享一篇阿里P8大佬最近 ...
2024-12-29如何构建数据分析整体框架? 要让数据分析发挥其最大效能,建立一个清晰、完善的整体框架至关重要。今天,就让我们一同深入探讨 ...
2024-12-27AI来了,数分人也可以很省力,今天给大家介绍7个AI+数据分析工具,建议收藏。 01酷表 EXCEL 网址:https://chatexcel.com/ 这是 ...
2024-12-26一个好的数据分析模型不仅能使分析具备条理性和逻辑性,而且还更具备结构化和体系化,并保证分析结果的有效性和准确性。好的数据 ...
2024-12-26当下,AI 的发展堪称狂飙猛进。从 ChatGPT 横空出世到各种大语言模型(LLM)接连上线,似乎每个人的朋友圈都在讨论 AI 会不会“ ...
2024-12-26数据分析师这个职业已经成为了职场中的“香饽饽”,无论是互联网公司还是传统行业,都离不开数据支持。想成为一名优秀的数据分析 ...
2024-12-26在数据驱动决策成为商业常态的今天,数据分析师这一职业正迎来前所未有的机遇与挑战。很多希望转行或初入职场的人士不禁询问:数 ...
2024-12-25