请问SCAD惩罚函数是怎么把参数估计为0的？我这里有Python代码

是 否 +2 论坛币

k人 参与回答

经管之家送您一份

应届毕业生专属福利!

求职就业群

赵安豆老师微信：zhaoandou666

经管之家联合CDA

送您一个全额奖学金名额~ !

立即领取

感谢您参与论坛问题回答

经管之家送您两个论坛币！

+2 论坛币

我用Python按照：https://ask.hellobi.com/blog/datasciencemeditation/7540，http://www.doc88.com/p-2089722981617.html。这两个链接实现门槛回归（没办法，没看到Python有什么库提供门槛回归方法），由于第2个链接里的论文根本没提到SCAD惩罚函数的样子，我就去搜索到了第一个链接。总而言之，现在我的确跑出来结果了，但是门槛值 c 的估计值竟然比实际存在的门槛变量的取值还大好多倍，这肯定不正确。。。
所以我怀疑是不是我对SCAD惩罚函数的理解有错误，总之我今天又谷歌，搜索到了http://vlambda.com/wz_wRtUGSimSJ.html。这篇文章，里面的SCAD惩罚函数又和第一个链接看到的不一样，我已经迷糊了

我的意思是，现在我搞不清楚SCAD是不是函数表达式弄错了，但不论是否弄错，我发现我不能理解SCAD如何把参数估计为0。参数估计为0就会将某段的回归全部系数估计为0，从而能直接告知这个分段不需要，这是第二个链接论文告诉我的。。。但反正，现在门槛 c 值就不对。。。门槛自变量的数据最大也才1W出头，而估计出的C竟然数倍于它们
以下是代码，其中pay_data_x为189*9的矩阵自变量数据，pay_data_y为因变量。M_c是分段，这里分了2段，先分别用pay_data_x的第2列数据的标准差和均值做分段门槛值的初始值（即第2列数据为门槛变量）
另外就是，是否有可能我对门槛值c的求导有错误？？

1. # 定义损失函数：光滑SCAD正则，参考《门槛回归模型中门槛值和回归参数的估计》
   
2. def loss_SCAD(beta_matrix, M_c):
   
3.     first_M_MSE = np.dot(pay_data_x, beta_matrix[:lt]) - pay_data_y.values.reshape((-1, 1))
   
4.     z = pay_data_x.iloc[:, 1].values.reshape((-1, 1))
   
5.     other_M_MSE = 0
   
6.     scad = 0
   
7. 
   
8.     for idx, c in enumerate(M_c):
   
9.         nor_distri_x = (z - c) / h
   
10.         nor_distri_CDF = st.norm.cdf(nor_distri_x)  # 光滑示性函数
    
11.         other_M_MSE += np.dot(pay_data_x, beta_matrix[lt * (idx + 1):lt * (idx + 2)]) * nor_distri_CDF
    
12.         for beta_param in beta_matrix[lt * (idx + 1):lt * (idx + 2)]:
    
13.             if abs(beta_param) < lam:
    
14.                 scad += lam ** 2 - (abs(beta_param) - lam) ** 2
    
15.             else:
    
16.                 scad += lam ** 2
    
17. 
    
18.     MSE = 0.5 * np.average((first_M_MSE + other_M_MSE) ** 2) + scad
    
19. 
    
20.     return MSE
    
21. 
    
22. 
    
23. # 更新梯度
    
24. def loss_derivative(beta_matrix, M_c):
    
25.     deriv_list = []
    
26.     other_M_MSE = 0
    
27.     nor_distri_CDF = 0
    
28. 
    
29.     first_M_MSE = np.dot(pay_data_x, beta_matrix[:lt]) - pay_data_y.values.reshape((-1, 1))
    
30.     z = pay_data_x.iloc[:, 1].values.reshape((-1, 1))
    
31. 
    
32.     for idx, c in enumerate(M_c):
    
33.         nor_distri_x = (z - c) / h
    
34.         nor_distri_CDF = st.norm.cdf(nor_distri_x)  # 光滑示性函数
    
35.         other_M_MSE += np.dot(pay_data_x, beta_matrix[lt * (idx + 1):lt * (idx + 2)]) * nor_distri_CDF
    
36. 
    
37.     square_part = first_M_MSE + other_M_MSE
    
38. 
    
39.     for idx, beta_param in enumerate(beta_matrix[:lt]):  # 更新梯度：第一段
    
40.         deriv_list.append(np.average(square_part * pay_data_x.iloc[:, idx].values.reshape((-1, 1))))
    
41. 
    
42.     for idx, c in enumerate(M_c):
    
43.         for i, beta_param in enumerate(beta_matrix[lt * (idx + 1):lt * (idx + 2)]):  # 更新梯度：分段
    
44.             if abs(beta_param) > lam:  # 更新梯度：SCAD函数
    
45.                 if a * lam > beta_param:
    
46.                     scad = (lam + (a * lam - beta_param) / ((a - 1) * lam))
    
47.                 else:
    
48.                     scad = 0
    
49.             else:
    
