1 # include <RcppArmadillo.h>
2 # include <cmath>
3 // [[Rcpp::depends(RcppArmadillo)]]
4 // [[Rcpp::plugins(cpp11)]]
5
6 // [[Rcpp::export]]
sgn(const double x)7 int sgn(const double x) {
8 return (x > 0) - (x < 0);
9 }
10
11 // [[Rcpp::export]]
mad(const arma::vec & x)12 double mad(const arma::vec& x) {
13 return 1.482602 * arma::median(arma::abs(x - arma::median(x)));
14 }
15
16 // Asymmetric huber regression adjusted to quantile tau for initialization
17 // [[Rcpp::export]]
updateHuber(const arma::mat & Z,const arma::vec & res,const double tau,arma::vec & der,arma::vec & grad,const int n,const double rob,const double n1)18 void updateHuber(const arma::mat& Z, const arma::vec& res, const double tau, arma::vec& der, arma::vec& grad, const int n, const double rob, const double n1) {
19 for (int i = 0; i < n; i++) {
20 double cur = res(i);
21 if (cur > rob) {
22 der(i) = -2 * tau * rob;
23 } else if (cur > 0) {
24 der(i) = -2 * tau * cur;
25 } else if (cur > -rob) {
26 der(i) = 2 * (tau - 1) * cur;
27 } else {
28 der(i) = 2 * (1 - tau) * rob;
29 }
30 }
31 grad = n1 * Z.t() * der;
32 }
33
34 // [[Rcpp::export]]
huberReg(const arma::mat & Z,const arma::vec & Y,const double tau,arma::vec & der,arma::vec & gradOld,arma::vec & gradNew,const int n,const int p,const double n1,const double tol=0.0001,const double constTau=1.345,const int iteMax=5000)35 arma::vec huberReg(const arma::mat& Z, const arma::vec& Y, const double tau, arma::vec& der, arma::vec& gradOld, arma::vec& gradNew, const int n, const int p,
36 const double n1, const double tol = 0.0001, const double constTau = 1.345, const int iteMax = 5000) {
37 double rob = constTau * mad(Y);
38 updateHuber(Z, Y, tau, der, gradOld, n, rob, n1);
39 arma::vec beta = -gradOld, betaDiff = -gradOld;
40 arma::vec res = Y - Z * beta;
41 rob = constTau * mad(res);
42 updateHuber(Z, res, tau, der, gradNew, n, rob, n1);
43 arma::vec gradDiff = gradNew - gradOld;
44 int ite = 1;
45 while (arma::norm(gradNew, "inf") > tol && ite <= iteMax) {
46 double alpha = 1.0;
47 double cross = arma::as_scalar(betaDiff.t() * gradDiff);
48 if (cross > 0) {
49 double a1 = cross / arma::as_scalar(gradDiff.t() * gradDiff);
50 double a2 = arma::as_scalar(betaDiff.t() * betaDiff) / cross;
51 alpha = std::min(std::min(a1, a2), 100.0);
52 }
53 gradOld = gradNew;
54 betaDiff = -alpha * gradNew;
55 beta += betaDiff;
56 res -= Z * betaDiff;
57 rob = constTau * mad(res);
58 updateHuber(Z, res, tau, der, gradNew, n, rob, n1);
59 gradDiff = gradNew - gradOld;
60 ite++;
61 }
62 return beta;
63 }
64
65 // [[Rcpp::export]]
standardize(arma::mat X,const arma::rowvec & mx,const arma::vec & sx1,const int p)66 arma::mat standardize(arma::mat X, const arma::rowvec& mx, const arma::vec& sx1, const int p) {
67 for (int i = 0; i < p; i++) {
68 X.col(i) = (X.col(i) - mx(i)) * sx1(i);
69 }
70 return X;
71 }
72
73 // Different kernels for low-dimensional conquer
74 // [[Rcpp::export]]
updateGauss(const arma::mat & Z,const arma::vec & res,arma::vec & der,arma::vec & grad,const double tau,const double n1,const double h1)75 void updateGauss(const arma::mat& Z, const arma::vec& res, arma::vec& der, arma::vec& grad, const double tau, const double n1, const double h1) {
76 der = arma::normcdf(-res * h1) - tau;
77 grad = n1 * Z.t() * der;
78 }
79
80 // [[Rcpp::export]]
updateGaussWeight(const arma::mat & Z,const arma::vec & weight,const arma::vec & res,arma::vec & der,arma::vec & grad,const double tau,const double n1,const double h1)81 void updateGaussWeight(const arma::mat& Z, const arma::vec& weight, const arma::vec& res, arma::vec& der, arma::vec& grad, const double tau,
82 const double n1, const double h1) {
83 der = weight % (arma::normcdf(-res * h1) - tau);
84 grad = n1 * Z.t() * der;
85 }
86
87 // [[Rcpp::export]]
updateLogistic(const arma::mat & Z,const arma::vec & res,arma::vec & der,arma::vec & grad,const double tau,const double n1,const double h1)88 void updateLogistic(const arma::mat& Z, const arma::vec& res, arma::vec& der, arma::vec& grad, const double tau, const double n1, const double h1) {
89 der = 1.0 / (1.0 + arma::exp(res * h1)) - tau;
90 grad = n1 * Z.t() * der;
91 }
92
93 // [[Rcpp::export]]
updateUnif(const arma::mat & Z,const arma::vec & res,arma::vec & der,arma::vec & grad,const int n,const double tau,const double h,const double n1,const double h1)94 void updateUnif(const arma::mat& Z, const arma::vec& res, arma::vec& der, arma::vec& grad, const int n, const double tau, const double h,
95 const double n1, const double h1) {
96 for (int i = 0; i < n; i++) {
97 double cur = res(i);
98 if (cur <= -h) {
99 der(i) = 1 - tau;
100 } else if (cur < h) {
101 der(i) = 0.5 - tau - 0.5 * h1 * cur;
102 } else {
103 der(i) = -tau;
104 }
105 }
106 grad = n1 * Z.t() * der;
107 }
108
109 // [[Rcpp::export]]
updatePara(const arma::mat & Z,const arma::vec & res,arma::vec & der,arma::vec & grad,const int n,const double tau,const double h,const double n1,const double h1,const double h3)110 void updatePara(const arma::mat& Z, const arma::vec& res, arma::vec& der, arma::vec& grad, const int n, const double tau, const double h,
111 const double n1, const double h1, const double h3) {
112 for (int i = 0; i < n; i++) {
113 double cur = res(i);
114 if (cur <= -h) {
115 der(i) = 1 - tau;
116 } else if (cur < h) {
117 der(i) = 0.5 - tau - 0.75 * h1 * cur + 0.25 * h3 * cur * cur * cur;
118 } else {
119 der(i) = -tau;
120 }
121 }
122 grad = n1 * Z.t() * der;
123 }
124
125 // [[Rcpp::export]]
updateTrian(const arma::mat & Z,const arma::vec & res,arma::vec & der,arma::vec & grad,const int n,const double tau,const double h,const double n1,const double h1,const double h2)126 void updateTrian(const arma::mat& Z, const arma::vec& res, arma::vec& der, arma::vec& grad, const int n, const double tau, const double h,
127 const double n1, const double h1, const double h2) {
128 for (int i = 0; i < n; i++) {
129 double cur = res(i);
130 if (cur <= -h) {
131 der(i) = 1 - tau;
132 } else if (cur < 0) {
133 der(i) = 0.5 - tau - h1 * cur - 0.5 * h2 * cur * cur;
134 } else if (cur < h) {
135 der(i) = 0.5 - tau - h1 * cur + 0.5 * h2 * cur * cur;
136 } else {
137 der(i) = -tau;
138 }
139 }
140 grad = n1 * Z.t() * der;
141 }
142
143 // Low-dimensional conquer: estimation
144 // [[Rcpp::export]]
smqrGauss(const arma::mat & X,arma::vec Y,const double tau=0.5,double h=0.0,const double constTau=1.345,const double tol=0.0001,const int iteMax=5000)145 Rcpp::List smqrGauss(const arma::mat& X, arma::vec Y, const double tau = 0.5, double h = 0.0, const double constTau = 1.345,
146 const double tol = 0.0001, const int iteMax = 5000) {
147 const int n = X.n_rows;
148 const int p = X.n_cols;
149 if (h <= 0.05) {
150 h = std::max(std::pow((std::log(n) + p) / n, 0.4), 0.05);
151 }
152 const double n1 = 1.0 / n;
153 const double h1 = 1.0 / h;
154 arma::rowvec mx = arma::mean(X, 0);
155 arma::vec sx1 = 1.0 / arma::stddev(X, 0, 0).t();
156 arma::mat Z = arma::join_rows(arma::ones(n), standardize(X, mx, sx1, p));
157 double my = arma::mean(Y);
158 Y -= my;
159 arma::vec der(n);
160 arma::vec gradOld(p + 1), gradNew(p + 1);
161 arma::vec beta = huberReg(Z, Y, tau, der, gradOld, gradNew, n, p, n1, tol, constTau, iteMax);
162 arma::vec quant = {tau};
163 beta(0) = arma::as_scalar(arma::quantile(Y - Z.cols(1, p) * beta.rows(1, p), quant));
164 arma::vec res = Y - Z * beta;
165 updateGauss(Z, res, der, gradOld, tau, n1, h1);
166 beta -= gradOld;
167 arma::vec betaDiff = -gradOld;
168 res -= Z * betaDiff;
169 updateGauss(Z, res, der, gradNew, tau, n1, h1);
170 arma::vec gradDiff = gradNew - gradOld;
171 int ite = 1;
172 while (arma::norm(gradNew, "inf") > tol && ite <= iteMax) {
173 double alpha = 1.0;
174 double cross = arma::as_scalar(betaDiff.t() * gradDiff);
175 if (cross > 0) {
176 double a1 = cross / arma::as_scalar(gradDiff.t() * gradDiff);
177 double a2 = arma::as_scalar(betaDiff.t() * betaDiff) / cross;
178 alpha = std::min(std::min(a1, a2), 100.0);
179 }
180 gradOld = gradNew;
181 betaDiff = -alpha * gradNew;
182 beta += betaDiff;
183 res -= Z * betaDiff;
184 updateGauss(Z, res, der, gradNew, tau, n1, h1);
185 gradDiff = gradNew - gradOld;
186 ite++;
187 }
188 beta.rows(1, p) %= sx1;
189 beta(0) += my - arma::as_scalar(mx * beta.rows(1, p));
190 return Rcpp::List::create(Rcpp::Named("coeff") = beta, Rcpp::Named("ite") = ite, Rcpp::Named("residual") = res, Rcpp::Named("bandwidth") = h);
191 }
192
193 // [[Rcpp::export]]
smqrGaussNsd(const arma::mat & Z,const arma::vec & Y,const double tau=0.5,double h=0.0,const double constTau=1.345,const double tol=0.0001,const int iteMax=5000)194 Rcpp::List smqrGaussNsd(const arma::mat& Z, const arma::vec& Y, const double tau = 0.5, double h = 0.0, const double constTau = 1.345,
195 const double tol = 0.0001, const int iteMax = 5000) {
196 const int n = Z.n_rows;
197 const int p = Z.n_cols - 1;
198 if (h <= 0.05) {
199 h = std::max(std::pow((std::log(n) + p) / n, 0.4), 0.05);
200 }
201 const double n1 = 1.0 / n;
202 const double h1 = 1.0 / h;
203 arma::vec der(n);
204 arma::vec gradOld(p + 1), gradNew(p + 1);
205 arma::vec beta = huberReg(Z, Y, tau, der, gradOld, gradNew, n, p, n1, tol, constTau, iteMax);
206 arma::vec quant = {tau};
207 beta(0) = arma::as_scalar(arma::quantile(Y - Z.cols(1, p) * beta.rows(1, p), quant));
208 arma::vec res = Y - Z * beta;
209 updateGauss(Z, res, der, gradOld, tau, n1, h1);
210 beta -= gradOld;
211 arma::vec betaDiff = -gradOld;
212 res -= Z * betaDiff;
213 updateGauss(Z, res, der, gradNew, tau, n1, h1);
214 arma::vec gradDiff = gradNew - gradOld;
215 int ite = 1;
216 while (arma::norm(gradNew, "inf") > tol && ite <= iteMax) {
217 double alpha = 1.0;
218 double cross = arma::as_scalar(betaDiff.t() * gradDiff);
219 if (cross > 0) {
220 double a1 = cross / arma::as_scalar(gradDiff.t() * gradDiff);
221 double a2 = arma::as_scalar(betaDiff.t() * betaDiff) / cross;
222 alpha = std::min(std::min(a1, a2), 100.0);
223 }
224 gradOld = gradNew;
225 betaDiff = -alpha * gradNew;
226 beta += betaDiff;
227 res -= Z * betaDiff;
228 updateGauss(Z, res, der, gradNew, tau, n1, h1);
229 gradDiff = gradNew - gradOld;
230 ite++;
231 }
232 return Rcpp::List::create(Rcpp::Named("coeff") = beta, Rcpp::Named("ite") = ite, Rcpp::Named("residual") = res, Rcpp::Named("bandwidth") = h);
233 }
234
235 // [[Rcpp::export]]
smqrGaussIni(const arma::mat & X,arma::vec Y,const arma::vec & betaHat,const int p,const double tau=0.5,double h=0.0,const double tol=0.0001,const int iteMax=5000)236 arma::vec smqrGaussIni(const arma::mat& X, arma::vec Y, const arma::vec& betaHat, const int p, const double tau = 0.5, double h = 0.0,
237 const double tol = 0.0001, const int iteMax = 5000) {
238 const int n = X.n_rows;
239 if (h <= 0.05) {
240 h = std::max(std::pow((std::log(n) + p) / n, 0.4), 0.05);
241 }
242 const double n1 = 1.0 / n;
243 const double h1 = 1.0 / h;
244 arma::rowvec mx = arma::mean(X, 0);
245 arma::vec sx1 = 1.0 / arma::stddev(X, 0, 0).t();
246 arma::mat Z = arma::join_rows(arma::ones(n), standardize(X, mx, sx1, p));
247 double my = arma::mean(Y);
248 Y -= my;
249 arma::vec der(n);
250 arma::vec gradOld(p + 1), gradNew(p + 1);
251 arma::vec beta = betaHat;
252 arma::vec res = Y - Z * beta;
253 updateGauss(Z, res, der, gradOld, tau, n1, h1);
254 beta -= gradOld;
255 arma::vec betaDiff = -gradOld;
256 res -= Z * betaDiff;
257 updateGauss(Z, res, der, gradNew, tau, n1, h1);
258 arma::vec gradDiff = gradNew - gradOld;
259 int ite = 1;
260 while (arma::norm(gradNew, "inf") > tol && ite <= iteMax) {
261 double alpha = 1.0;
262 double cross = arma::as_scalar(betaDiff.t() * gradDiff);
263 if (cross > 0) {
264 double a1 = cross / arma::as_scalar(gradDiff.t() * gradDiff);
265 double a2 = arma::as_scalar(betaDiff.t() * betaDiff) / cross;
266 alpha = std::min(std::min(a1, a2), 100.0);
267 }
268 gradOld = gradNew;
269 betaDiff = -alpha * gradNew;
270 beta += betaDiff;
271 res -= Z * betaDiff;
272 updateGauss(Z, res, der, gradNew, tau, n1, h1);
273 gradDiff = gradNew - gradOld;
274 ite++;
275 }
276 beta.rows(1, p) %= sx1;
277 beta(0) += my - arma::as_scalar(mx * beta.rows(1, p));
278 return beta;
279 }
280
281 // [[Rcpp::export]]
smqrGaussIniWeight(const arma::mat & X,arma::vec Y,const arma::vec & weight,const arma::vec & betaHat,const int p,const double tau=0.5,double h=0.0,const double tol=0.0001,const int iteMax=5000)282 arma::vec smqrGaussIniWeight(const arma::mat& X, arma::vec Y, const arma::vec& weight, const arma::vec& betaHat, const int p, const double tau = 0.5,
283 double h = 0.0, const double tol = 0.0001, const int iteMax = 5000) {
284 const int n = X.n_rows;
285 if (h <= 0.05) {
286 h = std::max(std::pow((std::log(n) + p) / n, 0.4), 0.05);
287 }
288 const double n1 = 1.0 / n;
289 const double h1 = 1.0 / h;
290 arma::rowvec mx = arma::mean(X, 0);
291 arma::vec sx1 = 1.0 / arma::stddev(X, 0, 0).t();
292 arma::mat Z = arma::join_rows(arma::ones(n), standardize(X, mx, sx1, p));
293 double my = arma::mean(Y);
294 Y -= my;
295 arma::vec der(n);
296 arma::vec gradOld(p + 1), gradNew(p + 1);
297 arma::vec beta = betaHat;
298 arma::vec res = Y - Z * beta;
299 updateGaussWeight(Z, weight, res, der, gradOld, tau, n1, h1);
300 beta -= gradOld;
301 arma::vec betaDiff = -gradOld;
302 res -= Z * betaDiff;
303 updateGaussWeight(Z, weight, res, der, gradNew, tau, n1, h1);
304 arma::vec gradDiff = gradNew - gradOld;
305 int ite = 1;
306 while (arma::norm(gradNew, "inf") > tol && ite <= iteMax) {
307 double alpha = 1.0;
308 double cross = arma::as_scalar(betaDiff.t() * gradDiff);
309 if (cross > 0) {
310 double a1 = cross / arma::as_scalar(gradDiff.t() * gradDiff);
311 double a2 = arma::as_scalar(betaDiff.t() * betaDiff) / cross;
312 alpha = std::min(std::min(a1, a2), 100.0);
313 }
314 gradOld = gradNew;
315 betaDiff = -alpha * gradNew;
316 beta += betaDiff;
317 res -= Z * betaDiff;
318 updateGaussWeight(Z, weight, res, der, gradNew, tau, n1, h1);
319 gradDiff = gradNew - gradOld;
320 ite++;
321 }
322 beta.rows(1, p) %= sx1;
323 beta(0) += my - arma::as_scalar(mx * beta.rows(1, p));
324 return beta;
325 }
326
327 // [[Rcpp::export]]
smqrLogistic(const arma::mat & X,arma::vec Y,const double tau=0.5,double h=0.0,const double constTau=1.345,const double tol=0.0001,const int iteMax=5000)328 Rcpp::List smqrLogistic(const arma::mat& X, arma::vec Y, const double tau = 0.5, double h = 0.0, const double constTau = 1.345,
329 const double tol = 0.0001, const int iteMax = 5000) {
330 const int n = X.n_rows;
331 const int p = X.n_cols;
332 if (h <= 0.05) {
333 h = std::max(std::pow((std::log(n) + p) / n, 0.4), 0.05);
334 }
335 const double n1 = 1.0 / n;
336 const double h1 = 1.0 / h;
337 arma::rowvec mx = arma::mean(X, 0);
338 arma::vec sx1 = 1.0 / arma::stddev(X, 0, 0).t();
339 arma::mat Z = arma::join_rows(arma::ones(n), standardize(X, mx, sx1, p));
340 double my = arma::mean(Y);
341 Y -= my;
342 arma::vec der(n);
343 arma::vec gradOld(p + 1), gradNew(p + 1);
344 arma::vec beta = huberReg(Z, Y, tau, der, gradOld, gradNew, n, p, n1, tol, constTau, iteMax);
345 arma::vec quant = {tau};
346 beta(0) = arma::as_scalar(arma::quantile(Y - Z.cols(1, p) * beta.rows(1, p), quant));
347 arma::vec res = Y - Z * beta;
348 updateLogistic(Z, res, der, gradOld, tau, n1, h1);
349 beta -= gradOld;
350 arma::vec betaDiff = -gradOld;
351 res -= Z * betaDiff;
352 updateLogistic(Z, res, der, gradNew, tau, n1, h1);
353 arma::vec gradDiff = gradNew - gradOld;
354 int ite = 1;
355 while (arma::norm(gradNew, "inf") > tol && ite <= iteMax) {
356 double alpha = 1.0;
357 double cross = arma::as_scalar(betaDiff.t() * gradDiff);
358 if (cross > 0) {
359 double a1 = cross / arma::as_scalar(gradDiff.t() * gradDiff);
360 double a2 = arma::as_scalar(betaDiff.t() * betaDiff) / cross;
361 alpha = std::min(std::min(a1, a2), 100.0);
362 }
363 gradOld = gradNew;
364 betaDiff = -alpha * gradNew;
365 beta += betaDiff;
366 res -= Z * betaDiff;
367 updateLogistic(Z, res, der, gradNew, tau, n1, h1);
368 gradDiff = gradNew - gradOld;
369 ite++;
370 }
371 beta.rows(1, p) %= sx1;
372 beta(0) += my - arma::as_scalar(mx * beta.rows(1, p));
373 return Rcpp::List::create(Rcpp::Named("coeff") = beta, Rcpp::Named("ite") = ite, Rcpp::Named("residual") = res, Rcpp::Named("bandwidth") = h);
374 }
375
376 // [[Rcpp::export]]
smqrLogisticNsd(const arma::mat & Z,const arma::vec & Y,const double tau=0.5,double h=0.0,const double constTau=1.345,const double tol=0.0001,const int iteMax=5000)377 Rcpp::List smqrLogisticNsd(const arma::mat& Z, const arma::vec& Y, const double tau = 0.5, double h = 0.0, const double constTau = 1.345,
378 const double tol = 0.0001, const int iteMax = 5000) {
379 const int n = Z.n_rows;
380 const int p = Z.n_cols - 1;
381 if (h <= 0.05) {
382 h = std::max(std::pow((std::log(n) + p) / n, 0.4), 0.05);
383 }
384 const double n1 = 1.0 / n;
385 const double h1 = 1.0 / h;
386 arma::vec der(n);
387 arma::vec gradOld(p + 1), gradNew(p + 1);
388 arma::vec beta = huberReg(Z, Y, tau, der, gradOld, gradNew, n, p, n1, tol, constTau, iteMax);
389 arma::vec quant = {tau};
390 beta(0) = arma::as_scalar(arma::quantile(Y - Z.cols(1, p) * beta.rows(1, p), quant));
391 arma::vec res = Y - Z * beta;
392 updateLogistic(Z, res, der, gradOld, tau, n1, h1);
393 beta -= gradOld;
394 arma::vec betaDiff = -gradOld;
395 res -= Z * betaDiff;
396 updateLogistic(Z, res, der, gradNew, tau, n1, h1);
397 arma::vec gradDiff = gradNew - gradOld;
398 int ite = 1;
399 while (arma::norm(gradNew, "inf") > tol && ite <= iteMax) {
400 double alpha = 1.0;
401 double cross = arma::as_scalar(betaDiff.t() * gradDiff);
402 if (cross > 0) {
403 double a1 = cross / arma::as_scalar(gradDiff.t() * gradDiff);
404 double a2 = arma::as_scalar(betaDiff.t() * betaDiff) / cross;
405 alpha = std::min(std::min(a1, a2), 100.0);
406 }
407 gradOld = gradNew;
408 betaDiff = -alpha * gradNew;
409 beta += betaDiff;
410 res -= Z * betaDiff;
411 updateLogistic(Z, res, der, gradNew, tau, n1, h1);
412 gradDiff = gradNew - gradOld;
413 ite++;
414 }
415 return Rcpp::List::create(Rcpp::Named("coeff") = beta, Rcpp::Named("ite") = ite, Rcpp::Named("residual") = res, Rcpp::Named("bandwidth") = h);
416 }
417
418 // [[Rcpp::export]]
smqrLogisticIni(const arma::mat & X,arma::vec Y,const arma::vec & betaHat,const int p,const double tau=0.5,double h=0.0,const double tol=0.0001,const int iteMax=5000)419 arma::vec smqrLogisticIni(const arma::mat& X, arma::vec Y, const arma::vec& betaHat, const int p, const double tau = 0.5, double h = 0.0,
420 const double tol = 0.0001, const int iteMax = 5000) {
421 const int n = X.n_rows;
422 if (h <= 0.05) {
423 h = std::max(std::pow((std::log(n) + p) / n, 0.4), 0.05);
424 }
425 const double n1 = 1.0 / n;
426 const double h1 = 1.0 / h;
427 arma::rowvec mx = arma::mean(X, 0);
428 arma::vec sx1 = 1.0 / arma::stddev(X, 0, 0).t();
429 arma::mat Z = arma::join_rows(arma::ones(n), standardize(X, mx, sx1, p));
430 double my = arma::mean(Y);
431 Y -= my;
432 arma::vec der(n);
433 arma::vec gradOld(p + 1), gradNew(p + 1);
434 arma::vec beta = betaHat;
435 arma::vec res = Y - Z * beta;
436 updateLogistic(Z, res, der, gradOld, tau, n1, h1);
437 beta -= gradOld;
438 arma::vec betaDiff = -gradOld;
439 res -= Z * betaDiff;
440 updateLogistic(Z, res, der, gradNew, tau, n1, h1);
441 arma::vec gradDiff = gradNew - gradOld;
442 int ite = 1;
443 while (arma::norm(gradNew, "inf") > tol && ite <= iteMax) {
444 double alpha = 1.0;
445 double cross = arma::as_scalar(betaDiff.t() * gradDiff);
446 if (cross > 0) {
447 double a1 = cross / arma::as_scalar(gradDiff.t() * gradDiff);
448 double a2 = arma::as_scalar(betaDiff.t() * betaDiff) / cross;
449 alpha = std::min(std::min(a1, a2), 100.0);
450 }
451 gradOld = gradNew;
452 betaDiff = -alpha * gradNew;
453 beta += betaDiff;
454 res -= Z * betaDiff;
455 updateLogistic(Z, res, der, gradNew, tau, n1, h1);
456 gradDiff = gradNew - gradOld;
457 ite++;
458 }
459 beta.rows(1, p) %= sx1;
460 beta(0) += my - arma::as_scalar(mx * beta.rows(1, p));
461 return beta;
462 }
463
464 // [[Rcpp::export]]
smqrUnif(const arma::mat & X,arma::vec Y,const double tau=0.5,double h=0.0,const double constTau=1.345,const double tol=0.0001,const int iteMax=5000)465 Rcpp::List smqrUnif(const arma::mat& X, arma::vec Y, const double tau = 0.5, double h = 0.0, const double constTau = 1.345,
466 const double tol = 0.0001, const int iteMax = 5000) {
467 const int n = X.n_rows;
468 const int p = X.n_cols;
469 if (h <= 0.05) {
470 h = std::max(std::pow((std::log(n) + p) / n, 0.4), 0.05);
471 }
472 const double n1 = 1.0 / n;
473 const double h1 = 1.0 / h;
474 arma::rowvec mx = arma::mean(X, 0);
475 arma::vec sx1 = 1.0 / arma::stddev(X, 0, 0).t();
476 arma::mat Z = arma::join_rows(arma::ones(n), standardize(X, mx, sx1, p));
477 double my = arma::mean(Y);
478 Y -= my;
479 arma::vec der(n);
480 arma::vec gradOld(p + 1), gradNew(p + 1);
481 arma::vec beta = huberReg(Z, Y, tau, der, gradOld, gradNew, n, p, n1, tol, constTau, iteMax);
482 arma::vec quant = {tau};
483 beta(0) = arma::as_scalar(arma::quantile(Y - Z.cols(1, p) * beta.rows(1, p), quant));
484 arma::vec res = Y - Z * beta;
485 updateUnif(Z, res, der, gradOld, n, tau, h, n1, h1);
486 beta -= gradOld;
487 arma::vec betaDiff = -gradOld;
488 res -= Z * betaDiff;
489 updateUnif(Z, res, der, gradNew, n, tau, h, n1, h1);
490 arma::vec gradDiff = gradNew - gradOld;
491 int ite = 1;
492 while (arma::norm(gradNew, "inf") > tol && ite <= iteMax) {
493 double alpha = 1.0;
494 double cross = arma::as_scalar(betaDiff.t() * gradDiff);
495 if (cross > 0) {
496 double a1 = cross / arma::as_scalar(gradDiff.t() * gradDiff);
497 double a2 = arma::as_scalar(betaDiff.t() * betaDiff) / cross;
498 alpha = std::min(std::min(a1, a2), 100.0);
499 }
500 gradOld = gradNew;
501 betaDiff = -alpha * gradNew;
502 beta += betaDiff;
503 res -= Z * betaDiff;
504 updateUnif(Z, res, der, gradNew, n, tau, h, n1, h1);
505 gradDiff = gradNew - gradOld;
506 ite++;
507 }
508 beta.rows(1, p) %= sx1;
509 beta(0) += my - arma::as_scalar(mx * beta.rows(1, p));
510 return Rcpp::List::create(Rcpp::Named("coeff") = beta, Rcpp::Named("ite") = ite, Rcpp::Named("residual") = res, Rcpp::Named("bandwidth") = h);
511 }
512
513 // [[Rcpp::export]]
smqrUnifNsd(const arma::mat & Z,const arma::vec & Y,const double tau=0.5,double h=0.0,const double constTau=1.345,const double tol=0.0001,const int iteMax=5000)514 Rcpp::List smqrUnifNsd(const arma::mat& Z, const arma::vec& Y, const double tau = 0.5, double h = 0.0, const double constTau = 1.345,
515 const double tol = 0.0001, const int iteMax = 5000) {
516 const int n = Z.n_rows;
517 const int p = Z.n_cols - 1;
518 if (h <= 0.05) {
519 h = std::max(std::pow((std::log(n) + p) / n, 0.4), 0.05);
520 }
521 const double n1 = 1.0 / n;
522 const double h1 = 1.0 / h;
523 arma::vec der(n);
524 arma::vec gradOld(p + 1), gradNew(p + 1);
525 arma::vec beta = huberReg(Z, Y, tau, der, gradOld, gradNew, n, p, n1, tol, constTau, iteMax);
526 arma::vec quant = {tau};
527 beta(0) = arma::as_scalar(arma::quantile(Y - Z.cols(1, p) * beta.rows(1, p), quant));
528 arma::vec res = Y - Z * beta;
529 updateUnif(Z, res, der, gradOld, n, tau, h, n1, h1);
530 beta -= gradOld;
531 arma::vec betaDiff = -gradOld;
532 res -= Z * betaDiff;
533 updateUnif(Z, res, der, gradNew, n, tau, h, n1, h1);
534 arma::vec gradDiff = gradNew - gradOld;
535 int ite = 1;
536 while (arma::norm(gradNew, "inf") > tol && ite <= iteMax) {
537 double alpha = 1.0;
538 double cross = arma::as_scalar(betaDiff.t() * gradDiff);
539 if (cross > 0) {
540 double a1 = cross / arma::as_scalar(gradDiff.t() * gradDiff);
541 double a2 = arma::as_scalar(betaDiff.t() * betaDiff) / cross;
542 alpha = std::min(std::min(a1, a2), 100.0);
543 }
544 gradOld = gradNew;
545 betaDiff = -alpha * gradNew;
546 beta += betaDiff;
547 res -= Z * betaDiff;
548 updateUnif(Z, res, der, gradNew, n, tau, h, n1, h1);
549 gradDiff = gradNew - gradOld;
550 ite++;
551 }
552 return Rcpp::List::create(Rcpp::Named("coeff") = beta, Rcpp::Named("ite") = ite, Rcpp::Named("residual") = res, Rcpp::Named("bandwidth") = h);
553 }
554
555 // [[Rcpp::export]]
smqrUnifIni(const arma::mat & X,arma::vec Y,const arma::vec & betaHat,const int p,const double tau=0.5,double h=0.0,const double tol=0.0001,const int iteMax=5000)556 arma::vec smqrUnifIni(const arma::mat& X, arma::vec Y, const arma::vec& betaHat, const int p, const double tau = 0.5, double h = 0.0,
557 const double tol = 0.0001, const int iteMax = 5000) {
558 const int n = X.n_rows;
559 if (h <= 0.05) {
560 h = std::max(std::pow((std::log(n) + p) / n, 0.4), 0.05);
561 }
562 const double n1 = 1.0 / n;
563 const double h1 = 1.0 / h;
564 arma::rowvec mx = arma::mean(X, 0);
565 arma::vec sx1 = 1.0 / arma::stddev(X, 0, 0).t();
566 arma::mat Z = arma::join_rows(arma::ones(n), standardize(X, mx, sx1, p));
567 double my = arma::mean(Y);
568 Y -= my;
569 arma::vec der(n);
570 arma::vec gradOld(p + 1), gradNew(p + 1);
571 arma::vec beta = betaHat;
572 arma::vec res = Y - Z * beta;
573 updateUnif(Z, res, der, gradOld, n, tau, h, n1, h1);
574 beta -= gradOld;
575 arma::vec betaDiff = -gradOld;
576 res -= Z * betaDiff;
577 updateUnif(Z, res, der, gradNew, n, tau, h, n1, h1);
578 arma::vec gradDiff = gradNew - gradOld;
579 int ite = 1;
580 while (arma::norm(gradNew, "inf") > tol && ite <= iteMax) {
581 double alpha = 1.0;
582 double cross = arma::as_scalar(betaDiff.t() * gradDiff);
583 if (cross > 0) {
584 double a1 = cross / arma::as_scalar(gradDiff.t() * gradDiff);
585 double a2 = arma::as_scalar(betaDiff.t() * betaDiff) / cross;
586 alpha = std::min(std::min(a1, a2), 100.0);
587 }
588 gradOld = gradNew;
589 betaDiff = -alpha * gradNew;
590 beta += betaDiff;
591 res -= Z * betaDiff;
592 updateUnif(Z, res, der, gradNew, n, tau, h, n1, h1);
593 gradDiff = gradNew - gradOld;
594 ite++;
595 }
596 beta.rows(1, p) %= sx1;
597 beta(0) += my - arma::as_scalar(mx * beta.rows(1, p));
598 return beta;
599 }
600
601 // [[Rcpp::export]]
smqrPara(const arma::mat & X,arma::vec Y,const double tau=0.5,double h=0.0,const double constTau=1.345,const double tol=0.0001,const int iteMax=5000)602 Rcpp::List smqrPara(const arma::mat& X, arma::vec Y, const double tau = 0.5, double h = 0.0, const double constTau = 1.345,
603 const double tol = 0.0001, const int iteMax = 5000) {
604 const int n = X.n_rows;
605 const int p = X.n_cols;
606 if (h <= 0.05) {
607 h = std::max(std::pow((std::log(n) + p) / n, 0.4), 0.05);
608 }
609 const double n1 = 1.0 / n;
610 const double h1 = 1.0 / h, h3 = 1.0 / (h * h * h);
611 arma::rowvec mx = arma::mean(X, 0);
612 arma::vec sx1 = 1.0 / arma::stddev(X, 0, 0).t();
613 arma::mat Z = arma::join_rows(arma::ones(n), standardize(X, mx, sx1, p));
614 double my = arma::mean(Y);
615 Y -= my;
616 arma::vec der(n);
617 arma::vec gradOld(p + 1), gradNew(p + 1);
618 arma::vec beta = huberReg(Z, Y, tau, der, gradOld, gradNew, n, p, n1, tol, constTau, iteMax);
619 arma::vec quant = {tau};
620 beta(0) = arma::as_scalar(arma::quantile(Y - Z.cols(1, p) * beta.rows(1, p), quant));
621 arma::vec res = Y - Z * beta;
622 updatePara(Z, res, der, gradOld, n, tau, h, n1, h1, h3);
623 beta -= gradOld;
624 arma::vec betaDiff = -gradOld;
625 res -= Z * betaDiff;
626 updatePara(Z, res, der, gradNew, n, tau, h, n1, h1, h3);
627 arma::vec gradDiff = gradNew - gradOld;
628 int ite = 1;
629 while (arma::norm(gradNew, "inf") > tol && ite <= iteMax) {
630 double alpha = 1.0;
631 double cross = arma::as_scalar(betaDiff.t() * gradDiff);
632 if (cross > 0) {
633 double a1 = cross / arma::as_scalar(gradDiff.t() * gradDiff);
634 double a2 = arma::as_scalar(betaDiff.t() * betaDiff) / cross;
635 alpha = std::min(std::min(a1, a2), 100.0);
636 }
637 gradOld = gradNew;
638 betaDiff = -alpha * gradNew;
639 beta += betaDiff;
640 res -= Z * betaDiff;
641 updatePara(Z, res, der, gradNew, n, tau, h, n1, h1, h3);
642 gradDiff = gradNew - gradOld;
643 ite++;
644 }
645 beta.rows(1, p) %= sx1;
646 beta(0) += my - arma::as_scalar(mx * beta.rows(1, p));
647 return Rcpp::List::create(Rcpp::Named("coeff") = beta, Rcpp::Named("ite") = ite, Rcpp::Named("residual") = res, Rcpp::Named("bandwidth") = h);
648 }
649
650 // [[Rcpp::export]]
smqrParaNsd(const arma::mat & Z,const arma::vec & Y,const double tau=0.5,double h=0.0,const double constTau=1.345,const double tol=0.0001,const int iteMax=5000)651 Rcpp::List smqrParaNsd(const arma::mat& Z, const arma::vec& Y, const double tau = 0.5, double h = 0.0, const double constTau = 1.345,
652 const double tol = 0.0001, const int iteMax = 5000) {
653 const int n = Z.n_rows;
654 const int p = Z.n_cols - 1;
655 if (h <= 0.05) {
656 h = std::max(std::pow((std::log(n) + p) / n, 0.4), 0.05);
657 }
658 const double n1 = 1.0 / n;
659 const double h1 = 1.0 / h, h3 = 1.0 / (h * h * h);
660 arma::vec der(n);
661 arma::vec gradOld(p + 1), gradNew(p + 1);
662 arma::vec beta = huberReg(Z, Y, tau, der, gradOld, gradNew, n, p, n1, tol, constTau, iteMax);
663 arma::vec quant = {tau};
664 beta(0) = arma::as_scalar(arma::quantile(Y - Z.cols(1, p) * beta.rows(1, p), quant));
665 arma::vec res = Y - Z * beta;
666 updatePara(Z, res, der, gradOld, n, tau, h, n1, h1, h3);
667 beta -= gradOld;
668 arma::vec betaDiff = -gradOld;
669 res -= Z * betaDiff;
670 updatePara(Z, res, der, gradNew, n, tau, h, n1, h1, h3);
671 arma::vec gradDiff = gradNew - gradOld;
672 int ite = 1;
673 while (arma::norm(gradNew, "inf") > tol && ite <= iteMax) {
674 double alpha = 1.0;
675 double cross = arma::as_scalar(betaDiff.t() * gradDiff);
676 if (cross > 0) {
677 double a1 = cross / arma::as_scalar(gradDiff.t() * gradDiff);
678 double a2 = arma::as_scalar(betaDiff.t() * betaDiff) / cross;
679 alpha = std::min(std::min(a1, a2), 100.0);
680 }
681 gradOld = gradNew;
682 betaDiff = -alpha * gradNew;
683 beta += betaDiff;
684 res -= Z * betaDiff;
685 updatePara(Z, res, der, gradNew, n, tau, h, n1, h1, h3);
686 gradDiff = gradNew - gradOld;
687 ite++;
688 }
689 return Rcpp::List::create(Rcpp::Named("coeff") = beta, Rcpp::Named("ite") = ite, Rcpp::Named("residual") = res, Rcpp::Named("bandwidth") = h);
690 }
691
692 // [[Rcpp::export]]
smqrParaIni(const arma::mat & X,arma::vec Y,const arma::vec & betaHat,const int p,const double tau=0.5,double h=0.0,const double tol=0.0001,const int iteMax=5000)693 arma::vec smqrParaIni(const arma::mat& X, arma::vec Y, const arma::vec& betaHat, const int p, const double tau = 0.5, double h = 0.0,
694 const double tol = 0.0001, const int iteMax = 5000) {
695 const int n = X.n_rows;
696 if (h <= 0.05) {
697 h = std::max(std::pow((std::log(n) + p) / n, 0.4), 0.05);
698 }
699 const double n1 = 1.0 / n;
700 const double h1 = 1.0 / h, h3 = 1.0 / (h * h * h);
701 arma::rowvec mx = arma::mean(X, 0);
702 arma::vec sx1 = 1.0 / arma::stddev(X, 0, 0).t();
703 arma::mat Z = arma::join_rows(arma::ones(n), standardize(X, mx, sx1, p));
704 double my = arma::mean(Y);
705 Y -= my;
706 arma::vec der(n);
707 arma::vec gradOld(p + 1), gradNew(p + 1);
708 arma::vec beta = betaHat;
709 arma::vec res = Y - Z * beta;
710 updatePara(Z, res, der, gradOld, n, tau, h, n1, h1, h3);
711 beta -= gradOld;
712 arma::vec betaDiff = -gradOld;
713 res -= Z * betaDiff;
714 updatePara(Z, res, der, gradNew, n, tau, h, n1, h1, h3);
715 arma::vec gradDiff = gradNew - gradOld;
716 int ite = 1;
717 while (arma::norm(gradNew, "inf") > tol && ite <= iteMax) {
718 double alpha = 1.0;
719 double cross = arma::as_scalar(betaDiff.t() * gradDiff);
720 if (cross > 0) {
721 double a1 = cross / arma::as_scalar(gradDiff.t() * gradDiff);
722 double a2 = arma::as_scalar(betaDiff.t() * betaDiff) / cross;
723 alpha = std::min(std::min(a1, a2), 100.0);
724 }
725 gradOld = gradNew;
726 betaDiff = -alpha * gradNew;
727 beta += betaDiff;
728 res -= Z * betaDiff;
729 updatePara(Z, res, der, gradNew, n, tau, h, n1, h1, h3);
730 gradDiff = gradNew - gradOld;
731 ite++;
732 }
733 beta.rows(1, p) %= sx1;
734 beta(0) += my - arma::as_scalar(mx * beta.rows(1, p));
735 return beta;
736 }
737
738 // [[Rcpp::export]]
smqrTrian(const arma::mat & X,arma::vec Y,const double tau=0.5,double h=0.0,const double constTau=1.345,const double tol=0.0001,const int iteMax=5000)739 Rcpp::List smqrTrian(const arma::mat& X, arma::vec Y, const double tau = 0.5, double h = 0.0, const double constTau = 1.345,
740 const double tol = 0.0001, const int iteMax = 5000) {
741 const int n = X.n_rows;
742 const int p = X.n_cols;
743 if (h <= 0.05) {
744 h = std::max(std::pow((std::log(n) + p) / n, 0.4), 0.05);
745 }
746 const double n1 = 1.0 / n;
747 const double h1 = 1.0 / h, h2 = 1.0 / (h * h);
748 arma::rowvec mx = arma::mean(X, 0);
749 arma::vec sx1 = 1.0 / arma::stddev(X, 0, 0).t();
750 arma::mat Z = arma::join_rows(arma::ones(n), standardize(X, mx, sx1, p));
751 double my = arma::mean(Y);
752 Y -= my;
753 arma::vec der(n);
754 arma::vec gradOld(p + 1), gradNew(p + 1);
755 arma::vec beta = huberReg(Z, Y, tau, der, gradOld, gradNew, n, p, n1, tol, constTau, iteMax);
756 arma::vec quant = {tau};
757 beta(0) = arma::as_scalar(arma::quantile(Y - Z.cols(1, p) * beta.rows(1, p), quant));
758 arma::vec res = Y - Z * beta;
759 updateTrian(Z, res, der, gradOld, n, tau, h, n1, h1, h2);
760 beta -= gradOld;
761 arma::vec betaDiff = -gradOld;
762 res -= Z * betaDiff;
763 updateTrian(Z, res, der, gradNew, n, tau, h, n1, h1, h2);
764 arma::vec gradDiff = gradNew - gradOld;
765 int ite = 1;
766 while (arma::norm(gradNew, "inf") > tol && ite <= iteMax) {
767 double alpha = 1.0;
768 double cross = arma::as_scalar(betaDiff.t() * gradDiff);
769 if (cross > 0) {
770 double a1 = cross / arma::as_scalar(gradDiff.t() * gradDiff);
771 double a2 = arma::as_scalar(betaDiff.t() * betaDiff) / cross;
772 alpha = std::min(std::min(a1, a2), 100.0);
773 }
774 gradOld = gradNew;
775 betaDiff = -alpha * gradNew;
776 beta += betaDiff;
777 res -= Z * betaDiff;
778 updateTrian(Z, res, der, gradNew, n, tau, h, n1, h1, h2);
779 gradDiff = gradNew - gradOld;
780 ite++;
781 }
782 beta.rows(1, p) %= sx1;
783 beta(0) += my - arma::as_scalar(mx * beta.rows(1, p));
784 return Rcpp::List::create(Rcpp::Named("coeff") = beta, Rcpp::Named("ite") = ite, Rcpp::Named("residual") = res, Rcpp::Named("bandwidth") = h);
785 }
786
787 // [[Rcpp::export]]
smqrTrianNsd(const arma::mat & Z,const arma::vec & Y,const double tau=0.5,double h=0.0,const double constTau=1.345,const double tol=0.0001,const int iteMax=5000)788 Rcpp::List smqrTrianNsd(const arma::mat& Z, const arma::vec& Y, const double tau = 0.5, double h = 0.0, const double constTau = 1.345,
789 const double tol = 0.0001, const int iteMax = 5000) {
790 const int n = Z.n_rows;
791 const int p = Z.n_cols - 1;
792 if (h <= 0.05) {
793 h = std::max(std::pow((std::log(n) + p) / n, 0.4), 0.05);
794 }
795 const double n1 = 1.0 / n;
796 const double h1 = 1.0 / h, h2 = 1.0 / (h * h);
797 arma::vec der(n);
798 arma::vec gradOld(p + 1), gradNew(p + 1);
799 arma::vec beta = huberReg(Z, Y, tau, der, gradOld, gradNew, n, p, n1, tol, constTau, iteMax);
800 arma::vec quant = {tau};
801 beta(0) = arma::as_scalar(arma::quantile(Y - Z.cols(1, p) * beta.rows(1, p), quant));
802 arma::vec res = Y - Z * beta;
803 updateTrian(Z, res, der, gradOld, n, tau, h, n1, h1, h2);
804 beta -= gradOld;
805 arma::vec betaDiff = -gradOld;
806 res -= Z * betaDiff;
807 updateTrian(Z, res, der, gradNew, n, tau, h, n1, h1, h2);
808 arma::vec gradDiff = gradNew - gradOld;
809 int ite = 1;
810 while (arma::norm(gradNew, "inf") > tol && ite <= iteMax) {
811 double alpha = 1.0;
812 double cross = arma::as_scalar(betaDiff.t() * gradDiff);
813 if (cross > 0) {
814 double a1 = cross / arma::as_scalar(gradDiff.t() * gradDiff);
815 double a2 = arma::as_scalar(betaDiff.t() * betaDiff) / cross;
816 alpha = std::min(std::min(a1, a2), 100.0);
817 }
818 gradOld = gradNew;
819 betaDiff = -alpha * gradNew;
820 beta += betaDiff;
821 res -= Z * betaDiff;
822 updateTrian(Z, res, der, gradNew, n, tau, h, n1, h1, h2);
823 gradDiff = gradNew - gradOld;
824 ite++;
825 }
826 return Rcpp::List::create(Rcpp::Named("coeff") = beta, Rcpp::Named("ite") = ite, Rcpp::Named("residual") = res, Rcpp::Named("bandwidth") = h);
827 }
828
829 // [[Rcpp::export]]
smqrTrianIni(const arma::mat & X,arma::vec Y,const arma::vec & betaHat,const int p,const double tau=0.5,double h=0.0,const double tol=0.0001,const int iteMax=5000)830 arma::vec smqrTrianIni(const arma::mat& X, arma::vec Y, const arma::vec& betaHat, const int p, const double tau = 0.5, double h = 0.0,
831 const double tol = 0.0001, const int iteMax = 5000) {
832 const int n = X.n_rows;
833 if (h <= 0.05) {
834 h = std::max(std::pow((std::log(n) + p) / n, 0.4), 0.05);
835 }
836 const double n1 = 1.0 / n;
837 const double h1 = 1.0 / h, h2 = 1.0 / (h * h);
838 arma::rowvec mx = arma::mean(X, 0);
839 arma::vec sx1 = 1.0 / arma::stddev(X, 0, 0).t();
840 arma::mat Z = arma::join_rows(arma::ones(n), standardize(X, mx, sx1, p));
841 double my = arma::mean(Y);
842 Y -= my;
843 arma::vec der(n);
844 arma::vec gradOld(p + 1), gradNew(p + 1);
845 arma::vec beta = betaHat;
846 arma::vec res = Y - Z * beta;
847 updateTrian(Z, res, der, gradOld, n, tau, h, n1, h1, h2);
848 beta -= gradOld;
849 arma::vec betaDiff = -gradOld;
850 res -= Z * betaDiff;
851 updateTrian(Z, res, der, gradNew, n, tau, h, n1, h1, h2);
852 arma::vec gradDiff = gradNew - gradOld;
853 int ite = 1;
854 while (arma::norm(gradNew, "inf") > tol && ite <= iteMax) {
855 double alpha = 1.0;
856 double cross = arma::as_scalar(betaDiff.t() * gradDiff);
857 if (cross > 0) {
858 double a1 = cross / arma::as_scalar(gradDiff.t() * gradDiff);
859 double a2 = arma::as_scalar(betaDiff.t() * betaDiff) / cross;
860 alpha = std::min(std::min(a1, a2), 100.0);
861 }
862 gradOld = gradNew;
863 betaDiff = -alpha * gradNew;
864 beta += betaDiff;
865 res -= Z * betaDiff;
866 updateTrian(Z, res, der, gradNew, n, tau, h, n1, h1, h2);
867 gradDiff = gradNew - gradOld;
868 ite++;
869 }
870 beta.rows(1, p) %= sx1;
871 beta(0) += my - arma::as_scalar(mx * beta.rows(1, p));
872 return beta;
873 }
874
875 // Global conquer process with a quantile grid
876 // [[Rcpp::export]]
smqrGaussProc(const arma::mat & X,arma::vec Y,const arma::vec tauSeq,double h=0.0,const double constTau=1.345,const double tol=0.0001,const int iteMax=5000)877 Rcpp::List smqrGaussProc(const arma::mat& X, arma::vec Y, const arma::vec tauSeq, double h = 0.0, const double constTau = 1.345,
878 const double tol = 0.0001, const int iteMax = 5000) {
879 const int n = X.n_rows;
880 const int p = X.n_cols;
881 const int m = tauSeq.size();
882 if (h <= 0.05) {
883 h = std::max(std::pow((std::log(n) + p) / n, 0.4), 0.05);
884 }
885 const double n1 = 1.0 / n;
886 const double h1 = 1.0 / h;
887 arma::rowvec mx = arma::mean(X, 0);
888 arma::vec sx1 = 1.0 / arma::stddev(X, 0, 0).t();
889 arma::mat Z = arma::join_rows(arma::ones(n), standardize(X, mx, sx1, p));
890 double my = arma::mean(Y);
891 Y -= my;
892 arma::vec der(n);
893 arma::vec gradOld(p + 1), gradNew(p + 1);
894 arma::mat betaProc(p + 1, m);
895 for (int i = 0; i < m; i++) {
896 double tau = tauSeq(i);
897 arma::vec beta = huberReg(Z, Y, tau, der, gradOld, gradNew, n, p, n1, tol, constTau, iteMax);
898 arma::vec quant = {tau};
899 beta(0) = arma::as_scalar(arma::quantile(Y - Z.cols(1, p) * beta.rows(1, p), quant));
900 arma::vec res = Y - Z * beta;
901 updateGauss(Z, res, der, gradOld, tau, n1, h1);
902 beta -= gradOld;
903 arma::vec betaDiff = -gradOld;
904 res -= Z * betaDiff;
905 updateGauss(Z, res, der, gradNew, tau, n1, h1);
906 arma::vec gradDiff = gradNew - gradOld;
907 int ite = 1;
908 while (arma::norm(gradNew, "inf") > tol && ite <= iteMax) {
909 double alpha = 1.0;
910 double cross = arma::as_scalar(betaDiff.t() * gradDiff);
911 if (cross > 0) {
912 double a1 = cross / arma::as_scalar(gradDiff.t() * gradDiff);
913 double a2 = arma::as_scalar(betaDiff.t() * betaDiff) / cross;
914 alpha = std::min(std::min(a1, a2), 100.0);
915 }
916 gradOld = gradNew;
917 betaDiff = -alpha * gradNew;
918 beta += betaDiff;
919 res -= Z * betaDiff;
920 updateGauss(Z, res, der, gradNew, tau, n1, h1);
921 gradDiff = gradNew - gradOld;
922 ite++;
923 }
924 betaProc.col(i) = beta;
925 }
926 betaProc.rows(1, p).each_col() %= sx1;
927 betaProc.row(0) += my - mx * betaProc.rows(1, p);
928 return Rcpp::List::create(Rcpp::Named("coeff") = betaProc, Rcpp::Named("bandwidth") = h);
929 }
930
931 // [[Rcpp::export]]
smqrLogisticProc(const arma::mat & X,arma::vec Y,const arma::vec tauSeq,double h=0.0,const double constTau=1.345,const double tol=0.0001,const int iteMax=5000)932 Rcpp::List smqrLogisticProc(const arma::mat& X, arma::vec Y, const arma::vec tauSeq, double h = 0.0, const double constTau = 1.345,
933 const double tol = 0.0001, const int iteMax = 5000) {
934 const int n = X.n_rows;
935 const int p = X.n_cols;
936 const int m = tauSeq.size();
937 if (h <= 0.05) {
938 h = std::max(std::pow((std::log(n) + p) / n, 0.4), 0.05);
939 }
940 const double n1 = 1.0 / n;
941 const double h1 = 1.0 / h;
942 arma::rowvec mx = arma::mean(X, 0);
943 arma::vec sx1 = 1.0 / arma::stddev(X, 0, 0).t();
944 arma::mat Z = arma::join_rows(arma::ones(n), standardize(X, mx, sx1, p));
945 double my = arma::mean(Y);
946 Y -= my;
947 arma::vec der(n);
948 arma::vec gradOld(p + 1), gradNew(p + 1);
949 arma::mat betaProc(p + 1, m);
950 for (int i = 0; i < m; i++) {
951 double tau = tauSeq(i);
952 arma::vec beta = huberReg(Z, Y, tau, der, gradOld, gradNew, n, p, n1, tol, constTau, iteMax);
953 arma::vec quant = {tau};
954 beta(0) = arma::as_scalar(arma::quantile(Y - Z.cols(1, p) * beta.rows(1, p), quant));
955 arma::vec res = Y - Z * beta;
956 updateLogistic(Z, res, der, gradOld, tau, n1, h1);
957 beta -= gradOld;
958 arma::vec betaDiff = -gradOld;
959 res -= Z * betaDiff;
960 updateLogistic(Z, res, der, gradNew, tau, n1, h1);
961 arma::vec gradDiff = gradNew - gradOld;
962 int ite = 1;
963 while (arma::norm(gradNew, "inf") > tol && ite <= iteMax) {
964 double alpha = 1.0;
965 double cross = arma::as_scalar(betaDiff.t() * gradDiff);
966 if (cross > 0) {
967 double a1 = cross / arma::as_scalar(gradDiff.t() * gradDiff);
968 double a2 = arma::as_scalar(betaDiff.t() * betaDiff) / cross;
969 alpha = std::min(std::min(a1, a2), 100.0);
970 }
971 gradOld = gradNew;
972 betaDiff = -alpha * gradNew;
973 beta += betaDiff;
974 res -= Z * betaDiff;
975 updateLogistic(Z, res, der, gradNew, tau, n1, h1);
976 gradDiff = gradNew - gradOld;
977 ite++;
978 }
979 betaProc.col(i) = beta;
980 }
981 betaProc.rows(1, p).each_col() %= sx1;
982 betaProc.row(0) += my - mx * betaProc.rows(1, p);
983 return Rcpp::List::create(Rcpp::Named("coeff") = betaProc, Rcpp::Named("bandwidth") = h);
984 }
985
986 // [[Rcpp::export]]
smqrUnifProc(const arma::mat & X,arma::vec Y,const arma::vec tauSeq,double h=0.0,const double constTau=1.345,const double tol=0.0001,const int iteMax=5000)987 Rcpp::List smqrUnifProc(const arma::mat& X, arma::vec Y, const arma::vec tauSeq, double h = 0.0, const double constTau = 1.345,
988 const double tol = 0.0001, const int iteMax = 5000) {
989 const int n = X.n_rows;
990 const int p = X.n_cols;
991 const int m = tauSeq.size();
992 if (h <= 0.05) {
993 h = std::max(std::pow((std::log(n) + p) / n, 0.4), 0.05);
994 }
995 const double n1 = 1.0 / n;
996 const double h1 = 1.0 / h;
997 arma::rowvec mx = arma::mean(X, 0);
998 arma::vec sx1 = 1.0 / arma::stddev(X, 0, 0).t();
999 arma::mat Z = arma::join_rows(arma::ones(n), standardize(X, mx, sx1, p));
1000 double my = arma::mean(Y);
1001 Y -= my;
1002 arma::vec der(n);
1003 arma::vec gradOld(p + 1), gradNew(p + 1);
1004 arma::mat betaProc(p + 1, m);
1005 for (int i = 0; i < m; i++) {
1006 double tau = tauSeq(i);
1007 arma::vec beta = huberReg(Z, Y, tau, der, gradOld, gradNew, n, p, n1, tol, constTau, iteMax);
1008 arma::vec quant = {tau};
1009 beta(0) = arma::as_scalar(arma::quantile(Y - Z.cols(1, p) * beta.rows(1, p), quant));
1010 arma::vec res = Y - Z * beta;
1011 updateUnif(Z, res, der, gradOld, n, tau, h, n1, h1);
1012 beta -= gradOld;
1013 arma::vec betaDiff = -gradOld;
1014 res -= Z * betaDiff;
1015 updateUnif(Z, res, der, gradOld, n, tau, h, n1, h1);
1016 arma::vec gradDiff = gradNew - gradOld;
1017 int ite = 1;
1018 while (arma::norm(gradNew, "inf") > tol && ite <= iteMax) {
1019 double alpha = 1.0;
1020 double cross = arma::as_scalar(betaDiff.t() * gradDiff);
1021 if (cross > 0) {
1022 double a1 = cross / arma::as_scalar(gradDiff.t() * gradDiff);
1023 double a2 = arma::as_scalar(betaDiff.t() * betaDiff) / cross;
1024 alpha = std::min(std::min(a1, a2), 100.0);
1025 }
1026 gradOld = gradNew;
1027 betaDiff = -alpha * gradNew;
1028 beta += betaDiff;
1029 res -= Z * betaDiff;
1030 updateUnif(Z, res, der, gradOld, n, tau, h, n1, h1);
1031 gradDiff = gradNew - gradOld;
1032 ite++;
1033 }
1034 betaProc.col(i) = beta;
1035 }
1036 betaProc.rows(1, p).each_col() %= sx1;
1037 betaProc.row(0) += my - mx * betaProc.rows(1, p);
1038 return Rcpp::List::create(Rcpp::Named("coeff") = betaProc, Rcpp::Named("bandwidth") = h);
1039 }
1040
1041 // [[Rcpp::export]]
smqrParaProc(const arma::mat & X,arma::vec Y,const arma::vec tauSeq,double h=0.0,const double constTau=1.345,const double tol=0.0001,const int iteMax=5000)1042 Rcpp::List smqrParaProc(const arma::mat& X, arma::vec Y, const arma::vec tauSeq, double h = 0.0, const double constTau = 1.345,
1043 const double tol = 0.0001, const int iteMax = 5000) {
1044 const int n = X.n_rows;
1045 const int p = X.n_cols;
1046 const int m = tauSeq.size();
1047 if (h <= 0.05) {
1048 h = std::max(std::pow((std::log(n) + p) / n, 0.4), 0.05);
1049 }
1050 const double n1 = 1.0 / n;
1051 const double h1 = 1.0 / h, h3 = 1.0 / (h * h * h);
1052 arma::rowvec mx = arma::mean(X, 0);
1053 arma::vec sx1 = 1.0 / arma::stddev(X, 0, 0).t();
1054 arma::mat Z = arma::join_rows(arma::ones(n), standardize(X, mx, sx1, p));
1055 double my = arma::mean(Y);
1056 Y -= my;
1057 arma::vec der(n);
1058 arma::vec gradOld(p + 1), gradNew(p + 1);
1059 arma::mat betaProc(p + 1, m);
1060 for (int i = 0; i < m; i++) {
1061 double tau = tauSeq(i);
1062 arma::vec beta = huberReg(Z, Y, tau, der, gradOld, gradNew, n, p, n1, tol, constTau, iteMax);
1063 arma::vec quant = {tau};
1064 beta(0) = arma::as_scalar(arma::quantile(Y - Z.cols(1, p) * beta.rows(1, p), quant));
1065 arma::vec res = Y - Z * beta;
1066 updatePara(Z, res, der, gradOld, n, tau, h, n1, h1, h3);
1067 beta -= gradOld;
1068 arma::vec betaDiff = -gradOld;
1069 res -= Z * betaDiff;
1070 updatePara(Z, res, der, gradOld, n, tau, h, n1, h1, h3);
1071 arma::vec gradDiff = gradNew - gradOld;
1072 int ite = 1;
1073 while (arma::norm(gradNew, "inf") > tol && ite <= iteMax) {
1074 double alpha = 1.0;
1075 double cross = arma::as_scalar(betaDiff.t() * gradDiff);
1076 if (cross > 0) {
1077 double a1 = cross / arma::as_scalar(gradDiff.t() * gradDiff);
1078 double a2 = arma::as_scalar(betaDiff.t() * betaDiff) / cross;
1079 alpha = std::min(std::min(a1, a2), 100.0);
1080 }
1081 gradOld = gradNew;
1082 betaDiff = -alpha * gradNew;
1083 beta += betaDiff;
1084 res -= Z * betaDiff;
1085 updatePara(Z, res, der, gradOld, n, tau, h, n1, h1, h3);
1086 gradDiff = gradNew - gradOld;
1087 ite++;
1088 }
1089 betaProc.col(i) = beta;
1090 }
1091 betaProc.rows(1, p).each_col() %= sx1;
1092 betaProc.row(0) += my - mx * betaProc.rows(1, p);
1093 return Rcpp::List::create(Rcpp::Named("coeff") = betaProc, Rcpp::Named("bandwidth") = h);
1094 }
1095
1096 // [[Rcpp::export]]
smqrTrianProc(const arma::mat & X,arma::vec Y,const arma::vec tauSeq,double h=0.0,const double constTau=1.345,const double tol=0.0001,const int iteMax=5000)1097 Rcpp::List smqrTrianProc(const arma::mat& X, arma::vec Y, const arma::vec tauSeq, double h = 0.0, const double constTau = 1.345,
1098 const double tol = 0.0001, const int iteMax = 5000) {
1099 const int n = X.n_rows;
1100 const int p = X.n_cols;
1101 const int m = tauSeq.size();
1102 if (h <= 0.05) {
1103 h = std::max(std::pow((std::log(n) + p) / n, 0.4), 0.05);
1104 }
1105 const double n1 = 1.0 / n;
1106 const double h1 = 1.0 / h, h2 = 1.0 / (h * h);
1107 arma::rowvec mx = arma::mean(X, 0);
1108 arma::vec sx1 = 1.0 / arma::stddev(X, 0, 0).t();
1109 arma::mat Z = arma::join_rows(arma::ones(n), standardize(X, mx, sx1, p));
1110 double my = arma::mean(Y);
1111 Y -= my;
1112 arma::vec der(n);
1113 arma::vec gradOld(p + 1), gradNew(p + 1);
1114 arma::mat betaProc(p + 1, m);
1115 for (int i = 0; i < m; i++) {
1116 double tau = tauSeq(i);
1117 arma::vec beta = huberReg(Z, Y, tau, der, gradOld, gradNew, n, p, n1, tol, constTau, iteMax);
1118 arma::vec quant = {tau};
1119 beta(0) = arma::as_scalar(arma::quantile(Y - Z.cols(1, p) * beta.rows(1, p), quant));
1120 arma::vec res = Y - Z * beta;
1121 updateTrian(Z, res, der, gradOld, n, tau, h, n1, h1, h2);
1122 beta -= gradOld;
1123 arma::vec betaDiff = -gradOld;
1124 res -= Z * betaDiff;
1125 updateTrian(Z, res, der, gradOld, n, tau, h, n1, h1, h2);
1126 arma::vec gradDiff = gradNew - gradOld;
1127 int ite = 1;
1128 while (arma::norm(gradNew, "inf") > tol && ite <= iteMax) {
1129 double alpha = 1.0;
1130 double cross = arma::as_scalar(betaDiff.t() * gradDiff);
1131 if (cross > 0) {
1132 double a1 = cross / arma::as_scalar(gradDiff.t() * gradDiff);
1133 double a2 = arma::as_scalar(betaDiff.t() * betaDiff) / cross;
1134 alpha = std::min(std::min(a1, a2), 100.0);
1135 }
1136 gradOld = gradNew;
1137 betaDiff = -alpha * gradNew;
1138 beta += betaDiff;
1139 res -= Z * betaDiff;
1140 updateTrian(Z, res, der, gradOld, n, tau, h, n1, h1, h2);
1141 gradDiff = gradNew - gradOld;
1142 ite++;
1143 }
1144 betaProc.col(i) = beta;
1145 }
1146 betaProc.rows(1, p).each_col() %= sx1;
1147 betaProc.row(0) += my - mx * betaProc.rows(1, p);
1148 return Rcpp::List::create(Rcpp::Named("coeff") = betaProc, Rcpp::Named("bandwidth") = h);
1149 }
1150
1151 // Conquer with bootstrap inference
1152 // [[Rcpp::export]]
smqrGaussInf(const arma::mat & X,const arma::vec & Y,const arma::vec & betaHat,const int n,const int p,double h=0.0,const double tau=0.5,const int B=1000,const double tol=0.0001,const int iteMax=5000)1153 arma::mat smqrGaussInf(const arma::mat& X, const arma::vec& Y, const arma::vec& betaHat, const int n, const int p, double h = 0.0, const double tau = 0.5,
1154 const int B = 1000, const double tol = 0.0001, const int iteMax = 5000) {
1155 arma::mat rst(p + 1, B);
1156 if (h <= 0.05) {
1157 h = std::max(std::pow((std::log(n) + p) / n, 0.4), 0.05);
1158 }
1159 for (int b = 0; b < B; b++) {
1160 arma::uvec idx = arma::find(arma::randi(n, arma::distr_param(0, 1)) == 1);
1161 arma::mat mbX = X.rows(idx);
1162 arma::mat mbY = Y.rows(idx);
1163 rst.col(b) = smqrGaussIni(mbX, mbY, betaHat, p, tau, h, tol, iteMax);
1164 }
1165 return rst;
1166 }
1167
1168 // [[Rcpp::export]]
smqrGaussInfWeight(const arma::mat & X,const arma::vec & Y,const arma::vec & betaHat,const int n,const int p,const double tau=0.5,const int B=1000,const double tol=0.0001,const int iteMax=5000)1169 arma::mat smqrGaussInfWeight(const arma::mat& X, const arma::vec& Y, const arma::vec& betaHat, const int n, const int p, const double tau = 0.5,
1170 const int B = 1000, const double tol = 0.0001, const int iteMax = 5000) {
1171 arma::mat rst(p + 1, B);
1172 double h = std::max(std::pow((std::log(n) + p) / n, 0.4), 0.05);
1173 for (int b = 0; b < B; b++) {
1174 arma::vec weight = 2 * arma::randu(n);
1175 rst.col(b) = smqrGaussIniWeight(X, Y, weight, betaHat, p, tau, h, tol, iteMax);
1176 }
1177 return rst;
1178 }
1179
1180 // [[Rcpp::export]]
smqrLogisticInf(const arma::mat & X,const arma::vec & Y,const arma::vec & betaHat,const int n,const int p,double h=0.0,const double tau=0.5,const int B=1000,const double tol=0.0001,const int iteMax=5000)1181 arma::mat smqrLogisticInf(const arma::mat& X, const arma::vec& Y, const arma::vec& betaHat, const int n, const int p, double h = 0.0, const double tau = 0.5,
1182 const int B = 1000, const double tol = 0.0001, const int iteMax = 5000) {
1183 arma::mat rst(p + 1, B);
1184 if (h <= 0.05) {
1185 h = std::max(std::pow((std::log(n) + p) / n, 0.4), 0.05);
1186 }
1187 for (int b = 0; b < B; b++) {
1188 arma::uvec idx = arma::find(arma::randi(n, arma::distr_param(0, 1)) == 1);
1189 arma::mat mbX = X.rows(idx);
1190 arma::mat mbY = Y.rows(idx);
1191 rst.col(b) = smqrLogisticIni(mbX, mbY, betaHat, p, tau, h, tol, iteMax);
1192 }
1193 return rst;
1194 }
1195
1196 // [[Rcpp::export]]
smqrUnifInf(const arma::mat & X,const arma::vec & Y,const arma::vec & betaHat,const int n,const int p,double h=0.0,const double tau=0.5,const int B=1000,const double tol=0.0001,const int iteMax=5000)1197 arma::mat smqrUnifInf(const arma::mat& X, const arma::vec& Y, const arma::vec& betaHat, const int n, const int p, double h = 0.0, const double tau = 0.5,
1198 const int B = 1000, const double tol = 0.0001, const int iteMax = 5000) {
1199 arma::mat rst(p + 1, B);
1200 if (h <= 0.05) {
1201 h = std::max(std::pow((std::log(n) + p) / n, 0.4), 0.05);
1202 }
1203 for (int b = 0; b < B; b++) {
1204 arma::uvec idx = arma::find(arma::randi(n, arma::distr_param(0, 1)) == 1);
1205 arma::mat mbX = X.rows(idx);
1206 arma::mat mbY = Y.rows(idx);
1207 rst.col(b) = smqrUnifIni(mbX, mbY, betaHat, p, tau, h, tol, iteMax);
1208 }
1209 return rst;
1210 }
1211
1212 // [[Rcpp::export]]
smqrParaInf(const arma::mat & X,const arma::vec & Y,const arma::vec & betaHat,const int n,const int p,double h=0.0,const double tau=0.5,const int B=1000,const double tol=0.0001,const int iteMax=5000)1213 arma::mat smqrParaInf(const arma::mat& X, const arma::vec& Y, const arma::vec& betaHat, const int n, const int p, double h = 0.0, const double tau = 0.5,
1214 const int B = 1000, const double tol = 0.0001, const int iteMax = 5000) {
1215 arma::mat rst(p + 1, B);
1216 if (h <= 0.05) {
1217 h = std::max(std::pow((std::log(n) + p) / n, 0.4), 0.05);
1218 }
1219 for (int b = 0; b < B; b++) {
1220 arma::uvec idx = arma::find(arma::randi(n, arma::distr_param(0, 1)) == 1);
1221 arma::mat mbX = X.rows(idx);
1222 arma::mat mbY = Y.rows(idx);
1223 rst.col(b) = smqrParaIni(mbX, mbY, betaHat, p, tau, h, tol, iteMax);
1224 }
1225 return rst;
1226 }
1227
1228 // [[Rcpp::export]]
smqrTrianInf(const arma::mat & X,const arma::vec & Y,const arma::vec & betaHat,const int n,const int p,double h=0.0,const double tau=0.5,const int B=1000,const double tol=0.0001,const int iteMax=5000)1229 arma::mat smqrTrianInf(const arma::mat& X, const arma::vec& Y, const arma::vec& betaHat, const int n, const int p, double h = 0.0, const double tau = 0.5,
1230 const int B = 1000, const double tol = 0.0001, const int iteMax = 5000) {
1231 arma::mat rst(p + 1, B);
1232 if (h <= 0.05) {
1233 h = std::max(std::pow((std::log(n) + p) / n, 0.4), 0.05);
1234 }
1235 for (int b = 0; b < B; b++) {
1236 arma::uvec idx = arma::find(arma::randi(n, arma::distr_param(0, 1)) == 1);
1237 arma::mat mbX = X.rows(idx);
1238 arma::mat mbY = Y.rows(idx);
1239 rst.col(b) = smqrTrianIni(mbX, mbY, betaHat, p, tau, h, tol, iteMax);
1240 }
1241 return rst;
1242 }
1243
1244 // The following code is high-dimensional conquer via an iterative local majorize-minimize algorithm
1245 // [[Rcpp::export]]
softThresh(const arma::vec & x,const arma::vec & lambda,const int p)1246 arma::vec softThresh(const arma::vec& x, const arma::vec& lambda, const int p) {
1247 return arma::sign(x) % arma::max(arma::abs(x) - lambda, arma::zeros(p + 1));
1248 }
1249
1250 // [[Rcpp::export]]
cmptLambdaLasso(const double lambda,const int p)1251 arma::vec cmptLambdaLasso(const double lambda, const int p) {
1252 arma::vec rst = lambda * arma::ones(p + 1);
1253 rst(0) = 0;
1254 return rst;
1255 }
1256
1257 // [[Rcpp::export]]
cmptLambdaSCAD(const arma::vec & beta,const double lambda,const int p,const double para=3.7)1258 arma::vec cmptLambdaSCAD(const arma::vec& beta, const double lambda, const int p, const double para = 3.7) {
1259 arma::vec rst = arma::zeros(p + 1);
1260 for (int i = 1; i <= p; i++) {
1261 double abBeta = std::abs(beta(i));
1262 if (abBeta <= lambda) {
1263 rst(i) = lambda;
1264 } else if (abBeta <= para * lambda) {
1265 rst(i) = (para * lambda - abBeta) / (para - 1);
1266 }
1267 }
1268 return rst;
1269 }
1270
1271 // [[Rcpp::export]]
cmptLambdaMCP(const arma::vec & beta,const double lambda,const int p,const double para=3.0)1272 arma::vec cmptLambdaMCP(const arma::vec& beta, const double lambda, const int p, const double para = 3.0) {
1273 arma::vec rst = arma::zeros(p + 1);
1274 for (int i = 1; i <= p; i++) {
1275 double abBeta = std::abs(beta(i));
1276 if (abBeta <= para * lambda) {
1277 rst(i) = lambda - abBeta / para;
1278 }
1279 }
1280 return rst;
1281 }
1282
1283 // Expectile regression (asymmetric l_2 loss) as an initial value for high-dim regression
1284 // [[Rcpp::export]]
lossL2(const arma::mat & Z,const arma::vec & Y,const arma::vec & beta,const double n1,const double tau)1285 double lossL2(const arma::mat& Z, const arma::vec& Y, const arma::vec& beta, const double n1, const double tau) {
1286 arma::vec res = Y - Z * beta;
1287 double rst = 0.0;
1288 for (int i = 0; i < Y.size(); i++) {
1289 rst += (res(i) > 0) ? (tau * res(i) * res(i)) : ((1 - tau) * res(i) * res(i));
1290 }
1291 return 0.5 * n1 * rst;
1292 }
1293
1294 // [[Rcpp::export]]
updateL2(const arma::mat & Z,const arma::vec & Y,const arma::vec & beta,arma::vec & grad,const double n1,const double tau)1295 double updateL2(const arma::mat& Z, const arma::vec& Y, const arma::vec& beta, arma::vec& grad, const double n1, const double tau) {
1296 arma::vec res = Y - Z * beta;
1297 double rst = 0.0;
1298 grad = arma::zeros(grad.size());
1299 for (int i = 0; i < Y.size(); i++) {
1300 double temp = res(i) > 0 ? tau : (1 - tau);
1301 grad -= temp * res(i) * Z.row(i).t();
1302 rst += temp * res(i) * res(i);
1303 }
1304 grad *= n1;
1305 return 0.5 * n1 * rst;
1306 }
1307
1308 // Smoothed quantile loss with different kernels
1309 // [[Rcpp::export]]
lossGaussHd(const arma::mat & Z,const arma::vec & Y,const arma::vec & beta,const double tau,const double h,const double h1,const double h2)1310 double lossGaussHd(const arma::mat& Z, const arma::vec& Y, const arma::vec& beta, const double tau, const double h, const double h1, const double h2) {
1311 arma::vec res = Y - Z * beta;
1312 arma::vec temp = 0.3989423 * h * arma::exp(-0.5 * h2 * arma::square(res)) + tau * res - res % arma::normcdf(-h1 * res);
1313 return arma::mean(temp);
1314 }
1315
1316 // [[Rcpp::export]]
updateGaussHd(const arma::mat & Z,const arma::vec & Y,const arma::vec & beta,arma::vec & grad,const double tau,const double n1,const double h,const double h1,const double h2)1317 double updateGaussHd(const arma::mat& Z, const arma::vec& Y, const arma::vec& beta, arma::vec& grad, const double tau, const double n1, const double h,
1318 const double h1, const double h2) {
1319 arma::vec res = Y - Z * beta;
1320 arma::vec der = arma::normcdf(-h1 * res) - tau;
1321 grad = n1 * Z.t() * der;
1322 arma::vec temp = 0.3989423 * h * arma::exp(-0.5 * h2 * arma::square(res)) + tau * res - res % arma::normcdf(-h1 * res);
1323 return arma::mean(temp);
1324 }
1325
1326 // [[Rcpp::export]]
lossLogisticHd(const arma::mat & Z,const arma::vec & Y,const arma::vec & beta,const double tau,const double h,const double h1)1327 double lossLogisticHd(const arma::mat& Z, const arma::vec& Y, const arma::vec& beta, const double tau, const double h, const double h1) {
1328 arma::vec res = Y - Z * beta;
1329 arma::vec temp = tau * res + h * arma::log(1.0 + arma::exp(-h1 * res));
1330 return arma::mean(temp);
1331 }
1332
1333 // [[Rcpp::export]]
updateLogisticHd(const arma::mat & Z,const arma::vec & Y,const arma::vec & beta,arma::vec & grad,const double tau,const double n1,const double h,const double h1)1334 double updateLogisticHd(const arma::mat& Z, const arma::vec& Y, const arma::vec& beta, arma::vec& grad, const double tau, const double n1, const double h,
1335 const double h1) {
1336 arma::vec res = Y - Z * beta;
1337 arma::vec der = 1.0 / (1.0 + arma::exp(res * h1)) - tau;
1338 grad = n1 * Z.t() * der;
1339 arma::vec temp = tau * res + h * arma::log(1.0 + arma::exp(-h1 * res));
1340 return arma::mean(temp);
1341 }
1342
1343 // [[Rcpp::export]]
lossUnifHd(const arma::mat & Z,const arma::vec & Y,const arma::vec & beta,const double tau,const double h,const double h1)1344 double lossUnifHd(const arma::mat& Z, const arma::vec& Y, const arma::vec& beta, const double tau, const double h, const double h1) {
1345 arma::vec res = Y - Z * beta;
1346 arma::vec temp = (tau - 0.5) * res;
1347 for (int i = 0; i < res.size(); i++) {
1348 double cur = std::abs(res(i));
1349 temp(i) += cur <= h ? (0.25 * h1 * cur * cur + 0.25 * h) : 0.5 * cur;
1350 }
1351 return arma::mean(temp);
1352 }
1353
1354 // [[Rcpp::export]]
updateUnifHd(const arma::mat & Z,const arma::vec & Y,const arma::vec & beta,arma::vec & grad,const double tau,const double n1,const double h,const double h1)1355 double updateUnifHd(const arma::mat& Z, const arma::vec& Y, const arma::vec& beta, arma::vec& grad, const double tau, const double n1, const double h,
1356 const double h1) {
1357 arma::vec res = Y - Z * beta;
1358 arma::vec temp = (tau - 0.5) * res;
1359 arma::vec der(res.size());
1360 for (int i = 0; i < res.size(); i++) {
1361 double cur = res(i);
1362 if (cur <= -h) {
1363 der(i) = 1 - tau;
1364 temp(i) -= 0.5 * cur;
1365 } else if (cur < h) {
1366 der(i) = 0.5 - tau - 0.5 * h1 * cur;
1367 temp(i) += 0.25 * h1 * cur * cur + 0.25 * h;
1368 } else {
1369 der(i) = -tau;
1370 temp(i) += 0.5 * cur;
1371 }
1372 }
1373 grad = n1 * Z.t() * der;
1374 return arma::mean(temp);
1375 }
1376
1377 // [[Rcpp::export]]
lossParaHd(const arma::mat & Z,const arma::vec & Y,const arma::vec & beta,const double tau,const double h,const double h1,const double h3)1378 double lossParaHd(const arma::mat& Z, const arma::vec& Y, const arma::vec& beta, const double tau, const double h, const double h1, const double h3) {
1379 arma::vec res = Y - Z * beta;
1380 arma::vec temp = (tau - 0.5) * res;
1381 for (int i = 0; i < res.size(); i++) {
1382 double cur = std::abs(res(i));
1383 temp(i) += cur <= h ? (0.375 * h1 * cur * cur - 0.0625 * h3 * cur * cur * cur * cur + 0.1875 * h) : 0.5 * cur;
1384 }
1385 return arma::mean(temp);
1386 }
1387
1388 // [[Rcpp::export]]
updateParaHd(const arma::mat & Z,const arma::vec & Y,const arma::vec & beta,arma::vec & grad,const double tau,const double n1,const double h,const double h1,const double h3)1389 double updateParaHd(const arma::mat& Z, const arma::vec& Y, const arma::vec& beta, arma::vec& grad, const double tau, const double n1, const double h,
1390 const double h1, const double h3) {
1391 arma::vec res = Y - Z * beta;
1392 arma::vec temp = (tau - 0.5) * res;
1393 arma::vec der(res.size());
1394 for (int i = 0; i < res.size(); i++) {
1395 double cur = res(i);
1396 if (cur <= -h) {
1397 der(i) = 1 - tau;
1398 temp(i) -= 0.5 * cur;
1399 } else if (cur < h) {
1400 der(i) = 0.5 - tau - 0.75 * h1 * cur + 0.25 * h3 * cur * cur * cur;
1401 temp(i) += 0.375 * h1 * cur * cur - 0.0625 * h3 * cur * cur * cur * cur + 0.1875 * h;
1402 } else {
1403 der(i) = -tau;
1404 temp(i) += 0.5 * cur;
1405 }
1406 }
1407 grad = n1 * Z.t() * der;
1408 return arma::mean(temp);
1409 }
1410
1411 // [[Rcpp::export]]
lossTrianHd(const arma::mat & Z,const arma::vec & Y,const arma::vec & beta,const double tau,const double h,const double h1,const double h2)1412 double lossTrianHd(const arma::mat& Z, const arma::vec& Y, const arma::vec& beta, const double tau, const double h, const double h1, const double h2) {
1413 arma::vec res = Y - Z * beta;
1414 arma::vec temp = (tau - 0.5) * res;
1415 for (int i = 0; i < res.size(); i++) {
1416 double cur = std::abs(res(i));
1417 temp(i) += cur <= h ? (0.5 * h1 * cur * cur - 0.1666667 * h2 * cur * cur * cur + 0.1666667 * h) : 0.5 * cur;
1418 }
1419 return arma::mean(temp);
1420 }
1421
1422 // [[Rcpp::export]]
updateTrianHd(const arma::mat & Z,const arma::vec & Y,const arma::vec & beta,arma::vec & grad,const double tau,const double n1,const double h,const double h1,const double h2)1423 double updateTrianHd(const arma::mat& Z, const arma::vec& Y, const arma::vec& beta, arma::vec& grad, const double tau, const double n1, const double h,
1424 const double h1, const double h2) {
1425 arma::vec res = Y - Z * beta;
1426 arma::vec temp = (tau - 0.5) * res;
1427 arma::vec der(res.size());
1428 for (int i = 0; i < res.size(); i++) {
1429 double cur = res(i);
1430 if (cur <= -h) {
1431 der(i) = 1 - tau;
1432 temp(i) -= 0.5 * cur;
1433 } else if (cur < 0) {
1434 der(i) = 0.5 - tau - h1 * cur - 0.5 * h2 * cur * cur;
1435 temp(i) += 0.5 * h1 * cur * cur + 0.1666667 * h2 * cur * cur * cur + 0.1666667 * h;
1436 } else if (cur < h) {
1437 der(i) = 0.5 - tau - h1 * cur + 0.5 * h2 * cur * cur;
1438 temp(i) += 0.5 * h1 * cur * cur - 0.1666667 * h2 * cur * cur * cur + 0.1666667 * h;
1439 } else {
1440 der(i) = -tau;
1441 temp(i) += 0.5 * cur;
1442 }
1443 }
1444 grad = n1 * Z.t() * der;
1445 return arma::mean(temp);
1446 }
1447
1448 // LAMM core code for different loss functions, update beta, return phi
1449 // [[Rcpp::export]]
lammL2(const arma::mat & Z,const arma::vec & Y,const arma::vec & Lambda,arma::vec & beta,const double tau,const double phi,const double gamma,const int p,const double n1)1450 double lammL2(const arma::mat& Z, const arma::vec& Y, const arma::vec& Lambda, arma::vec& beta, const double tau, const double phi, const double gamma,
1451 const int p, const double n1) {
1452 double phiNew = phi;
1453 arma::vec betaNew(p + 1);
1454 arma::vec grad(p + 1);
1455 double loss = updateL2(Z, Y, beta, grad, n1, tau);
1456 while (true) {
1457 arma::vec first = beta - grad / phiNew;
1458 arma::vec second = Lambda / phiNew;
1459 betaNew = softThresh(first, second, p);
1460 double fVal = lossL2(Z, Y, betaNew, n1, tau);
1461 arma::vec diff = betaNew - beta;
1462 double psiVal = loss + arma::as_scalar(grad.t() * diff) + 0.5 * phiNew * arma::as_scalar(diff.t() * diff);
1463 if (fVal <= psiVal) {
1464 break;
1465 }
1466 phiNew *= gamma;
1467 }
1468 beta = betaNew;
1469 return phiNew;
1470 }
1471
1472 // [[Rcpp::export]]
lasso(const arma::mat & Z,const arma::vec & Y,const double lambda,const double tau,const int p,const double n1,const double phi0=0.1,const double gamma=1.2,const double epsilon=0.01,const int iteMax=500)1473 arma::vec lasso(const arma::mat& Z, const arma::vec& Y, const double lambda, const double tau, const int p, const double n1, const double phi0 = 0.1,
1474 const double gamma = 1.2, const double epsilon = 0.01, const int iteMax = 500) {
1475 arma::vec beta = arma::zeros(p + 1);
1476 arma::vec betaNew = arma::zeros(p + 1);
1477 arma::vec Lambda = cmptLambdaLasso(lambda, p);
1478 double phi = phi0;
1479 int ite = 0;
1480 while (ite <= iteMax) {
1481 ite++;
1482 phi = lammL2(Z, Y, Lambda, betaNew, tau, phi, gamma, p, n1);
1483 phi = std::max(phi0, phi / gamma);
1484 if (arma::norm(betaNew - beta, "inf") <= epsilon) {
1485 break;
1486 }
1487 beta = betaNew;
1488 }
1489 return betaNew;
1490 }
1491
1492 // [[Rcpp::export]]
lammSmqrGauss(const arma::mat & Z,const arma::vec & Y,const arma::vec & Lambda,arma::vec & beta,const double phi,const double tau,const double gamma,const int p,const double h,const double n1,const double h1,const double h2)1493 double lammSmqrGauss(const arma::mat& Z, const arma::vec& Y, const arma::vec& Lambda, arma::vec& beta, const double phi, const double tau,
1494 const double gamma, const int p, const double h, const double n1, const double h1, const double h2) {
1495 double phiNew = phi;
1496 arma::vec betaNew(p + 1);
1497 arma::vec grad(p + 1);
1498 double loss = updateGaussHd(Z, Y, beta, grad, tau, n1, h, h1, h2);
1499 while (true) {
1500 arma::vec first = beta - grad / phiNew;
1501 arma::vec second = Lambda / phiNew;
1502 betaNew = softThresh(first, second, p);
1503 double fVal = lossGaussHd(Z, Y, betaNew, tau, h, h1, h2);
1504 arma::vec diff = betaNew - beta;
1505 double psiVal = loss + arma::as_scalar(grad.t() * diff) + 0.5 * phiNew * arma::as_scalar(diff.t() * diff);
1506 if (fVal <= psiVal) {
1507 break;
1508 }
1509 phiNew *= gamma;
1510 }
1511 beta = betaNew;
1512 return phiNew;
1513 }
1514
1515 // [[Rcpp::export]]
lammSmqrLogistic(const arma::mat & Z,const arma::vec & Y,const arma::vec & Lambda,arma::vec & beta,const double phi,const double tau,const double gamma,const int p,const double h,const double n1,const double h1)1516 double lammSmqrLogistic(const arma::mat& Z, const arma::vec& Y, const arma::vec& Lambda, arma::vec& beta, const double phi, const double tau,
1517 const double gamma, const int p, const double h, const double n1, const double h1) {
1518 double phiNew = phi;
1519 arma::vec betaNew(p + 1);
1520 arma::vec grad(p + 1);
1521 double loss = updateLogisticHd(Z, Y, beta, grad, tau, n1, h, h1);
1522 while (true) {
1523 arma::vec first = beta - grad / phiNew;
1524 arma::vec second = Lambda / phiNew;
1525 betaNew = softThresh(first, second, p);
1526 double fVal = lossLogisticHd(Z, Y, betaNew, tau, h, h1);
1527 arma::vec diff = betaNew - beta;
1528 double psiVal = loss + arma::as_scalar(grad.t() * diff) + 0.5 * phiNew * arma::as_scalar(diff.t() * diff);
1529 if (fVal <= psiVal) {
1530 break;
1531 }
1532 phiNew *= gamma;
1533 }
1534 beta = betaNew;
1535 return phiNew;
1536 }
1537
1538 // [[Rcpp::export]]
lammSmqrUnif(const arma::mat & Z,const arma::vec & Y,const arma::vec & Lambda,arma::vec & beta,const double phi,const double tau,const double gamma,const int p,const double h,const double n1,const double h1)1539 double lammSmqrUnif(const arma::mat& Z, const arma::vec& Y, const arma::vec& Lambda, arma::vec& beta, const double phi, const double tau,
1540 const double gamma, const int p, const double h, const double n1, const double h1) {
1541 double phiNew = phi;
1542 arma::vec betaNew(p + 1);
1543 arma::vec grad(p + 1);
1544 double loss = updateUnifHd(Z, Y, beta, grad, tau, n1, h, h1);
1545 while (true) {
1546 arma::vec first = beta - grad / phiNew;
1547 arma::vec second = Lambda / phiNew;
1548 betaNew = softThresh(first, second, p);
1549 double fVal = lossUnifHd(Z, Y, betaNew, tau, h, h1);
1550 arma::vec diff = betaNew - beta;
1551 double psiVal = loss + arma::as_scalar(grad.t() * diff) + 0.5 * phiNew * arma::as_scalar(diff.t() * diff);
1552 if (fVal <= psiVal) {
1553 break;
1554 }
1555 phiNew *= gamma;
1556 }
1557 beta = betaNew;
1558 return phiNew;
1559 }
1560
1561 // [[Rcpp::export]]
lammSmqrPara(const arma::mat & Z,const arma::vec & Y,const arma::vec & Lambda,arma::vec & beta,const double phi,const double tau,const double gamma,const int p,const double h,const double n1,const double h1,const double h3)1562 double lammSmqrPara(const arma::mat& Z, const arma::vec& Y, const arma::vec& Lambda, arma::vec& beta, const double phi, const double tau,
1563 const double gamma, const int p, const double h, const double n1, const double h1, const double h3) {
1564 double phiNew = phi;
1565 arma::vec betaNew(p + 1);
1566 arma::vec grad(p + 1);
1567 double loss = updateParaHd(Z, Y, beta, grad, tau, n1, h, h1, h3);
1568 while (true) {
1569 arma::vec first = beta - grad / phiNew;
1570 arma::vec second = Lambda / phiNew;
1571 betaNew = softThresh(first, second, p);
1572 double fVal = lossParaHd(Z, Y, betaNew, tau, h, h1, h3);
1573 arma::vec diff = betaNew - beta;
1574 double psiVal = loss + arma::as_scalar(grad.t() * diff) + 0.5 * phiNew * arma::as_scalar(diff.t() * diff);
1575 if (fVal <= psiVal) {
1576 break;
1577 }
1578 phiNew *= gamma;
1579 }
1580 beta = betaNew;
1581 return phiNew;
1582 }
1583
1584 // [[Rcpp::export]]
lammSmqrTrian(const arma::mat & Z,const arma::vec & Y,const arma::vec & Lambda,arma::vec & beta,const double phi,const double tau,const double gamma,const int p,const double h,const double n1,const double h1,const double h2)1585 double lammSmqrTrian(const arma::mat& Z, const arma::vec& Y, const arma::vec& Lambda, arma::vec& beta, const double phi, const double tau,
1586 const double gamma, const int p, const double h, const double n1, const double h1, const double h2) {
1587 double phiNew = phi;
1588 arma::vec betaNew(p + 1);
1589 arma::vec grad(p + 1);
1590 double loss = updateTrianHd(Z, Y, beta, grad, tau, n1, h, h1, h2);
1591 while (true) {
1592 arma::vec first = beta - grad / phiNew;
1593 arma::vec second = Lambda / phiNew;
1594 betaNew = softThresh(first, second, p);
1595 double fVal = lossTrianHd(Z, Y, betaNew, tau, h, h1, h2);
1596 arma::vec diff = betaNew - beta;
1597 double psiVal = loss + arma::as_scalar(grad.t() * diff) + 0.5 * phiNew * arma::as_scalar(diff.t() * diff);
1598 if (fVal <= psiVal) {
1599 break;
1600 }
1601 phiNew *= gamma;
1602 }
1603 beta = betaNew;
1604 return phiNew;
1605 }
1606
1607 // High-dim conquer with a standardized design matrix and a given lambda
1608 // [[Rcpp::export]]
smqrLassoGauss(const arma::mat & Z,const arma::vec & Y,const double lambda,const arma::vec & sx1,const double tau,const int p,const double n1,const double h,const double h1,const double h2,const double phi0=0.01,const double gamma=1.2,const double epsilon=0.001,const int iteMax=500)1609 arma::vec smqrLassoGauss(const arma::mat& Z, const arma::vec& Y, const double lambda, const arma::vec& sx1, const double tau, const int p, const double n1,
1610 const double h, const double h1, const double h2, const double phi0 = 0.01, const double gamma = 1.2, const double epsilon = 0.001,
1611 const int iteMax = 500) {
1612 arma::vec beta = lasso(Z, Y, lambda, tau, p, n1, phi0, gamma, epsilon, iteMax);
1613 arma::vec betaNew = beta;
1614 arma::vec Lambda = cmptLambdaLasso(lambda, p);
1615 double phi = phi0;
1616 int ite = 0;
1617 while (ite <= iteMax) {
1618 ite++;
1619 phi = lammSmqrGauss(Z, Y, Lambda, betaNew, phi, tau, gamma, p, h, n1, h1, h2);
1620 phi = std::max(phi0, phi / gamma);
1621 if (arma::norm(betaNew - beta, "inf") <= epsilon) {
1622 break;
1623 }
1624 beta = betaNew;
1625 }
1626 return betaNew;
1627 }
1628
1629 // [[Rcpp::export]]
smqrLassoLogistic(const arma::mat & Z,const arma::vec & Y,const double lambda,const arma::vec & sx1,const double tau,const int p,const double n1,const double h,const double h1,const double phi0=0.01,const double gamma=1.2,const double epsilon=0.001,const int iteMax=500)1630 arma::vec smqrLassoLogistic(const arma::mat& Z, const arma::vec& Y, const double lambda, const arma::vec& sx1, const double tau, const int p, const double n1,
1631 const double h, const double h1, const double phi0 = 0.01, const double gamma = 1.2, const double epsilon = 0.001,
1632 const int iteMax = 500) {
1633 arma::vec beta = lasso(Z, Y, lambda, tau, p, n1, phi0, gamma, epsilon, iteMax);
1634 arma::vec betaNew = beta;
1635 arma::vec Lambda = cmptLambdaLasso(lambda, p);
1636 double phi = phi0;
1637 int ite = 0;
1638 while (ite <= iteMax) {
1639 ite++;
1640 phi = lammSmqrLogistic(Z, Y, Lambda, betaNew, phi, tau, gamma, p, h, n1, h1);
1641 phi = std::max(phi0, phi / gamma);
1642 if (arma::norm(betaNew - beta, "inf") <= epsilon) {
1643 break;
1644 }
1645 beta = betaNew;
1646 }
1647 return betaNew;
1648 }
1649
1650 // [[Rcpp::export]]
smqrLassoUnif(const arma::mat & Z,const arma::vec & Y,const double lambda,const arma::vec & sx1,const double tau,const int p,const double n1,const double h,const double h1,const double phi0=0.01,const double gamma=1.2,const double epsilon=0.001,const int iteMax=500)1651 arma::vec smqrLassoUnif(const arma::mat& Z, const arma::vec& Y, const double lambda, const arma::vec& sx1, const double tau, const int p, const double n1,
1652 const double h, const double h1, const double phi0 = 0.01, const double gamma = 1.2, const double epsilon = 0.001,
1653 const int iteMax = 500) {
1654 arma::vec beta = lasso(Z, Y, lambda, tau, p, n1, phi0, gamma, epsilon, iteMax);
1655 arma::vec betaNew = beta;
1656 arma::vec Lambda = cmptLambdaLasso(lambda, p);
1657 double phi = phi0;
1658 int ite = 0;
1659 while (ite <= iteMax) {
1660 ite++;
1661 phi = lammSmqrUnif(Z, Y, Lambda, betaNew, phi, tau, gamma, p, h, n1, h1);
1662 phi = std::max(phi0, phi / gamma);
1663 if (arma::norm(betaNew - beta, "inf") <= epsilon) {
1664 break;
1665 }
1666 beta = betaNew;
1667 }
1668 return betaNew;
1669 }
1670
1671 // [[Rcpp::export]]
smqrLassoPara(const arma::mat & Z,const arma::vec & Y,const double lambda,const arma::vec & sx1,const double tau,const int p,const double n1,const double h,const double h1,const double h3,const double phi0=0.01,const double gamma=1.2,const double epsilon=0.001,const int iteMax=500)1672 arma::vec smqrLassoPara(const arma::mat& Z, const arma::vec& Y, const double lambda, const arma::vec& sx1, const double tau, const int p, const double n1,
1673 const double h, const double h1, const double h3, const double phi0 = 0.01, const double gamma = 1.2, const double epsilon = 0.001,
1674 const int iteMax = 500) {
1675 arma::vec beta = lasso(Z, Y, lambda, tau, p, n1, phi0, gamma, epsilon, iteMax);
1676 arma::vec betaNew = beta;
1677 arma::vec Lambda = cmptLambdaLasso(lambda, p);
1678 double phi = phi0;
1679 int ite = 0;
1680 while (ite <= iteMax) {
1681 ite++;
1682 phi = lammSmqrPara(Z, Y, Lambda, betaNew, phi, tau, gamma, p, h, n1, h1, h3);
1683 phi = std::max(phi0, phi / gamma);
1684 if (arma::norm(betaNew - beta, "inf") <= epsilon) {
1685 break;
1686 }
1687 beta = betaNew;
1688 }
1689 return betaNew;
1690 }
1691
1692 // [[Rcpp::export]]
smqrLassoTrian(const arma::mat & Z,const arma::vec & Y,const double lambda,const arma::vec & sx1,const double tau,const int p,const double n1,const double h,const double h1,const double h2,const double phi0=0.01,const double gamma=1.2,const double epsilon=0.001,const int iteMax=500)1693 arma::vec smqrLassoTrian(const arma::mat& Z, const arma::vec& Y, const double lambda, const arma::vec& sx1, const double tau, const int p, const double n1,
1694 const double h, const double h1, const double h2, const double phi0 = 0.01, const double gamma = 1.2, const double epsilon = 0.001,
1695 const int iteMax = 500) {
1696 arma::vec beta = lasso(Z, Y, lambda, tau, p, n1, phi0, gamma, epsilon, iteMax);
1697 arma::vec betaNew = beta;
1698 arma::vec Lambda = cmptLambdaLasso(lambda, p);
1699 double phi = phi0;
1700 int ite = 0;
1701 while (ite <= iteMax) {
1702 ite++;
1703 phi = lammSmqrTrian(Z, Y, Lambda, betaNew, phi, tau, gamma, p, h, n1, h1, h2);
1704 phi = std::max(phi0, phi / gamma);
1705 if (arma::norm(betaNew - beta, "inf") <= epsilon) {
1706 break;
1707 }
1708 beta = betaNew;
1709 }
1710 return betaNew;
1711 }
1712
1713 // [[Rcpp::export]]
smqrScadGauss(const arma::mat & Z,const arma::vec & Y,const double lambda,const arma::vec & sx1,const double tau,const int p,const double n1,const double h,const double h1,const double h2,const double phi0=0.01,const double gamma=1.2,const double epsilon=0.001,const int iteMax=500,const int iteTight=3,const double para=3.7)1714 arma::vec smqrScadGauss(const arma::mat& Z, const arma::vec& Y, const double lambda, const arma::vec& sx1, const double tau, const int p, const double n1,
1715 const double h, const double h1, const double h2, const double phi0 = 0.01, const double gamma = 1.2, const double epsilon = 0.001,
1716 const int iteMax = 500, const int iteTight = 3, const double para = 3.7) {
1717 arma::vec beta = lasso(Z, Y, lambda, tau, p, n1, phi0, gamma, epsilon, iteMax);
1718 arma::vec betaNew = beta;
1719 // Contraction
1720 arma::vec Lambda = cmptLambdaSCAD(beta, lambda, p, para);
1721 double phi = phi0;
1722 int ite = 0;
1723 while (ite <= iteMax) {
1724 ite++;
1725 phi = lammSmqrGauss(Z, Y, Lambda, betaNew, phi, tau, gamma, p, h, n1, h1, h2);
1726 phi = std::max(phi0, phi / gamma);
1727 if (arma::norm(betaNew - beta, "inf") <= epsilon) {
1728 break;
1729 }
1730 beta = betaNew;
1731 }
1732 int iteT = 1;
1733 // Tightening
1734 arma::vec beta0(p + 1);
1735 while (iteT <= iteTight) {
1736 iteT++;
1737 beta = betaNew;
1738 beta0 = betaNew;
1739 Lambda = cmptLambdaSCAD(beta, lambda, p, para);
1740 phi = phi0;
1741 ite = 0;
1742 while (ite <= iteMax) {
1743 ite++;
1744 phi = lammSmqrGauss(Z, Y, Lambda, betaNew, phi, tau, gamma, p, h, n1, h1, h2);
1745 phi = std::max(phi0, phi / gamma);
1746 if (arma::norm(betaNew - beta, "inf") <= epsilon) {
1747 break;
1748 }
1749 beta = betaNew;
1750 }
1751 if (arma::norm(betaNew - beta0, "inf") <= epsilon) {
1752 break;
1753 }
1754 }
1755 return betaNew;
1756 }
1757
1758 // [[Rcpp::export]]
smqrScadLogistic(const arma::mat & Z,const arma::vec & Y,const double lambda,const arma::vec & sx1,const double tau,const int p,const double n1,const double h,const double h1,const double phi0=0.01,const double gamma=1.2,const double epsilon=0.001,const int iteMax=500,const int iteTight=3,const double para=3.7)1759 arma::vec smqrScadLogistic(const arma::mat& Z, const arma::vec& Y, const double lambda, const arma::vec& sx1, const double tau, const int p, const double n1,
1760 const double h, const double h1, const double phi0 = 0.01, const double gamma = 1.2, const double epsilon = 0.001,
1761 const int iteMax = 500, const int iteTight = 3, const double para = 3.7) {
1762 arma::vec beta = lasso(Z, Y, lambda, tau, p, n1, phi0, gamma, epsilon, iteMax);
1763 arma::vec betaNew = beta;
1764 // Contraction
1765 arma::vec Lambda = cmptLambdaSCAD(beta, lambda, p, para);
1766 double phi = phi0;
1767 int ite = 0;
1768 while (ite <= iteMax) {
1769 ite++;
1770 phi = lammSmqrLogistic(Z, Y, Lambda, betaNew, phi, tau, gamma, p, h, n1, h1);
1771 phi = std::max(phi0, phi / gamma);
1772 if (arma::norm(betaNew - beta, "inf") <= epsilon) {
1773 break;
1774 }
1775 beta = betaNew;
1776 }
1777 int iteT = 1;
1778 // Tightening
1779 arma::vec beta0(p + 1);
1780 while (iteT <= iteTight) {
1781 iteT++;
1782 beta = betaNew;
1783 beta0 = betaNew;
1784 Lambda = cmptLambdaSCAD(beta, lambda, p, para);
1785 phi = phi0;
1786 ite = 0;
1787 while (ite <= iteMax) {
1788 ite++;
1789 phi = lammSmqrLogistic(Z, Y, Lambda, betaNew, phi, tau, gamma, p, h, n1, h1);
1790 phi = std::max(phi0, phi / gamma);
1791 if (arma::norm(betaNew - beta, "inf") <= epsilon) {
1792 break;
1793 }
1794 beta = betaNew;
1795 }
1796 if (arma::norm(betaNew - beta0, "inf") <= epsilon) {
1797 break;
1798 }
1799 }
1800 return betaNew;
1801 }
1802
1803 // [[Rcpp::export]]
smqrScadUnif(const arma::mat & Z,const arma::vec & Y,const double lambda,const arma::vec & sx1,const double tau,const int p,const double n1,const double h,const double h1,const double phi0=0.01,const double gamma=1.2,const double epsilon=0.001,const int iteMax=500,const int iteTight=3,const double para=3.7)1804 arma::vec smqrScadUnif(const arma::mat& Z, const arma::vec& Y, const double lambda, const arma::vec& sx1, const double tau, const int p, const double n1,
1805 const double h, const double h1, const double phi0 = 0.01, const double gamma = 1.2, const double epsilon = 0.001,
1806 const int iteMax = 500, const int iteTight = 3, const double para = 3.7) {
1807 arma::vec beta = lasso(Z, Y, lambda, tau, p, n1, phi0, gamma, epsilon, iteMax);
1808 arma::vec betaNew = beta;
1809 // Contraction
1810 arma::vec Lambda = cmptLambdaSCAD(beta, lambda, p, para);
1811 double phi = phi0;
1812 int ite = 0;
1813 while (ite <= iteMax) {
1814 ite++;
1815 phi = lammSmqrUnif(Z, Y, Lambda, betaNew, phi, tau, gamma, p, h, n1, h1);
1816 phi = std::max(phi0, phi / gamma);
1817 if (arma::norm(betaNew - beta, "inf") <= epsilon) {
1818 break;
1819 }
1820 beta = betaNew;
1821 }
1822 int iteT = 1;
1823 // Tightening
1824 arma::vec beta0(p + 1);
1825 while (iteT <= iteTight) {
1826 iteT++;
1827 beta = betaNew;
1828 beta0 = betaNew;
1829 Lambda = cmptLambdaSCAD(beta, lambda, p, para);
1830 phi = phi0;
1831 ite = 0;
1832 while (ite <= iteMax) {
1833 ite++;
1834 phi = lammSmqrUnif(Z, Y, Lambda, betaNew, phi, tau, gamma, p, h, n1, h1);
1835 phi = std::max(phi0, phi / gamma);
1836 if (arma::norm(betaNew - beta, "inf") <= epsilon) {
1837 break;
1838 }
1839 beta = betaNew;
1840 }
1841 if (arma::norm(betaNew - beta0, "inf") <= epsilon) {
1842 break;
1843 }
1844 }
1845 return betaNew;
1846 }
1847
1848 // [[Rcpp::export]]
smqrScadPara(const arma::mat & Z,const arma::vec & Y,const double lambda,const arma::vec & sx1,const double tau,const int p,const double n1,const double h,const double h1,const double h3,const double phi0=0.01,const double gamma=1.2,const double epsilon=0.001,const int iteMax=500,const int iteTight=3,const double para=3.7)1849 arma::vec smqrScadPara(const arma::mat& Z, const arma::vec& Y, const double lambda, const arma::vec& sx1, const double tau, const int p, const double n1,
1850 const double h, const double h1, const double h3, const double phi0 = 0.01, const double gamma = 1.2, const double epsilon = 0.001,
1851 const int iteMax = 500, const int iteTight = 3, const double para = 3.7) {
1852 arma::vec beta = lasso(Z, Y, lambda, tau, p, n1, phi0, gamma, epsilon, iteMax);
1853 arma::vec betaNew = beta;
1854 // Contraction
1855 arma::vec Lambda = cmptLambdaSCAD(beta, lambda, p, para);
1856 double phi = phi0;
1857 int ite = 0;
1858 while (ite <= iteMax) {
1859 ite++;
1860 phi = lammSmqrPara(Z, Y, Lambda, betaNew, phi, tau, gamma, p, h, n1, h1, h3);
1861 phi = std::max(phi0, phi / gamma);
1862 if (arma::norm(betaNew - beta, "inf") <= epsilon) {
1863 break;
1864 }
1865 beta = betaNew;
1866 }
1867 int iteT = 1;
1868 // Tightening
1869 arma::vec beta0(p + 1);
1870 while (iteT <= iteTight) {
1871 iteT++;
1872 beta = betaNew;
1873 beta0 = betaNew;
1874 Lambda = cmptLambdaSCAD(beta, lambda, p, para);
1875 phi = phi0;
1876 ite = 0;
1877 while (ite <= iteMax) {
1878 ite++;
1879 phi = lammSmqrPara(Z, Y, Lambda, betaNew, phi, tau, gamma, p, h, n1, h1, h3);
1880 phi = std::max(phi0, phi / gamma);
1881 if (arma::norm(betaNew - beta, "inf") <= epsilon) {
1882 break;
1883 }
1884 beta = betaNew;
1885 }
1886 if (arma::norm(betaNew - beta0, "inf") <= epsilon) {
1887 break;
1888 }
1889 }
1890 return betaNew;
1891 }
1892
1893 // [[Rcpp::export]]
smqrScadTrian(const arma::mat & Z,const arma::vec & Y,const double lambda,const arma::vec & sx1,const double tau,const int p,const double n1,const double h,const double h1,const double h2,const double phi0=0.01,const double gamma=1.2,const double epsilon=0.001,const int iteMax=500,const int iteTight=3,const double para=3.7)1894 arma::vec smqrScadTrian(const arma::mat& Z, const arma::vec& Y, const double lambda, const arma::vec& sx1, const double tau, const int p, const double n1,
1895 const double h, const double h1, const double h2, const double phi0 = 0.01, const double gamma = 1.2, const double epsilon = 0.001,
1896 const int iteMax = 500, const int iteTight = 3, const double para = 3.7) {
1897 arma::vec beta = lasso(Z, Y, lambda, tau, p, n1, phi0, gamma, epsilon, iteMax);
1898 arma::vec betaNew = beta;
1899 // Contraction
1900 arma::vec Lambda = cmptLambdaSCAD(beta, lambda, p, para);
1901 double phi = phi0;
1902 int ite = 0;
1903 while (ite <= iteMax) {
1904 ite++;
1905 phi = lammSmqrTrian(Z, Y, Lambda, betaNew, phi, tau, gamma, p, h, n1, h1, h2);
1906 phi = std::max(phi0, phi / gamma);
1907 if (arma::norm(betaNew - beta, "inf") <= epsilon) {
1908 break;
1909 }
1910 beta = betaNew;
1911 }
1912 int iteT = 1;
1913 // Tightening
1914 arma::vec beta0(p + 1);
1915 while (iteT <= iteTight) {
1916 iteT++;
1917 beta = betaNew;
1918 beta0 = betaNew;
1919 Lambda = cmptLambdaSCAD(beta, lambda, p, para);
1920 phi = phi0;
1921 ite = 0;
1922 while (ite <= iteMax) {
1923 ite++;
1924 phi = lammSmqrTrian(Z, Y, Lambda, betaNew, phi, tau, gamma, p, h, n1, h1, h2);
1925 phi = std::max(phi0, phi / gamma);
1926 if (arma::norm(betaNew - beta, "inf") <= epsilon) {
1927 break;
1928 }
1929 beta = betaNew;
1930 }
1931 if (arma::norm(betaNew - beta0, "inf") <= epsilon) {
1932 break;
1933 }
1934 }
1935 return betaNew;
1936 }
1937
1938 // [[Rcpp::export]]
smqrMcpGauss(const arma::mat & Z,const arma::vec & Y,const double lambda,const arma::vec & sx1,const double tau,const int p,const double n1,const double h,const double h1,const double h2,const double phi0=0.01,const double gamma=1.2,const double epsilon=0.001,const int iteMax=500,const int iteTight=3,const double para=3)1939 arma::vec smqrMcpGauss(const arma::mat& Z, const arma::vec& Y, const double lambda, const arma::vec& sx1, const double tau, const int p, const double n1,
1940 const double h, const double h1, const double h2, const double phi0 = 0.01, const double gamma = 1.2, const double epsilon = 0.001,
1941 const int iteMax = 500, const int iteTight = 3, const double para = 3) {
1942 arma::vec beta = lasso(Z, Y, lambda, tau, p, n1, phi0, gamma, epsilon, iteMax);
1943 arma::vec betaNew = beta;
1944 // Contraction
1945 arma::vec Lambda = cmptLambdaMCP(beta, lambda, p, para);
1946 double phi = phi0;
1947 int ite = 0;
1948 while (ite <= iteMax) {
1949 ite++;
1950 phi = lammSmqrGauss(Z, Y, Lambda, betaNew, phi, tau, gamma, p, h, n1, h1, h2);
1951 phi = std::max(phi0, phi / gamma);
1952 if (arma::norm(betaNew - beta, "inf") <= epsilon) {
1953 break;
1954 }
1955 beta = betaNew;
1956 }
1957 int iteT = 1;
1958 // Tightening
1959 arma::vec beta0(p + 1);
1960 while (iteT <= iteTight) {
1961 iteT++;
1962 beta = betaNew;
1963 beta0 = betaNew;
1964 Lambda = cmptLambdaMCP(beta, lambda, p, para);
1965 phi = phi0;
1966 ite = 0;
1967 while (ite <= iteMax) {
1968 ite++;
1969 phi = lammSmqrGauss(Z, Y, Lambda, betaNew, phi, tau, gamma, p, h, n1, h1, h2);
1970 phi = std::max(phi0, phi / gamma);
1971 if (arma::norm(betaNew - beta, "inf") <= epsilon) {
1972 break;
1973 }
1974 beta = betaNew;
1975 }
1976 if (arma::norm(betaNew - beta0, "inf") <= epsilon) {
1977 break;
1978 }
1979 }
1980 return betaNew;
1981 }
1982
1983 // [[Rcpp::export]]
smqrMcpLogistic(const arma::mat & Z,const arma::vec & Y,const double lambda,const arma::vec & sx1,const double tau,const int p,const double n1,const double h,const double h1,const double phi0=0.01,const double gamma=1.2,const double epsilon=0.001,const int iteMax=500,const int iteTight=3,const double para=3)1984 arma::vec smqrMcpLogistic(const arma::mat& Z, const arma::vec& Y, const double lambda, const arma::vec& sx1, const double tau, const int p, const double n1,
1985 const double h, const double h1, const double phi0 = 0.01, const double gamma = 1.2, const double epsilon = 0.001,
1986 const int iteMax = 500, const int iteTight = 3, const double para = 3) {
1987 arma::vec beta = lasso(Z, Y, lambda, tau, p, n1, phi0, gamma, epsilon, iteMax);
1988 arma::vec betaNew = beta;
1989 // Contraction
1990 arma::vec Lambda = cmptLambdaMCP(beta, lambda, p, para);
1991 double phi = phi0;
1992 int ite = 0;
1993 while (ite <= iteMax) {
1994 ite++;
1995 phi = lammSmqrLogistic(Z, Y, Lambda, betaNew, phi, tau, gamma, p, h, n1, h1);
1996 phi = std::max(phi0, phi / gamma);
1997 if (arma::norm(betaNew - beta, "inf") <= epsilon) {
1998 break;
1999 }
2000 beta = betaNew;
2001 }
2002 int iteT = 1;
2003 // Tightening
2004 arma::vec beta0(p + 1);
2005 while (iteT <= iteTight) {
2006 iteT++;
2007 beta = betaNew;
2008 beta0 = betaNew;
2009 Lambda = cmptLambdaMCP(beta, lambda, p, para);
2010 phi = phi0;
2011 ite = 0;
2012 while (ite <= iteMax) {
2013 ite++;
2014 phi = lammSmqrLogistic(Z, Y, Lambda, betaNew, phi, tau, gamma, p, h, n1, h1);
2015 phi = std::max(phi0, phi / gamma);
2016 if (arma::norm(betaNew - beta, "inf") <= epsilon) {
2017 break;
2018 }
2019 beta = betaNew;
2020 }
2021 if (arma::norm(betaNew - beta0, "inf") <= epsilon) {
2022 break;
2023 }
2024 }
2025 return betaNew;
2026 }
2027
2028 // [[Rcpp::export]]
smqrMcpUnif(const arma::mat & Z,const arma::vec & Y,const double lambda,const arma::vec & sx1,const double tau,const int p,const double n1,const double h,const double h1,const double phi0=0.01,const double gamma=1.2,const double epsilon=0.001,const int iteMax=500,const int iteTight=3,const double para=3)2029 arma::vec smqrMcpUnif(const arma::mat& Z, const arma::vec& Y, const double lambda, const arma::vec& sx1, const double tau, const int p, const double n1,
2030 const double h, const double h1, const double phi0 = 0.01, const double gamma = 1.2, const double epsilon = 0.001,
2031 const int iteMax = 500, const int iteTight = 3, const double para = 3) {
2032 arma::vec beta = lasso(Z, Y, lambda, tau, p, n1, phi0, gamma, epsilon, iteMax);
2033 arma::vec betaNew = beta;
2034 // Contraction
2035 arma::vec Lambda = cmptLambdaMCP(beta, lambda, p, para);
2036 double phi = phi0;
2037 int ite = 0;
2038 while (ite <= iteMax) {
2039 ite++;
2040 phi = lammSmqrUnif(Z, Y, Lambda, betaNew, phi, tau, gamma, p, h, n1, h1);
2041 phi = std::max(phi0, phi / gamma);
2042 if (arma::norm(betaNew - beta, "inf") <= epsilon) {
2043 break;
2044 }
2045 beta = betaNew;
2046 }
2047 int iteT = 1;
2048 // Tightening
2049 arma::vec beta0(p + 1);
2050 while (iteT <= iteTight) {
2051 iteT++;
2052 beta = betaNew;
2053 beta0 = betaNew;
2054 Lambda = cmptLambdaMCP(beta, lambda, p, para);
2055 phi = phi0;
2056 ite = 0;
2057 while (ite <= iteMax) {
2058 ite++;
2059 phi = lammSmqrUnif(Z, Y, Lambda, betaNew, phi, tau, gamma, p, h, n1, h1);
2060 phi = std::max(phi0, phi / gamma);
2061 if (arma::norm(betaNew - beta, "inf") <= epsilon) {
2062 break;
2063 }
2064 beta = betaNew;
2065 }
2066 if (arma::norm(betaNew - beta0, "inf") <= epsilon) {
2067 break;
2068 }
2069 }
2070 return betaNew;
2071 }
2072
2073 // [[Rcpp::export]]
smqrMcpPara(const arma::mat & Z,const arma::vec & Y,const double lambda,const arma::vec & sx1,const double tau,const int p,const double n1,const double h,const double h1,const double h3,const double phi0=0.01,const double gamma=1.2,const double epsilon=0.001,const int iteMax=500,const int iteTight=3,const double para=3)2074 arma::vec smqrMcpPara(const arma::mat& Z, const arma::vec& Y, const double lambda, const arma::vec& sx1, const double tau, const int p, const double n1,
2075 const double h, const double h1, const double h3, const double phi0 = 0.01, const double gamma = 1.2, const double epsilon = 0.001,
2076 const int iteMax = 500, const int iteTight = 3, const double para = 3) {
2077 arma::vec beta = lasso(Z, Y, lambda, tau, p, n1, phi0, gamma, epsilon, iteMax);
2078 arma::vec betaNew = beta;
2079 // Contraction
2080 arma::vec Lambda = cmptLambdaMCP(beta, lambda, p, para);
2081 double phi = phi0;
2082 int ite = 0;
2083 while (ite <= iteMax) {
2084 ite++;
2085 phi = lammSmqrPara(Z, Y, Lambda, betaNew, phi, tau, gamma, p, h, n1, h1, h3);
2086 phi = std::max(phi0, phi / gamma);
2087 if (arma::norm(betaNew - beta, "inf") <= epsilon) {
2088 break;
2089 }
2090 beta = betaNew;
2091 }
2092 int iteT = 1;
2093 // Tightening
2094 arma::vec beta0(p + 1);
2095 while (iteT <= iteTight) {
2096 iteT++;
2097 beta = betaNew;
2098 beta0 = betaNew;
2099 Lambda = cmptLambdaMCP(beta, lambda, p, para);
2100 phi = phi0;
2101 ite = 0;
2102 while (ite <= iteMax) {
2103 ite++;
2104 phi = lammSmqrPara(Z, Y, Lambda, betaNew, phi, tau, gamma, p, h, n1, h1, h3);
2105 phi = std::max(phi0, phi / gamma);
2106 if (arma::norm(betaNew - beta, "inf") <= epsilon) {
2107 break;
2108 }
2109 beta = betaNew;
2110 }
2111 if (arma::norm(betaNew - beta0, "inf") <= epsilon) {
2112 break;
2113 }
2114 }
2115 return betaNew;
2116 }
2117
2118 // [[Rcpp::export]]
smqrMcpTrian(const arma::mat & Z,const arma::vec & Y,const double lambda,const arma::vec & sx1,const double tau,const int p,const double n1,const double h,const double h1,const double h2,const double phi0=0.01,const double gamma=1.2,const double epsilon=0.001,const int iteMax=500,const int iteTight=3,const double para=3)2119 arma::vec smqrMcpTrian(const arma::mat& Z, const arma::vec& Y, const double lambda, const arma::vec& sx1, const double tau, const int p, const double n1,
2120 const double h, const double h1, const double h2, const double phi0 = 0.01, const double gamma = 1.2, const double epsilon = 0.001,
2121 const int iteMax = 500, const int iteTight = 3, const double para = 3) {
2122 arma::vec beta = lasso(Z, Y, lambda, tau, p, n1, phi0, gamma, epsilon, iteMax);
2123 arma::vec betaNew = beta;
2124 // Contraction
2125 arma::vec Lambda = cmptLambdaMCP(beta, lambda, p, para);
2126 double phi = phi0;
2127 int ite = 0;
2128 while (ite <= iteMax) {
2129 ite++;
2130 phi = lammSmqrTrian(Z, Y, Lambda, betaNew, phi, tau, gamma, p, h, n1, h1, h2);
2131 phi = std::max(phi0, phi / gamma);
2132 if (arma::norm(betaNew - beta, "inf") <= epsilon) {
2133 break;
2134 }
2135 beta = betaNew;
2136 }
2137 int iteT = 1;
2138 // Tightening
2139 arma::vec beta0(p + 1);
2140 while (iteT <= iteTight) {
2141 iteT++;
2142 beta = betaNew;
2143 beta0 = betaNew;
2144 Lambda = cmptLambdaMCP(beta, lambda, p, para);
2145 phi = phi0;
2146 ite = 0;
2147 while (ite <= iteMax) {
2148 ite++;
2149 phi = lammSmqrTrian(Z, Y, Lambda, betaNew, phi, tau, gamma, p, h, n1, h1, h2);
2150 phi = std::max(phi0, phi / gamma);
2151 if (arma::norm(betaNew - beta, "inf") <= epsilon) {
2152 break;
2153 }
2154 beta = betaNew;
2155 }
2156 if (arma::norm(betaNew - beta0, "inf") <= epsilon) {
2157 break;
2158 }
2159 }
2160 return betaNew;
2161 }
2162
2163 // high-dim conquer with a specified lambda
2164 // [[Rcpp::export]]
conquerHdGauss(const arma::mat & X,arma::vec Y,const double lambda,const double tau,const double h,const int type=1,const double phi0=0.01,const double gamma=1.2,const double epsilon=0.001,const int iteMax=500,const int iteTight=3,const double para=3.7)2165 arma::vec conquerHdGauss(const arma::mat& X, arma::vec Y, const double lambda, const double tau, const double h, const int type = 1, const double phi0 = 0.01,
2166 const double gamma = 1.2, const double epsilon = 0.001, const int iteMax = 500, const int iteTight = 3, const double para = 3.7) {
2167 const int n = X.n_rows, p = X.n_cols;
2168 const double h1 = 1.0 / h, h2 = 1.0 / (h * h);
2169 arma::rowvec mx = arma::mean(X, 0);
2170 arma::vec sx1 = 1.0 / arma::stddev(X, 0, 0).t();
2171 arma::mat Z = arma::join_rows(arma::ones(n), standardize(X, mx, sx1, p));
2172 double my = arma::mean(Y);
2173 Y -= my;
2174 arma::vec betaHat(p + 1);
2175 if (type == 1) {
2176 betaHat = smqrLassoGauss(Z, Y, lambda, sx1, tau, p, 1.0 / n, h, h1, h2, phi0, gamma, epsilon, iteMax);
2177 } else if (type == 2) {
2178 betaHat = smqrScadGauss(Z, Y, lambda, sx1, tau, p, 1.0 / n, h, h1, h2, phi0, gamma, epsilon, iteMax, iteTight, para);
2179 } else {
2180 betaHat = smqrMcpGauss(Z, Y, lambda, sx1, tau, p, 1.0 / n, h, h1, h2, phi0, gamma, epsilon, iteMax, iteTight, para);
2181 }
2182 betaHat.rows(1, p) %= sx1;
2183 betaHat(0) += my - arma::as_scalar(mx * betaHat.rows(1, p));
2184 return betaHat;
2185 }
2186
2187 // [[Rcpp::export]]
conquerHdLogistic(const arma::mat & X,arma::vec Y,const double lambda,const double tau,const double h,const int type=1,const double phi0=0.01,const double gamma=1.2,const double epsilon=0.001,const int iteMax=500,const int iteTight=3,const double para=3.7)2188 arma::vec conquerHdLogistic(const arma::mat& X, arma::vec Y, const double lambda, const double tau, const double h, const int type = 1, const double phi0 = 0.01,
2189 const double gamma = 1.2, const double epsilon = 0.001, const int iteMax = 500, const int iteTight = 3, const double para = 3.7) {
2190 const int n = X.n_rows, p = X.n_cols;
2191 const double h1 = 1.0 / h;
2192 arma::rowvec mx = arma::mean(X, 0);
2193 arma::vec sx1 = 1.0 / arma::stddev(X, 0, 0).t();
2194 arma::mat Z = arma::join_rows(arma::ones(n), standardize(X, mx, sx1, p));
2195 double my = arma::mean(Y);
2196 Y -= my;
2197 arma::vec betaHat(p + 1);
2198 if (type == 1) {
2199 betaHat = smqrLassoLogistic(Z, Y, lambda, sx1, tau, p, 1.0 / n, h, h1, phi0, gamma, epsilon, iteMax);
2200 } else if (type == 2) {
2201 betaHat = smqrScadLogistic(Z, Y, lambda, sx1, tau, p, 1.0 / n, h, h1, phi0, gamma, epsilon, iteMax, iteTight, para);
2202 } else {
2203 betaHat = smqrMcpLogistic(Z, Y, lambda, sx1, tau, p, 1.0 / n, h, h1, phi0, gamma, epsilon, iteMax, iteTight, para);
2204 }
2205 betaHat.rows(1, p) %= sx1;
2206 betaHat(0) += my - arma::as_scalar(mx * betaHat.rows(1, p));
2207 return betaHat;
2208 }
2209
2210 // [[Rcpp::export]]
conquerHdUnif(const arma::mat & X,arma::vec Y,const double lambda,const double tau,const double h,const int type=1,const double phi0=0.01,const double gamma=1.2,const double epsilon=0.001,const int iteMax=500,const int iteTight=3,const double para=3.7)2211 arma::vec conquerHdUnif(const arma::mat& X, arma::vec Y, const double lambda, const double tau, const double h, const int type = 1, const double phi0 = 0.01,
2212 const double gamma = 1.2, const double epsilon = 0.001, const int iteMax = 500, const int iteTight = 3, const double para = 3.7) {
2213 const int n = X.n_rows, p = X.n_cols;
2214 const double h1 = 1.0 / h;
2215 arma::rowvec mx = arma::mean(X, 0);
2216 arma::vec sx1 = 1.0 / arma::stddev(X, 0, 0).t();
2217 arma::mat Z = arma::join_rows(arma::ones(n), standardize(X, mx, sx1, p));
2218 double my = arma::mean(Y);
2219 Y -= my;
2220 arma::vec betaHat(p + 1);
2221 if (type == 1) {
2222 betaHat = smqrLassoUnif(Z, Y, lambda, sx1, tau, p, 1.0 / n, h, h1, phi0, gamma, epsilon, iteMax);
2223 } else if (type == 2) {
2224 betaHat = smqrScadUnif(Z, Y, lambda, sx1, tau, p, 1.0 / n, h, h1, phi0, gamma, epsilon, iteMax, iteTight, para);
2225 } else {
2226 betaHat = smqrMcpUnif(Z, Y, lambda, sx1, tau, p, 1.0 / n, h, h1, phi0, gamma, epsilon, iteMax, iteTight, para);
2227 }
2228 betaHat.rows(1, p) %= sx1;
2229 betaHat(0) += my - arma::as_scalar(mx * betaHat.rows(1, p));
2230 return betaHat;
2231 }
2232
2233 // [[Rcpp::export]]
conquerHdPara(const arma::mat & X,arma::vec Y,const double lambda,const double tau,const double h,const int type=1,const double phi0=0.01,const double gamma=1.2,const double epsilon=0.001,const int iteMax=500,const int iteTight=3,const double para=3.7)2234 arma::vec conquerHdPara(const arma::mat& X, arma::vec Y, const double lambda, const double tau, const double h, const int type = 1, const double phi0 = 0.01,
2235 const double gamma = 1.2, const double epsilon = 0.001, const int iteMax = 500, const int iteTight = 3, const double para = 3.7) {
2236 const int n = X.n_rows, p = X.n_cols;
2237 const double h1 = 1.0 / h, h3 = 1.0 / (h * h * h);
2238 arma::rowvec mx = arma::mean(X, 0);
2239 arma::vec sx1 = 1.0 / arma::stddev(X, 0, 0).t();
2240 arma::mat Z = arma::join_rows(arma::ones(n), standardize(X, mx, sx1, p));
2241 double my = arma::mean(Y);
2242 Y -= my;
2243 arma::vec betaHat(p + 1);
2244 if (type == 1) {
2245 betaHat = smqrLassoPara(Z, Y, lambda, sx1, tau, p, 1.0 / n, h, h1, h3, phi0, gamma, epsilon, iteMax);
2246 } else if (type == 2) {
2247 betaHat = smqrScadPara(Z, Y, lambda, sx1, tau, p, 1.0 / n, h, h1, h3, phi0, gamma, epsilon, iteMax, iteTight, para);
2248 } else {
2249 betaHat = smqrMcpPara(Z, Y, lambda, sx1, tau, p, 1.0 / n, h, h1, h3, phi0, gamma, epsilon, iteMax, iteTight, para);
2250 }
2251 betaHat.rows(1, p) %= sx1;
2252 betaHat(0) += my - arma::as_scalar(mx * betaHat.rows(1, p));
2253 return betaHat;
2254 }
2255
2256 // [[Rcpp::export]]
conquerHdTrian(const arma::mat & X,arma::vec Y,const double lambda,const double tau,const double h,const int type=1,const double phi0=0.01,const double gamma=1.2,const double epsilon=0.001,const int iteMax=500,const int iteTight=3,const double para=3.7)2257 arma::vec conquerHdTrian(const arma::mat& X, arma::vec Y, const double lambda, const double tau, const double h, const int type = 1, const double phi0 = 0.01,
2258 const double gamma = 1.2, const double epsilon = 0.001, const int iteMax = 500, const int iteTight = 3, const double para = 3.7) {
2259 const int n = X.n_rows, p = X.n_cols;
2260 const double h1 = 1.0 / h, h2 = 1.0 / (h * h);
2261 arma::rowvec mx = arma::mean(X, 0);
2262 arma::vec sx1 = 1.0 / arma::stddev(X, 0, 0).t();
2263 arma::mat Z = arma::join_rows(arma::ones(n), standardize(X, mx, sx1, p));
2264 double my = arma::mean(Y);
2265 Y -= my;
2266 arma::vec betaHat(p + 1);
2267 if (type == 1) {
2268 betaHat = smqrLassoTrian(Z, Y, lambda, sx1, tau, p, 1.0 / n, h, h1, h2, phi0, gamma, epsilon, iteMax);
2269 } else if (type == 2) {
2270 betaHat = smqrScadTrian(Z, Y, lambda, sx1, tau, p, 1.0 / n, h, h1, h2, phi0, gamma, epsilon, iteMax, iteTight, para);
2271 } else {
2272 betaHat = smqrMcpTrian(Z, Y, lambda, sx1, tau, p, 1.0 / n, h, h1, h2, phi0, gamma, epsilon, iteMax, iteTight, para);
2273 }
2274 betaHat.rows(1, p) %= sx1;
2275 betaHat(0) += my - arma::as_scalar(mx * betaHat.rows(1, p));
2276 return betaHat;
2277 }
2278
2279 // cross-validation, no warm-start, the range of lambda is guided by the simulation-based mathod of Belloni & Chernozhukov (2011), AOS
2280 // [[Rcpp::export]]
lossQr(const arma::mat & Z,const arma::vec & Y,const arma::vec & beta,const double tau)2281 double lossQr(const arma::mat& Z, const arma::vec& Y, const arma::vec& beta, const double tau) {
2282 arma::vec res = Y - Z * beta;
2283 double rst = 0.0;
2284 for (int i = 0; i < res.size(); i++) {
2285 rst += res(i) >= 0 ? tau * res(i) : (tau - 1) * res(i);
2286 }
2287 return rst;
2288 }
2289
2290 // [[Rcpp::export]]
cvSmqrLassoGauss(const arma::mat & X,arma::vec Y,const arma::vec & lambdaSeq,const arma::vec & folds,const double tau,const int kfolds,const double h,const double phi0=0.01,const double gamma=1.2,const double epsilon=0.001,const int iteMax=500)2291 Rcpp::List cvSmqrLassoGauss(const arma::mat& X, arma::vec Y, const arma::vec& lambdaSeq, const arma::vec& folds, const double tau, const int kfolds,
2292 const double h, const double phi0 = 0.01, const double gamma = 1.2, const double epsilon = 0.001, const int iteMax = 500) {
2293 const int n = X.n_rows, p = X.n_cols, nlambda = lambdaSeq.size();
2294 const double h1 = 1.0 / h, h2 = 1.0 / (h * h);
2295 arma::vec betaHat(p + 1);
2296 arma::vec mse = arma::zeros(nlambda);
2297 arma::rowvec mx = arma::mean(X, 0);
2298 arma::vec sx1 = 1.0 / arma::stddev(X, 0, 0).t();
2299 arma::mat Z = arma::join_rows(arma::ones(n), standardize(X, mx, sx1, p));
2300 double my = arma::mean(Y);
2301 Y -= my;
2302 for (int j = 1; j <= kfolds; j++) {
2303 arma::uvec idx = arma::find(folds == j);
2304 arma::uvec idxComp = arma::find(folds != j);
2305 double n1Train = 1.0 / idxComp.size();
2306 arma::mat trainZ = Z.rows(idxComp), testZ = Z.rows(idx);
2307 arma::vec trainY = Y.rows(idxComp), testY = Y.rows(idx);
2308 for (int i = 0; i < nlambda; i++) {
2309 betaHat = smqrLassoGauss(trainZ, trainY, lambdaSeq(i), sx1, tau, p, n1Train, h, h1, h2, phi0, gamma, epsilon, iteMax);
2310 mse(i) += arma::accu(lossQr(testZ, testY, betaHat, tau));
2311 }
2312 }
2313 arma::uword cvIdx = arma::index_min(mse);
2314 betaHat = smqrLassoGauss(Z, Y, lambdaSeq(cvIdx), sx1, tau, p, 1.0 / n, h, h1, h2, phi0, gamma, epsilon, iteMax);
2315 betaHat.rows(1, p) %= sx1;
2316 betaHat(0) += my - arma::as_scalar(mx * betaHat.rows(1, p));
2317 return Rcpp::List::create(Rcpp::Named("coeff") = betaHat, Rcpp::Named("lambda") = lambdaSeq(cvIdx));
2318 }
2319
2320 // [[Rcpp::export]]
cvSmqrScadGauss(const arma::mat & X,arma::vec Y,const arma::vec & lambdaSeq,const arma::vec & folds,const double tau,const int kfolds,const double h,const double phi0=0.01,const double gamma=1.2,const double epsilon=0.001,const int iteMax=500,const int iteTight=3,const double para=3.7)2321 Rcpp::List cvSmqrScadGauss(const arma::mat& X, arma::vec Y, const arma::vec& lambdaSeq, const arma::vec& folds, const double tau, const int kfolds,
2322 const double h, const double phi0 = 0.01, const double gamma = 1.2, const double epsilon = 0.001, const int iteMax = 500,
2323 const int iteTight = 3, const double para = 3.7) {
2324 const int n = X.n_rows, p = X.n_cols, nlambda = lambdaSeq.size();
2325 const double h1 = 1.0 / h, h2 = 1.0 / (h * h);
2326 arma::vec betaHat(p + 1);
2327 arma::vec mse = arma::zeros(nlambda);
2328 arma::rowvec mx = arma::mean(X, 0);
2329 arma::vec sx1 = 1.0 / arma::stddev(X, 0, 0).t();
2330 arma::mat Z = arma::join_rows(arma::ones(n), standardize(X, mx, sx1, p));
2331 double my = arma::mean(Y);
2332 Y -= my;
2333 for (int j = 1; j <= kfolds; j++) {
2334 arma::uvec idx = arma::find(folds == j);
2335 arma::uvec idxComp = arma::find(folds != j);
2336 double n1Train = 1.0 / idxComp.size();
2337 arma::mat trainZ = Z.rows(idxComp), testZ = Z.rows(idx);
2338 arma::vec trainY = Y.rows(idxComp), testY = Y.rows(idx);
2339 for (int i = 0; i < nlambda; i++) {
2340 betaHat = smqrScadGauss(trainZ, trainY, lambdaSeq(i), sx1, tau, p, n1Train, h, h1, h2, phi0, gamma, epsilon, iteMax, iteTight, para);
2341 mse(i) += arma::accu(lossQr(testZ, testY, betaHat, tau));
2342 }
2343 }
2344 arma::uword cvIdx = arma::index_min(mse);
2345 betaHat = smqrScadGauss(Z, Y, lambdaSeq(cvIdx), sx1, tau, p, 1.0 / n, h, h1, h2, phi0, gamma, epsilon, iteMax, iteTight, para);
2346 betaHat.rows(1, p) %= sx1;
2347 betaHat(0) += my - arma::as_scalar(mx * betaHat.rows(1, p));
2348 return Rcpp::List::create(Rcpp::Named("coeff") = betaHat, Rcpp::Named("lambda") = lambdaSeq(cvIdx));
2349 }
2350
2351 // [[Rcpp::export]]
cvSmqrMcpGauss(const arma::mat & X,arma::vec Y,const arma::vec & lambdaSeq,const arma::vec & folds,const double tau,const int kfolds,const double h,const double phi0=0.01,const double gamma=1.5,const double epsilon=0.001,const int iteMax=500,const int iteTight=3,const double para=3)2352 Rcpp::List cvSmqrMcpGauss(const arma::mat& X, arma::vec Y, const arma::vec& lambdaSeq, const arma::vec& folds, const double tau, const int kfolds,
2353 const double h, const double phi0 = 0.01, const double gamma = 1.5, const double epsilon = 0.001, const int iteMax = 500,
2354 const int iteTight = 3, const double para = 3) {
2355 const int n = X.n_rows, p = X.n_cols, nlambda = lambdaSeq.size();
2356 const double h1 = 1.0 / h, h2 = 1.0 / (h * h);
2357 arma::vec betaHat(p + 1);
2358 arma::vec mse = arma::zeros(nlambda);
2359 arma::rowvec mx = arma::mean(X, 0);
2360 arma::vec sx1 = 1.0 / arma::stddev(X, 0, 0).t();
2361 arma::mat Z = arma::join_rows(arma::ones(n), standardize(X, mx, sx1, p));
2362 double my = arma::mean(Y);
2363 Y -= my;
2364 for (int j = 1; j <= kfolds; j++) {
2365 arma::uvec idx = arma::find(folds == j);
2366 arma::uvec idxComp = arma::find(folds != j);
2367 double n1Train = 1.0 / idxComp.size();
2368 arma::mat trainZ = Z.rows(idxComp), testZ = Z.rows(idx);
2369 arma::vec trainY = Y.rows(idxComp), testY = Y.rows(idx);
2370 for (int i = 0; i < nlambda; i++) {
2371 betaHat = smqrMcpGauss(trainZ, trainY, lambdaSeq(i), sx1, tau, p, n1Train, h, h1, h2, phi0, gamma, epsilon, iteMax, iteTight, para);
2372 mse(i) += arma::accu(lossQr(testZ, testY, betaHat, tau));
2373 }
2374 }
2375 arma::uword cvIdx = arma::index_min(mse);
2376 betaHat = smqrMcpGauss(Z, Y, lambdaSeq(cvIdx), sx1, tau, p, 1.0 / n, h, h1, h2, phi0, gamma, epsilon, iteMax, iteTight, para);
2377 betaHat.rows(1, p) %= sx1;
2378 betaHat(0) += my - arma::as_scalar(mx * betaHat.rows(1, p));
2379 return Rcpp::List::create(Rcpp::Named("coeff") = betaHat, Rcpp::Named("lambda") = lambdaSeq(cvIdx));
2380 }
2381
2382 // [[Rcpp::export]]
cvSmqrLassoLogistic(const arma::mat & X,arma::vec Y,const arma::vec & lambdaSeq,const arma::vec & folds,const double tau,const int kfolds,const double h,const double phi0=0.01,const double gamma=1.2,const double epsilon=0.001,const int iteMax=500)2383 Rcpp::List cvSmqrLassoLogistic(const arma::mat& X, arma::vec Y, const arma::vec& lambdaSeq, const arma::vec& folds, const double tau, const int kfolds,
2384 const double h, const double phi0 = 0.01, const double gamma = 1.2, const double epsilon = 0.001, const int iteMax = 500) {
2385 const int n = X.n_rows, p = X.n_cols, nlambda = lambdaSeq.size();
2386 const double h1 = 1.0 / h;
2387 arma::vec betaHat(p + 1);
2388 arma::vec mse = arma::zeros(nlambda);
2389 arma::rowvec mx = arma::mean(X, 0);
2390 arma::vec sx1 = 1.0 / arma::stddev(X, 0, 0).t();
2391 arma::mat Z = arma::join_rows(arma::ones(n), standardize(X, mx, sx1, p));
2392 double my = arma::mean(Y);
2393 Y -= my;
2394 for (int j = 1; j <= kfolds; j++) {
2395 arma::uvec idx = arma::find(folds == j);
2396 arma::uvec idxComp = arma::find(folds != j);
2397 double n1Train = 1.0 / idxComp.size();
2398 arma::mat trainZ = Z.rows(idxComp), testZ = Z.rows(idx);
2399 arma::vec trainY = Y.rows(idxComp), testY = Y.rows(idx);
2400 for (int i = 0; i < nlambda; i++) {
2401 betaHat = smqrLassoLogistic(trainZ, trainY, lambdaSeq(i), sx1, tau, p, n1Train, h, h1, phi0, gamma, epsilon, iteMax);
2402 mse(i) += arma::accu(lossQr(testZ, testY, betaHat, tau));
2403 }
2404 }
2405 arma::uword cvIdx = arma::index_min(mse);
2406 betaHat = smqrLassoLogistic(Z, Y, lambdaSeq(cvIdx), sx1, tau, p, 1.0 / n, h, h1, phi0, gamma, epsilon, iteMax);
2407 betaHat.rows(1, p) %= sx1;
2408 betaHat(0) += my - arma::as_scalar(mx * betaHat.rows(1, p));
2409 return Rcpp::List::create(Rcpp::Named("coeff") = betaHat, Rcpp::Named("lambda") = lambdaSeq(cvIdx));
2410 }
2411
2412 // [[Rcpp::export]]
cvSmqrScadLogistic(const arma::mat & X,arma::vec Y,const arma::vec & lambdaSeq,const arma::vec & folds,const double tau,const int kfolds,const double h,const double phi0=0.01,const double gamma=1.2,const double epsilon=0.001,const int iteMax=500,const int iteTight=3,const double para=3.7)2413 Rcpp::List cvSmqrScadLogistic(const arma::mat& X, arma::vec Y, const arma::vec& lambdaSeq, const arma::vec& folds, const double tau, const int kfolds,
2414 const double h, const double phi0 = 0.01, const double gamma = 1.2, const double epsilon = 0.001, const int iteMax = 500,
2415 const int iteTight = 3, const double para = 3.7) {
2416 const int n = X.n_rows, p = X.n_cols, nlambda = lambdaSeq.size();
2417 const double h1 = 1.0 / h;
2418 arma::vec betaHat(p + 1);
2419 arma::vec mse = arma::zeros(nlambda);
2420 arma::rowvec mx = arma::mean(X, 0);
2421 arma::vec sx1 = 1.0 / arma::stddev(X, 0, 0).t();
2422 arma::mat Z = arma::join_rows(arma::ones(n), standardize(X, mx, sx1, p));
2423 double my = arma::mean(Y);
2424 Y -= my;
2425 for (int j = 1; j <= kfolds; j++) {
2426 arma::uvec idx = arma::find(folds == j);
2427 arma::uvec idxComp = arma::find(folds != j);
2428 double n1Train = 1.0 / idxComp.size();
2429 arma::mat trainZ = Z.rows(idxComp), testZ = Z.rows(idx);
2430 arma::vec trainY = Y.rows(idxComp), testY = Y.rows(idx);
2431 for (int i = 0; i < nlambda; i++) {
2432 betaHat = smqrScadLogistic(trainZ, trainY, lambdaSeq(i), sx1, tau, p, n1Train, h, h1, phi0, gamma, epsilon, iteMax, iteTight, para);
2433 mse(i) += arma::accu(lossQr(testZ, testY, betaHat, tau));
2434 }
2435 }
2436 arma::uword cvIdx = arma::index_min(mse);
2437 betaHat = smqrScadLogistic(Z, Y, lambdaSeq(cvIdx), sx1, tau, p, 1.0 / n, h, h1, phi0, gamma, epsilon, iteMax, iteTight, para);
2438 betaHat.rows(1, p) %= sx1;
2439 betaHat(0) += my - arma::as_scalar(mx * betaHat.rows(1, p));
2440 return Rcpp::List::create(Rcpp::Named("coeff") = betaHat, Rcpp::Named("lambda") = lambdaSeq(cvIdx));
2441 }
2442
2443 // [[Rcpp::export]]
cvSmqrMcpLogistic(const arma::mat & X,arma::vec Y,const arma::vec & lambdaSeq,const arma::vec & folds,const double tau,const int kfolds,const double h,const double phi0=0.01,const double gamma=1.5,const double epsilon=0.001,const int iteMax=500,const int iteTight=3,const double para=3)2444 Rcpp::List cvSmqrMcpLogistic(const arma::mat& X, arma::vec Y, const arma::vec& lambdaSeq, const arma::vec& folds, const double tau, const int kfolds,
2445 const double h, const double phi0 = 0.01, const double gamma = 1.5, const double epsilon = 0.001, const int iteMax = 500,
2446 const int iteTight = 3, const double para = 3) {
2447 const int n = X.n_rows, p = X.n_cols, nlambda = lambdaSeq.size();
2448 const double h1 = 1.0 / h;
2449 arma::vec betaHat(p + 1);
2450 arma::vec mse = arma::zeros(nlambda);
2451 arma::rowvec mx = arma::mean(X, 0);
2452 arma::vec sx1 = 1.0 / arma::stddev(X, 0, 0).t();
2453 arma::mat Z = arma::join_rows(arma::ones(n), standardize(X, mx, sx1, p));
2454 double my = arma::mean(Y);
2455 Y -= my;
2456 for (int j = 1; j <= kfolds; j++) {
2457 arma::uvec idx = arma::find(folds == j);
2458 arma::uvec idxComp = arma::find(folds != j);
2459 double n1Train = 1.0 / idxComp.size();
2460 arma::mat trainZ = Z.rows(idxComp), testZ = Z.rows(idx);
2461 arma::vec trainY = Y.rows(idxComp), testY = Y.rows(idx);
2462 for (int i = 0; i < nlambda; i++) {
2463 betaHat = smqrMcpLogistic(trainZ, trainY, lambdaSeq(i), sx1, tau, p, n1Train, h, h1, phi0, gamma, epsilon, iteMax, iteTight, para);
2464 mse(i) += arma::accu(lossQr(testZ, testY, betaHat, tau));
2465 }
2466 }
2467 arma::uword cvIdx = arma::index_min(mse);
2468 betaHat = smqrMcpLogistic(Z, Y, lambdaSeq(cvIdx), sx1, tau, p, 1.0 / n, h, h1, phi0, gamma, epsilon, iteMax, iteTight, para);
2469 betaHat.rows(1, p) %= sx1;
2470 betaHat(0) += my - arma::as_scalar(mx * betaHat.rows(1, p));
2471 return Rcpp::List::create(Rcpp::Named("coeff") = betaHat, Rcpp::Named("lambda") = lambdaSeq(cvIdx));
2472 }
2473
2474 // [[Rcpp::export]]
cvSmqrLassoUnif(const arma::mat & X,arma::vec Y,const arma::vec & lambdaSeq,const arma::vec & folds,const double tau,const int kfolds,const double h,const double phi0=0.01,const double gamma=1.2,const double epsilon=0.001,const int iteMax=500)2475 Rcpp::List cvSmqrLassoUnif(const arma::mat& X, arma::vec Y, const arma::vec& lambdaSeq, const arma::vec& folds, const double tau, const int kfolds,
2476 const double h, const double phi0 = 0.01, const double gamma = 1.2, const double epsilon = 0.001, const int iteMax = 500) {
2477 const int n = X.n_rows, p = X.n_cols, nlambda = lambdaSeq.size();
2478 const double h1 = 1.0 / h;
2479 arma::vec betaHat(p + 1);
2480 arma::vec mse = arma::zeros(nlambda);
2481 arma::rowvec mx = arma::mean(X, 0);
2482 arma::vec sx1 = 1.0 / arma::stddev(X, 0, 0).t();
2483 arma::mat Z = arma::join_rows(arma::ones(n), standardize(X, mx, sx1, p));
2484 double my = arma::mean(Y);
2485 Y -= my;
2486 for (int j = 1; j <= kfolds; j++) {
2487 arma::uvec idx = arma::find(folds == j);
2488 arma::uvec idxComp = arma::find(folds != j);
2489 double n1Train = 1.0 / idxComp.size();
2490 arma::mat trainZ = Z.rows(idxComp), testZ = Z.rows(idx);
2491 arma::vec trainY = Y.rows(idxComp), testY = Y.rows(idx);
2492 for (int i = 0; i < nlambda; i++) {
2493 betaHat = smqrLassoUnif(trainZ, trainY, lambdaSeq(i), sx1, tau, p, n1Train, h, h1, phi0, gamma, epsilon, iteMax);
2494 mse(i) += arma::accu(lossQr(testZ, testY, betaHat, tau));
2495 }
2496 }
2497 arma::uword cvIdx = arma::index_min(mse);
2498 betaHat = smqrLassoUnif(Z, Y, lambdaSeq(cvIdx), sx1, tau, p, 1.0 / n, h, h1, phi0, gamma, epsilon, iteMax);
2499 betaHat.rows(1, p) %= sx1;
2500 betaHat(0) += my - arma::as_scalar(mx * betaHat.rows(1, p));
2501 return Rcpp::List::create(Rcpp::Named("coeff") = betaHat, Rcpp::Named("lambda") = lambdaSeq(cvIdx));
2502 }
2503
2504 // [[Rcpp::export]]
cvSmqrScadUnif(const arma::mat & X,arma::vec Y,const arma::vec & lambdaSeq,const arma::vec & folds,const double tau,const int kfolds,const double h,const double phi0=0.01,const double gamma=1.2,const double epsilon=0.001,const int iteMax=500,const int iteTight=3,const double para=3.7)2505 Rcpp::List cvSmqrScadUnif(const arma::mat& X, arma::vec Y, const arma::vec& lambdaSeq, const arma::vec& folds, const double tau, const int kfolds,
2506 const double h, const double phi0 = 0.01, const double gamma = 1.2, const double epsilon = 0.001, const int iteMax = 500,
2507 const int iteTight = 3, const double para = 3.7) {
2508 const int n = X.n_rows, p = X.n_cols, nlambda = lambdaSeq.size();
2509 const double h1 = 1.0 / h;
2510 arma::vec betaHat(p + 1);
2511 arma::vec mse = arma::zeros(nlambda);
2512 arma::rowvec mx = arma::mean(X, 0);
2513 arma::vec sx1 = 1.0 / arma::stddev(X, 0, 0).t();
2514 arma::mat Z = arma::join_rows(arma::ones(n), standardize(X, mx, sx1, p));
2515 double my = arma::mean(Y);
2516 Y -= my;
2517 for (int j = 1; j <= kfolds; j++) {
2518 arma::uvec idx = arma::find(folds == j);
2519 arma::uvec idxComp = arma::find(folds != j);
2520 double n1Train = 1.0 / idxComp.size();
2521 arma::mat trainZ = Z.rows(idxComp), testZ = Z.rows(idx);
2522 arma::vec trainY = Y.rows(idxComp), testY = Y.rows(idx);
2523 for (int i = 0; i < nlambda; i++) {
2524 betaHat = smqrScadUnif(trainZ, trainY, lambdaSeq(i), sx1, tau, p, n1Train, h, h1, phi0, gamma, epsilon, iteMax, iteTight, para);
2525 mse(i) += arma::accu(lossQr(testZ, testY, betaHat, tau));
2526 }
2527 }
2528 arma::uword cvIdx = arma::index_min(mse);
2529 betaHat = smqrScadUnif(Z, Y, lambdaSeq(cvIdx), sx1, tau, p, 1.0 / n, h, h1, phi0, gamma, epsilon, iteMax, iteTight, para);
2530 betaHat.rows(1, p) %= sx1;
2531 betaHat(0) += my - arma::as_scalar(mx * betaHat.rows(1, p));
2532 return Rcpp::List::create(Rcpp::Named("coeff") = betaHat, Rcpp::Named("lambda") = lambdaSeq(cvIdx));
2533 }
2534
2535 // [[Rcpp::export]]
cvSmqrMcpUnif(const arma::mat & X,arma::vec Y,const arma::vec & lambdaSeq,const arma::vec & folds,const double tau,const int kfolds,const double h,const double phi0=0.01,const double gamma=1.5,const double epsilon=0.001,const int iteMax=500,const int iteTight=3,const double para=3)2536 Rcpp::List cvSmqrMcpUnif(const arma::mat& X, arma::vec Y, const arma::vec& lambdaSeq, const arma::vec& folds, const double tau, const int kfolds,
2537 const double h, const double phi0 = 0.01, const double gamma = 1.5, const double epsilon = 0.001, const int iteMax = 500,
2538 const int iteTight = 3, const double para = 3) {
2539 const int n = X.n_rows, p = X.n_cols, nlambda = lambdaSeq.size();
2540 const double h1 = 1.0 / h;
2541 arma::vec betaHat(p + 1);
2542 arma::vec mse = arma::zeros(nlambda);
2543 arma::rowvec mx = arma::mean(X, 0);
2544 arma::vec sx1 = 1.0 / arma::stddev(X, 0, 0).t();
2545 arma::mat Z = arma::join_rows(arma::ones(n), standardize(X, mx, sx1, p));
2546 double my = arma::mean(Y);
2547 Y -= my;
2548 for (int j = 1; j <= kfolds; j++) {
2549 arma::uvec idx = arma::find(folds == j);
2550 arma::uvec idxComp = arma::find(folds != j);
2551 double n1Train = 1.0 / idxComp.size();
2552 arma::mat trainZ = Z.rows(idxComp), testZ = Z.rows(idx);
2553 arma::vec trainY = Y.rows(idxComp), testY = Y.rows(idx);
2554 for (int i = 0; i < nlambda; i++) {
2555 betaHat = smqrMcpUnif(trainZ, trainY, lambdaSeq(i), sx1, tau, p, n1Train, h, h1, phi0, gamma, epsilon, iteMax, iteTight, para);
2556 mse(i) += arma::accu(lossQr(testZ, testY, betaHat, tau));
2557 }
2558 }
2559 arma::uword cvIdx = arma::index_min(mse);
2560 betaHat = smqrMcpUnif(Z, Y, lambdaSeq(cvIdx), sx1, tau, p, 1.0 / n, h, h1, phi0, gamma, epsilon, iteMax, iteTight, para);
2561 betaHat.rows(1, p) %= sx1;
2562 betaHat(0) += my - arma::as_scalar(mx * betaHat.rows(1, p));
2563 return Rcpp::List::create(Rcpp::Named("coeff") = betaHat, Rcpp::Named("lambda") = lambdaSeq(cvIdx));
2564 }
2565
2566 // [[Rcpp::export]]
cvSmqrLassoPara(const arma::mat & X,arma::vec Y,const arma::vec & lambdaSeq,const arma::vec & folds,const double tau,const int kfolds,const double h,const double phi0=0.01,const double gamma=1.2,const double epsilon=0.001,const int iteMax=500)2567 Rcpp::List cvSmqrLassoPara(const arma::mat& X, arma::vec Y, const arma::vec& lambdaSeq, const arma::vec& folds, const double tau, const int kfolds,
2568 const double h, const double phi0 = 0.01, const double gamma = 1.2, const double epsilon = 0.001, const int iteMax = 500) {
2569 const int n = X.n_rows, p = X.n_cols, nlambda = lambdaSeq.size();
2570 const double h1 = 1.0 / h, h3 = 1.0 / (h * h * h);
2571 arma::vec betaHat(p + 1);
2572 arma::vec mse = arma::zeros(nlambda);
2573 arma::rowvec mx = arma::mean(X, 0);
2574 arma::vec sx1 = 1.0 / arma::stddev(X, 0, 0).t();
2575 arma::mat Z = arma::join_rows(arma::ones(n), standardize(X, mx, sx1, p));
2576 double my = arma::mean(Y);
2577 Y -= my;
2578 for (int j = 1; j <= kfolds; j++) {
2579 arma::uvec idx = arma::find(folds == j);
2580 arma::uvec idxComp = arma::find(folds != j);
2581 double n1Train = 1.0 / idxComp.size();
2582 arma::mat trainZ = Z.rows(idxComp), testZ = Z.rows(idx);
2583 arma::vec trainY = Y.rows(idxComp), testY = Y.rows(idx);
2584 for (int i = 0; i < nlambda; i++) {
2585 betaHat = smqrLassoPara(trainZ, trainY, lambdaSeq(i), sx1, tau, p, n1Train, h, h1, h3, phi0, gamma, epsilon, iteMax);
2586 mse(i) += arma::accu(lossQr(testZ, testY, betaHat, tau));
2587 }
2588 }
2589 arma::uword cvIdx = arma::index_min(mse);
2590 betaHat = smqrLassoPara(Z, Y, lambdaSeq(cvIdx), sx1, tau, p, 1.0 / n, h, h1, h3, phi0, gamma, epsilon, iteMax);
2591 betaHat.rows(1, p) %= sx1;
2592 betaHat(0) += my - arma::as_scalar(mx * betaHat.rows(1, p));
2593 return Rcpp::List::create(Rcpp::Named("coeff") = betaHat, Rcpp::Named("lambda") = lambdaSeq(cvIdx));
2594 }
2595
2596 // [[Rcpp::export]]
cvSmqrScadPara(const arma::mat & X,arma::vec Y,const arma::vec & lambdaSeq,const arma::vec & folds,const double tau,const int kfolds,const double h,const double phi0=0.01,const double gamma=1.2,const double epsilon=0.001,const int iteMax=500,const int iteTight=3,const double para=3.7)2597 Rcpp::List cvSmqrScadPara(const arma::mat& X, arma::vec Y, const arma::vec& lambdaSeq, const arma::vec& folds, const double tau, const int kfolds,
2598 const double h, const double phi0 = 0.01, const double gamma = 1.2, const double epsilon = 0.001, const int iteMax = 500,
2599 const int iteTight = 3, const double para = 3.7) {
2600 const int n = X.n_rows, p = X.n_cols, nlambda = lambdaSeq.size();
2601 const double h1 = 1.0 / h, h3 = 1.0 / (h * h * h);
2602 arma::vec betaHat(p + 1);
2603 arma::vec mse = arma::zeros(nlambda);
2604 arma::rowvec mx = arma::mean(X, 0);
2605 arma::vec sx1 = 1.0 / arma::stddev(X, 0, 0).t();
2606 arma::mat Z = arma::join_rows(arma::ones(n), standardize(X, mx, sx1, p));
2607 double my = arma::mean(Y);
2608 Y -= my;
2609 for (int j = 1; j <= kfolds; j++) {
2610 arma::uvec idx = arma::find(folds == j);
2611 arma::uvec idxComp = arma::find(folds != j);
2612 double n1Train = 1.0 / idxComp.size();
2613 arma::mat trainZ = Z.rows(idxComp), testZ = Z.rows(idx);
2614 arma::vec trainY = Y.rows(idxComp), testY = Y.rows(idx);
2615 for (int i = 0; i < nlambda; i++) {
2616 betaHat = smqrScadPara(trainZ, trainY, lambdaSeq(i), sx1, tau, p, n1Train, h, h1, h3, phi0, gamma, epsilon, iteMax, iteTight, para);
2617 mse(i) += arma::accu(lossQr(testZ, testY, betaHat, tau));
2618 }
2619 }
2620 arma::uword cvIdx = arma::index_min(mse);
2621 betaHat = smqrScadPara(Z, Y, lambdaSeq(cvIdx), sx1, tau, p, 1.0 / n, h, h1, h3, phi0, gamma, epsilon, iteMax, iteTight, para);
2622 betaHat.rows(1, p) %= sx1;
2623 betaHat(0) += my - arma::as_scalar(mx * betaHat.rows(1, p));
2624 return Rcpp::List::create(Rcpp::Named("coeff") = betaHat, Rcpp::Named("lambda") = lambdaSeq(cvIdx));
2625 }
2626
2627 // [[Rcpp::export]]
cvSmqrMcpPara(const arma::mat & X,arma::vec Y,const arma::vec & lambdaSeq,const arma::vec & folds,const double tau,const int kfolds,const double h,const double phi0=0.01,const double gamma=1.5,const double epsilon=0.001,const int iteMax=500,const int iteTight=3,const double para=3)2628 Rcpp::List cvSmqrMcpPara(const arma::mat& X, arma::vec Y, const arma::vec& lambdaSeq, const arma::vec& folds, const double tau, const int kfolds,
2629 const double h, const double phi0 = 0.01, const double gamma = 1.5, const double epsilon = 0.001, const int iteMax = 500,
2630 const int iteTight = 3, const double para = 3) {
2631 const int n = X.n_rows, p = X.n_cols, nlambda = lambdaSeq.size();
2632 const double h1 = 1.0 / h, h3 = 1.0 / (h * h * h);
2633 arma::vec betaHat(p + 1);
2634 arma::vec mse = arma::zeros(nlambda);
2635 arma::rowvec mx = arma::mean(X, 0);
2636 arma::vec sx1 = 1.0 / arma::stddev(X, 0, 0).t();
2637 arma::mat Z = arma::join_rows(arma::ones(n), standardize(X, mx, sx1, p));
2638 double my = arma::mean(Y);
2639 Y -= my;
2640 for (int j = 1; j <= kfolds; j++) {
2641 arma::uvec idx = arma::find(folds == j);
2642 arma::uvec idxComp = arma::find(folds != j);
2643 double n1Train = 1.0 / idxComp.size();
2644 arma::mat trainZ = Z.rows(idxComp), testZ = Z.rows(idx);
2645 arma::vec trainY = Y.rows(idxComp), testY = Y.rows(idx);
2646 for (int i = 0; i < nlambda; i++) {
2647 betaHat = smqrMcpPara(trainZ, trainY, lambdaSeq(i), sx1, tau, p, n1Train, h, h1, h3, phi0, gamma, epsilon, iteMax, iteTight, para);
2648 mse(i) += arma::accu(lossQr(testZ, testY, betaHat, tau));
2649 }
2650 }
2651 arma::uword cvIdx = arma::index_min(mse);
2652 betaHat = smqrMcpPara(Z, Y, lambdaSeq(cvIdx), sx1, tau, p, 1.0 / n, h, h1, h3, phi0, gamma, epsilon, iteMax, iteTight, para);
2653 betaHat.rows(1, p) %= sx1;
2654 betaHat(0) += my - arma::as_scalar(mx * betaHat.rows(1, p));
2655 return Rcpp::List::create(Rcpp::Named("coeff") = betaHat, Rcpp::Named("lambda") = lambdaSeq(cvIdx));
2656 }
2657
2658 // [[Rcpp::export]]
cvSmqrLassoTrian(const arma::mat & X,arma::vec Y,const arma::vec & lambdaSeq,const arma::vec & folds,const double tau,const int kfolds,const double h,const double phi0=0.01,const double gamma=1.2,const double epsilon=0.001,const int iteMax=500)2659 Rcpp::List cvSmqrLassoTrian(const arma::mat& X, arma::vec Y, const arma::vec& lambdaSeq, const arma::vec& folds, const double tau, const int kfolds,
2660 const double h, const double phi0 = 0.01, const double gamma = 1.2, const double epsilon = 0.001, const int iteMax = 500) {
2661 const int n = X.n_rows, p = X.n_cols, nlambda = lambdaSeq.size();
2662 const double h1 = 1.0 / h, h2 = 1.0 / (h * h);
2663 arma::vec betaHat(p + 1);
2664 arma::vec mse = arma::zeros(nlambda);
2665 arma::rowvec mx = arma::mean(X, 0);
2666 arma::vec sx1 = 1.0 / arma::stddev(X, 0, 0).t();
2667 arma::mat Z = arma::join_rows(arma::ones(n), standardize(X, mx, sx1, p));
2668 double my = arma::mean(Y);
2669 Y -= my;
2670 for (int j = 1; j <= kfolds; j++) {
2671 arma::uvec idx = arma::find(folds == j);
2672 arma::uvec idxComp = arma::find(folds != j);
2673 double n1Train = 1.0 / idxComp.size();
2674 arma::mat trainZ = Z.rows(idxComp), testZ = Z.rows(idx);
2675 arma::vec trainY = Y.rows(idxComp), testY = Y.rows(idx);
2676 for (int i = 0; i < nlambda; i++) {
2677 betaHat = smqrLassoTrian(trainZ, trainY, lambdaSeq(i), sx1, tau, p, n1Train, h, h1, h2, phi0, gamma, epsilon, iteMax);
2678 mse(i) += arma::accu(lossQr(testZ, testY, betaHat, tau));
2679 }
2680 }
2681 arma::uword cvIdx = arma::index_min(mse);
2682 betaHat = smqrLassoTrian(Z, Y, lambdaSeq(cvIdx), sx1, tau, p, 1.0 / n, h, h1, h2, phi0, gamma, epsilon, iteMax);
2683 betaHat.rows(1, p) %= sx1;
2684 betaHat(0) += my - arma::as_scalar(mx * betaHat.rows(1, p));
2685 return Rcpp::List::create(Rcpp::Named("coeff") = betaHat, Rcpp::Named("lambda") = lambdaSeq(cvIdx));
2686 }
2687
2688 // [[Rcpp::export]]
cvSmqrScadTrian(const arma::mat & X,arma::vec Y,const arma::vec & lambdaSeq,const arma::vec & folds,const double tau,const int kfolds,const double h,const double phi0=0.01,const double gamma=1.2,const double epsilon=0.001,const int iteMax=500,const int iteTight=3,const double para=3.7)2689 Rcpp::List cvSmqrScadTrian(const arma::mat& X, arma::vec Y, const arma::vec& lambdaSeq, const arma::vec& folds, const double tau, const int kfolds,
2690 const double h, const double phi0 = 0.01, const double gamma = 1.2, const double epsilon = 0.001, const int iteMax = 500,
2691 const int iteTight = 3, const double para = 3.7) {
2692 const int n = X.n_rows, p = X.n_cols, nlambda = lambdaSeq.size();
2693 const double h1 = 1.0 / h, h2 = 1.0 / (h * h);
2694 arma::vec betaHat(p + 1);
2695 arma::vec mse = arma::zeros(nlambda);
2696 arma::rowvec mx = arma::mean(X, 0);
2697 arma::vec sx1 = 1.0 / arma::stddev(X, 0, 0).t();
2698 arma::mat Z = arma::join_rows(arma::ones(n), standardize(X, mx, sx1, p));
2699 double my = arma::mean(Y);
2700 Y -= my;
2701 for (int j = 1; j <= kfolds; j++) {
2702 arma::uvec idx = arma::find(folds == j);
2703 arma::uvec idxComp = arma::find(folds != j);
2704 double n1Train = 1.0 / idxComp.size();
2705 arma::mat trainZ = Z.rows(idxComp), testZ = Z.rows(idx);
2706 arma::vec trainY = Y.rows(idxComp), testY = Y.rows(idx);
2707 for (int i = 0; i < nlambda; i++) {
2708 betaHat = smqrScadTrian(trainZ, trainY, lambdaSeq(i), sx1, tau, p, n1Train, h, h1, h2, phi0, gamma, epsilon, iteMax, iteTight, para);
2709 mse(i) += arma::accu(lossQr(testZ, testY, betaHat, tau));
2710 }
2711 }
2712 arma::uword cvIdx = arma::index_min(mse);
2713 betaHat = smqrScadTrian(Z, Y, lambdaSeq(cvIdx), sx1, tau, p, 1.0 / n, h, h1, h2, phi0, gamma, epsilon, iteMax, iteTight, para);
2714 betaHat.rows(1, p) %= sx1;
2715 betaHat(0) += my - arma::as_scalar(mx * betaHat.rows(1, p));
2716 return Rcpp::List::create(Rcpp::Named("coeff") = betaHat, Rcpp::Named("lambda") = lambdaSeq(cvIdx));
2717 }
2718
2719 // [[Rcpp::export]]
cvSmqrMcpTrian(const arma::mat & X,arma::vec Y,const arma::vec & lambdaSeq,const arma::vec & folds,const double tau,const int kfolds,const double h,const double phi0=0.01,const double gamma=1.5,const double epsilon=0.001,const int iteMax=500,const int iteTight=3,const double para=3)2720 Rcpp::List cvSmqrMcpTrian(const arma::mat& X, arma::vec Y, const arma::vec& lambdaSeq, const arma::vec& folds, const double tau, const int kfolds,
2721 const double h, const double phi0 = 0.01, const double gamma = 1.5, const double epsilon = 0.001, const int iteMax = 500,
2722 const int iteTight = 3, const double para = 3) {
2723 const int n = X.n_rows, p = X.n_cols, nlambda = lambdaSeq.size();
2724 const double h1 = 1.0 / h, h2 = 1.0 / (h * h);
2725 arma::vec betaHat(p + 1);
2726 arma::vec mse = arma::zeros(nlambda);
2727 arma::rowvec mx = arma::mean(X, 0);
2728 arma::vec sx1 = 1.0 / arma::stddev(X, 0, 0).t();
2729 arma::mat Z = arma::join_rows(arma::ones(n), standardize(X, mx, sx1, p));
2730 double my = arma::mean(Y);
2731 Y -= my;
2732 for (int j = 1; j <= kfolds; j++) {
2733 arma::uvec idx = arma::find(folds == j);
2734 arma::uvec idxComp = arma::find(folds != j);
2735 double n1Train = 1.0 / idxComp.size();
2736 arma::mat trainZ = Z.rows(idxComp), testZ = Z.rows(idx);
2737 arma::vec trainY = Y.rows(idxComp), testY = Y.rows(idx);
2738 for (int i = 0; i < nlambda; i++) {
2739 betaHat = smqrMcpTrian(trainZ, trainY, lambdaSeq(i), sx1, tau, p, n1Train, h, h1, h2, phi0, gamma, epsilon, iteMax, iteTight, para);
2740 mse(i) += arma::accu(lossQr(testZ, testY, betaHat, tau));
2741 }
2742 }
2743 arma::uword cvIdx = arma::index_min(mse);
2744 betaHat = smqrMcpTrian(Z, Y, lambdaSeq(cvIdx), sx1, tau, p, 1.0 / n, h, h1, h2, phi0, gamma, epsilon, iteMax, iteTight, para);
2745 betaHat.rows(1, p) %= sx1;
2746 betaHat(0) += my - arma::as_scalar(mx * betaHat.rows(1, p));
2747 return Rcpp::List::create(Rcpp::Named("coeff") = betaHat, Rcpp::Named("lambda") = lambdaSeq(cvIdx));
2748 }
2749
2750
2751