50.                 scad = lam
    
51. 
    
52.             deriv_list.append(np.average(square_part * pay_data_x.iloc[:, i].values.reshape((-1, 1)) * nor_distri_CDF) + scad)
    
53. 
    
54.     for idx, c in enumerate(M_c):  # 更新梯度：门槛参数
    
55.         nor_distri_PDF = st.norm.pdf((z - c) / h)
    
56.         deriv_list.append(np.average(square_part * np.dot(pay_data_x, beta_matrix[lt * (idx + 1):lt * (idx + 2)]) * nor_distri_PDF / -h))
    
57. 
    
58.     return np.array(deriv_list).reshape((-1, 1))
    
59. 
    
60. 
    
61. # DFP迭代
    
62. def DFP(M_c, t=0.5, epsilon=1e-5, max_iter=1000, sigma=0.4):
    
63.     c = np.array(M_c).reshape((-1, 1))
    
64.     tmp_beta = np.ones((pay_data_x.shape[1] * (len(M_c) + 1), 1))
    
65.     H = np.eye(pay_data_x.shape[1] * (len(M_c) + 1) + len(M_c))  # 初始化伪黑塞矩阵的逆
    
66.     beta_c_last = np.concatenate((tmp_beta, c))
    
67.     iter = 1
    
68.     tmp_beta_idx = tmp_beta.shape[0]
    
69. 
    
70.     while iter < max_iter:
    
71.         iter += 1
    
72. 
    
73.         beta_last = beta_c_last[:tmp_beta_idx]
    
74.         c_last = beta_c_last[tmp_beta_idx:]
    
75.         deriv_last = loss_derivative(beta_last, c_last)
    
76. 
    
77.         if np.linalg.norm(deriv_last) <= epsilon or iter >= max_iter:
    
78.             return beta_last.reshape((-1, len(M_c) + 1)), c_last, iter, loss_SCAD(beta_last, c_last)
    
79. 
    
80.         dk = -np.dot(H, deriv_last)
    
81. 
    
82.         m, mk = 0, 0
    
83.         loss_Armijo_1 = loss_SCAD(beta_last, c_last)
    
84.         while m < 20:  # 搜索最佳学习率
    
85.             if loss_SCAD(beta_last + t ** m * dk[:tmp_beta_idx], c_last + t ** m * dk[tmp_beta_idx:]) \
    
86.                     < loss_Armijo_1 + sigma * t ** m * np.dot(deriv_last.T, dk):
    
87.                 mk = m
    
88.                 break
    
89.             m += 1
    
90. 
    
91.         beta_c_new = beta_c_last + t ** mk * dk
    
92.         beta_new = beta_c_new[:tmp_beta_idx]
    
93.         c_new = beta_c_new[tmp_beta_idx:]
    
94. 
    
95.         deriv_new = loss_derivative(beta_new, c_new)
    
96.         s = beta_c_new - np.concatenate((beta_last, c_last))  # 参数一阶差分
    
97.         y = deriv_new - deriv_last  # 梯度一阶差分
    
98.         h_y = np.dot(H, y)
    
99.         H += np.dot(s, s.T) / np.dot(s.T, y) - np.dot(h_y, h_y.T) / np.dot(np.dot(y.T, H), y)  # 更新伪黑塞矩阵的逆
    
100. 
     
101.         beta_c_last = beta_c_new
     
102. 
     
103. 
     
104. s = DFP(M_c=(pay_data_y.std(), pay_data_y.mean()))
     
105. print(s[0], s[1], s[2], s[3], sep='\n')

复制代码

分3段，故有9*3=27个参数，可以看到第2和第3段参数估计一模一样

[[ 2.21651758e+02  1.24799988e+04  1.24799988e+04]

[-2.70472227e-02  1.26663922e+08  1.26663922e+08]

[ 2.84548849e-02  3.16436052e+07  3.16436052e+07]

[ 3.82476892e-02  3.05133523e+07  3.05133523e+07]

[ 1.13002044e-02  2.92317240e+07  2.92317240e+07]

[ 8.44805488e-02  3.13574382e+07  3.13574382e+07]

[-1.01630661e-01  2.93887077e+07  2.93887077e+07]

[ 3.42105205e-02  3.08086230e+07  3.08086230e+07]

[ 7.84099476e+00  2.79222208e+07  2.79222208e+07]]

门槛值c：[[ 26632.39914284]， [113303.41470849]]

扫码加我 拉你入群

请注明：姓名-公司-职位

以便审核进群资格，未注明则拒绝

分享0 收藏0 回帖

我最近也在看这个方法。不想调程序，想自己编一个。你现在搞明白了吗？

阿卜拉克萨斯 发表于 2019-1-3 10:49
我用Python按照：https://ask.hellobi.com/blog/datasciencemeditation/7540，http://www.doc88.com/p-2089 ...

楼主看范老师的原文吧  原文里面写的挺清楚的   2001年的JASA

画一下惩罚函数导数的图像，可以看到在小于lam那一段是赋予大的权重，这样就能达到收缩系数了