1 /* -*- mode: c; kept-new-versions: 40; kept-old-versions: 40 -*-
2 * Indentation (etc) style: C-c . gnu */
3
4 /* file lmrob.c
5 * was roblm/src/roblm.c - version 0.6 by Matias Salibian-Barreras
6
7 * Includes the stable correct asymptotic variance estimators
8 * of Croux, Dhaene, Hoorelbeke
9 * Includes the fast-s algorithm
10 */
11
12 /* Robust MM regression estimates *
13 * ------------------------------ */
14
15 /* comment code *
16 *
17 * adapt other sampler <<<<<<<<<< R's random number generator !!!!
18
19 * replace abort for too many singular resamples by
20 * returning the number of singular ones
21 */
22
23 /* MM:
24 - Done: fast_s[_large_n]() both had FIXED seed (= 37),
25 and effectively discarded the seed_rand argument below
26 - Done: drop 'register' : today's compilers do optimize well!
27 - Done: using Calloc() / Free() instead of malloc()/free()
28 */
29
30 /* kollerma:
31 Added alternative psi functions callable via psifun, chifun and
32 wgtfun. ipsi is used to distinguish between the different types.
33 The ipsi argument works for the S-estimator as well as for the
34 MM-estimator.
35
36 - Added implementation of M-S algorithm.
37
38 - Modified subsampling behaviour: avoiding singular resamples by using
39 customized LU decomposition.
40
41 - Replaced C style matrices with Fortran style matrices, with as little
42 copying as possible.
43
44 - Using LAPACK's DGELS instead of local lu() decomposition.
45
46 - Code clean up: removed all subroutines that were unused.
47 */
48
49 #ifndef USE_FC_LEN_T
50 # define USE_FC_LEN_T
51 #endif
52 #include <Rconfig.h>
53 #include <R_ext/BLAS.h>
54 #include <R_ext/Lapack.h>
55 #ifndef FCONE
56 # define FCONE
57 #endif
58
59 #include <Rmath.h>
60 #include <R_ext/Applic.h>
61
62
63 #include "robustbase.h"
64 //-> <R.h>, <Rinternals.h> -> XLENGTH, R_xlen_t
65
66
67 /* these will also move to "lmrob.h" ---
68 * but first make many of these 'static' <<< FIXME!
69 */
70 void fast_s_large_n(double *X, double *y,
71 int *nn, int *pp, int *nRes, int *max_it_scale, double *res,
72 int *ggroups, int *nn_group,
73 int *K, int *max_k, double rel_tol, double inv_tol, double scale_tol, int *converged,
74 int *best_r, double *bb, double *rrhoc, int *iipsi,
75 double *bbeta, double *sscale, int trace_lev, int mts, Rboolean ss);
76
77 void fast_s(double *X, double *y,
78 int *nn, int *pp, int *nRes, int *max_it_scale, double *res,
79 int *K, int *max_k, double rel_tol, double inv_tol, double scale_tol, int *converged,
80 int *best_r, double *bb, double *rrhoc, int *iipsi,
81 double *bbeta, double *sscale, int trace_lev, int mts, Rboolean ss);
82
83 Rboolean rwls(const double X[], const double y[], int n, int p,
84 double *estimate, double *i_estimate,
85 double *resid, double *loss,
86 double scale, double epsilon,
87 int *max_it, double *rho_c, const int ipsi, int trace_lev);
88
89 static void sample_noreplace(int *x, int n, int k, int *ind_space);
90
91
92 double norm2 (double *x, int n);
93 double norm (double *x, int n);
94 double norm1(double *x, int n);
95 double norm_diff2 (double *x, double *y, int n);
96 double norm_diff (double *x, double *y, int n);
97 double norm1_diff(double *x, double *y, int n);
98
99 /* moved to robustbase.h
100 *
101 * double normcnst(const double c[], int ipsi);
102 * double rho_inf (const double c[], int ipsi);
103
104 * double rho(double x, const double c[], int ipsi);
105 * double psi(double x, const double c[], int ipsi);
106 * double psip(double x, const double c[], int ipsi);// psi'
107 * double psi2(double x, const double c[], int ipsi);// psi''
108 * double wgt(double x, const double c[], int ipsi);
109 */
110
111 double rho_huber(double x, const double c[]);
112 double psi_huber(double x, const double c[]);
113 double psip_huber(double x, const double c[]);
114 double psi2_huber(double x, const double c[]);
115 double wgt_huber(double x, const double c[]);
116
117 double rho_biwgt(double x, const double c[]);
118 double psi_biwgt(double x, const double c[]);
119 double psip_biwgt(double x, const double c[]);
120 double psi2_biwgt(double x, const double c[]);
121 double wgt_biwgt(double x, const double c[]);
122
123 double rho_gwgt(double x, const double c[]);
124 double psi_gwgt(double x, const double c[]);
125 double psip_gwgt(double x, const double c[]);
126 double wgt_gwgt(double x, const double c[]);
127
128 double rho_opt(double x, const double c[]);
129 double psi_opt(double x, const double c[]);
130 double psip_opt(double x, const double c[]);
131 double wgt_opt(double x, const double c[]);
132
133 double rho_hmpl(double x, const double c[]);
134 double psi_hmpl(double x, const double c[]);
135 double psip_hmpl(double x, const double c[]);
136 double psi2_hmpl(double x, const double c[]);
137 double wgt_hmpl(double x, const double c[]);
138
139 double rho_ggw(double x, const double c[]);
140 void psi_ggw_vec(double *x, int n, void *k);
141 double psi_ggw(double x, const double c[]);
142 double psip_ggw(double x, const double c[]);
143 double wgt_ggw(double x, const double c[]);
144
145 double rho_lqq(double x, const double c[]);
146 double psi_lqq(double x, const double c[]);
147 double psip_lqq(double x, const double c[]);
148 double psi2_lqq(double x, const double c[]);
149 double wgt_lqq(double x, const double c[]);
150
151 double sum_rho_sc(const double r[], double scale, int n, int p,
152 const double c[], int ipsi);
153 void get_weights_rhop(const double r[], double s, int n, const double rrhoc[], int ipsi,
154 /* --> */ double *w);
155 int refine_fast_s(const double X[], double *wx, const double y[], double *wy,
156 double *weights, int n, int p, double *res,
157 double *work, int lwork,
158 double *beta_cand,
159 int kk, Rboolean *conv, int max_k, double rel_tol,
160 int trace_lev,
161 double b, double *rrhoc, int ipsi, double initial_scale,
162 double *beta_ref, double *scale);
163
164 void m_s_subsample(double *X1, double *y, int n, int p1, int p2,
165 int nResample, int max_it_scale,
166 double rel_tol, double inv_tol, double scale_tol, double *bb,
167 double *rrhoc, int ipsi, double *sscale, int trace_lev,
168 double *b1, double *b2, double *t1, double *t2,
169 double *y_tilde, double *res, double *x1, double *x2,
170 int *NIT, int *K, int *KODE, double *SIGMA, double *BET0,
171 double *SC1, double *SC2, double *SC3, double *SC4, int mts, Rboolean ss);
172
173 Rboolean m_s_descent(double *X1, double *X2, double *y,
174 int n, int p1, int p2, int K_m_s, int max_k, int max_it_scale,
175 double rel_tol, double scale_tol, double *bb, double *rrhoc, int ipsi,
176 double *sscale, int trace_lev,
177 double *b1, double *b2, double *t1, double *t2,
178 double *y_tilde, double *res, double *res2, double *x1, double *x2,
179 int *NIT, int *K, int *KODE, double *SIGMA, double *BET0,
180 double *SC1, double *SC2, double *SC3, double *SC4);
181
182 int subsample(const double x[], const double y[], int n, int m,
183 double *beta, int *ind_space, int *idc, int *idr,
184 double *lu, double *v, int *p,
185 double *Dr, double *Dc, int rowequ, int colequ,
186 Rboolean sample, int mts, Rboolean ss, double tol_inv, Rboolean solve);
187
188 Rboolean fast_s_with_memory(double *X, double *y, double *res,
189 int *nn, int *pp, int *nRes, int *max_it_scale, int *K, int *max_k,
190 double rel_tol, double inv_tol, double scale_tol, int trace_lev,
191 int *best_r, double *bb, double *rrhoc, int *iipsi,
192 double **best_betas, double *best_scales, int mts, Rboolean ss);
193
194 /* for "tracing" only : */
195 void disp_mat(double **a, int n, int m);
196 void disp_vec(double *a, int n);
197 void disp_veci(int *a, int n);
198
199 double kthplace(double *, int, int);
200
201 int find_max(double *a, int n);
202
203 double find_scale(const double r[], double b, const double rrhoc[], int ipsi,
204 double initial_scale, int n, int p,
205 int* iter, // input: max_iter, output: #{iterations used}
206 double scale_tol,
207 Rboolean trace);
208
209 double median_abs(double *, int, double *);
210 double MAD(double *a, int n, double center, double *tmp, double *tmp2);
211
212 void zero_mat(double **a, int n, int m);
213
214
215 #define INIT_WLS(_X_, _y_, _n_, _p_) \
216 /* Determine optimal block size for work array*/ \
217 F77_CALL(dgels)("N", &_n_, &_p_, &one, _X_, &_n_, _y_, \
218 &_n_, &work0, &lwork, &info FCONE); \
219 if (info) { \
220 warning(_(" Problem determining optimal block size, using minimum")); \
221 lwork = 2*_p_; \
222 } else \
223 lwork = (int)work0; \
224 \
225 if (trace_lev >= 4) \
226 Rprintf(" Optimal block size for DGELS: %d\n", lwork); \
227 \
228 /* allocate */ \
229 work = (double *) Calloc(lwork, double); \
230 weights = (double *) Calloc(n, double);
231
232 #define CLEANUP_WLS \
233 Free(work); Free(weights);
234
235 #define CLEANUP_EQUILIBRATION \
236 Free(Dr); Free(Dc); Free(Xe);
237
238 #define CLEANUP_SUBSAMPLE \
239 Free(ind_space); Free(idc); Free(idr); Free(pivot); \
240 Free(lu); Free(v); \
241 CLEANUP_EQUILIBRATION;
242
243 #define FIT_WLS(_X_, _x_, _y_, _n_, _p_) \
244 /* add weights to _y_ and _x_ */ \
245 for (j=0; j<_n_; j++) { \
246 wtmp = sqrt(weights[j]); \
247 _y_[j] *= wtmp; \
248 for (k=0; k<_p_; k++) \
249 _x_[_n_*k+j] = _X_[_n_*k+j] * wtmp; \
250 } \
251 /* solve weighted least squares problem */ \
252 F77_CALL(dgels)("N", &_n_, &_p_, &one, _x_, &_n_, _y_, \
253 &_n_, work, &lwork, &info FCONE); \
254 if (info) { \
255 if (info < 0) { \
256 CLEANUP_WLS; \
257 error(_("DGELS: illegal argument in %i. argument."), info); \
258 } else { \
259 if (trace_lev >= 4) { \
260 Rprintf(" Robustness weights in failing step: "); \
261 disp_vec(weights, _n_); \
262 } \
263 CLEANUP_WLS; \
264 error(_("DGELS: weighted design matrix not of full rank (column %d).\nUse control parameter 'trace.lev = 4' to get diagnostic output."), info); \
265 } \
266 }
267
268 #define SETUP_EQUILIBRATION(_n_, _p_, _X_, _large_n_) \
269 /* equilibration of matrix _X_ */ \
270 /* solve (Dr X Dc) b = Dr y with beta = Dc b instead of */ \
271 /* X beta = y */ \
272 /* see Demmel (1997) APPLIED NUMERICAL LINEAR ALGEBRA */ \
273 /* Section 2.5.2 Equilibration */ \
274 double *Dr, *Dc, *Xe, rowcnd, colcnd, amax; \
275 int rowequ = 0 , colequ = 0; \
276 Dr = (double *) Calloc(_n_, double); \
277 Dc = (double *) Calloc(_p_, double); \
278 Xe = (double *) Calloc(_n_*_p_, double); \
279 COPY(_X_, Xe, _n_*_p_); \
280 F77_CALL(dgeequ)(&_n_, &_p_, Xe, &_n_, Dr, Dc, &rowcnd, \
281 &colcnd, &amax, &info); \
282 if (info) { \
283 if (info < 0) { \
284 CLEANUP_EQUILIBRATION; \
285 error(_("DGEEQ: illegal argument in %i. argument"), -1 * info); \
286 } else if (info > _n_) { \
287 if (_large_n_) { \
288 error(_("Fast S large n strategy failed. Use control parameter 'fast.s.large.n = Inf'.")); \
289 } else { \
290 error(_("DGEEQU: column %i of the design matrix is exactly zero."), info - _n_); \
291 } \
292 } else { \
293 /* FIXME: replace dgeequ by our own version */ \
294 /* that does not treat this as error */ \
295 warning(_(" Skipping design matrix equilibration (DGEEQU): row %i is exactly zero."), info); \
296 } \
297 } else { \
298 /* scale _X_ */ \
299 char equed; \
300 F77_CALL(dlaqge)(&_n_, &_p_, Xe, &_n_, Dr, Dc, &rowcnd, \
301 &colcnd, &amax, &equed FCONE); \
302 rowequ = equed == 'B' || equed == 'R'; \
303 colequ = equed == 'B' || equed == 'C'; \
304 }
305
306 #define SETUP_SUBSAMPLE(_n_, _p_, _X_, _large_n_) \
307 /* (Pointers to) Arrays - to be allocated */ \
308 int *ind_space, *idc, *idr, *pivot; \
309 double *lu, *v; \
310 ind_space = (int *) Calloc(_n_, int); \
311 idc = (int *) Calloc(_n_, int); \
312 idr = (int *) Calloc(_p_, int); \
313 pivot = (int *) Calloc(_p_-1, int); \
314 lu = (double *) Calloc(_p_*_p_, double); \
315 v = (double *) Calloc(_p_, double); \
316 SETUP_EQUILIBRATION(_n_, _p_, _X_, _large_n_);
317
318 #define COPY(from, to, len) Memcpy(to, from, len)
319 /* This assumes that 'p' is correctly defined, and 'j' can be used in caller: */
320 /* #define COPY(BETA_FROM, BETA_TO, _p_) \ */
321 /* for(j=0; j < _p_; j++) BETA_TO[j] = BETA_FROM[j]; */
322 /* In theory BLAS should be fast, but this seems slightly slower,
323 * particularly for non-optimized BLAS :*/
324 /* static int one = 1; */
325 /* #define COPY(BETA_FROM, BETA_TO, _p_) \ */
326 /* F77_CALL(dcopy)(&_p_, BETA_FROM, &one, BETA_TO, &one); */
327
328 #define EPS_SCALE 1e-10
329 #define INFI 1e+20
330
331 /* Called from R's lmrob.S() in ../R/lmrob.MM.R,
332 * help() in ../man/lmrob.S.Rd, this function computes an S-regression estimator
333 ~~~~~~~~~~~~~~~~~
334 */
R_lmrob_S(double * X,double * y,int * n,int * P,int * nRes,double * scale,double * beta_s,double * rrhoc,int * iipsi,double * bb,int * best_r,int * Groups,int * N_group,int * K_s,int * max_k,int * max_it_scale,double * rel_tol,double * inv_tol,double * scale_tol,int * converged,int * trace_lev,int * mts,int * ss,int * cutoff)335 void R_lmrob_S(double *X, double *y, int *n, int *P,
336 int *nRes, // = nResample ( = 500, by default)
337 double *scale, double *beta_s,
338 double *rrhoc, int *iipsi, double *bb,
339 int *best_r, int *Groups, int *N_group,
340 int *K_s, int *max_k, int *max_it_scale,
341 double *rel_tol, double *inv_tol,
342 double *scale_tol, // <- new, was hardwired to EPS_SCALE := 1e-10
343 int *converged, int *trace_lev, int *mts, int *ss, int *cutoff)
344 {
345 /* best_r = 't' of Salibian-Barrera_Yohai(2006),
346 * = no. of best candidates to be iterated further ("refined")
347 * = 2, by default
348 */
349 if (*nRes > 0) {
350 double *res = (double *) R_alloc(*n, sizeof(double)); // residuals
351 if (*n > *cutoff) {
352 if(*trace_lev > 0)
353 Rprintf("lmrob_S(n = %d, nRes = %d): fast_s_large_n():\n", *n, *nRes);
354 fast_s_large_n(X, y, n, P, nRes, max_it_scale, res,
355 Groups, N_group,
356 K_s, max_k, *rel_tol, *inv_tol, *scale_tol, converged,
357 best_r, bb, rrhoc, iipsi, beta_s, scale, *trace_lev, *mts, (Rboolean)*ss);
358 } else {
359 if(*trace_lev > 0)
360 Rprintf("lmrob_S(n = %d, nRes = %d): fast_s() [non-large n]:\n", *n, *nRes);
361 fast_s(X, y, n, P, nRes, max_it_scale, res,
362 K_s, max_k, *rel_tol, *inv_tol, *scale_tol, converged,
363 best_r, bb, rrhoc, iipsi, beta_s, scale, *trace_lev, *mts, *ss);
364 }
365 COPY(res, y, *n); // return the 'residuals' in 'y'
366 } else { // nRes[] <= 0 <==> 'only.scale = TRUE'
367 if(*trace_lev > 0)
368 Rprintf("lmrob_S(nRes = 0, n = %d): --> find_scale(*, scale=%g) only:",
369 *n, *scale);
370 *scale = find_scale(y, *bb, rrhoc, *iipsi, *scale, *n, *P,
371 max_it_scale, *scale_tol, *trace_lev >= 3);
372 if(*trace_lev > 0) Rprintf(" used %d iterations\n", *max_it_scale);
373 }
374 }
375
376 /* Called from R, this function computes an M-S-regression estimator */
377 // not only called from ../R/lmrob.M.S.R, but also ../inst/xtraR/m-s_fns.R
378 // ~~~~~~~~~~~~~~~~ ~~~~~~~~~~~~~~~~~~~~~~~
R_lmrob_M_S(double * X1,double * X2,double * y,double * res,int * nn,int * pp1,int * pp2,int * nRes,int * max_it_scale,double * scale,double * b1,double * b2,double * rho_c,int * ipsi,double * bb,int * K_m_s,int * max_k,double * rel_tol,double * inv_tol,double * scale_tol,int * converged,int * trace_lev,int * orthogonalize,int * subsample,int * descent,int * mts,int * ss)379 void R_lmrob_M_S(double *X1, double *X2, double *y, double *res,
380 int *nn, int *pp1, int *pp2, int *nRes, int *max_it_scale,
381 double *scale, double *b1, double *b2,
382 double *rho_c, int *ipsi, double *bb,
383 int *K_m_s, int *max_k,
384 double *rel_tol, double *inv_tol, double *scale_tol,
385 int *converged,
386 int *trace_lev,
387 int *orthogonalize, int *subsample, int *descent,
388 int *mts, int *ss)
389 {
390 /* Initialize (some of the) memory here,
391 * so that we have to do it only once */
392 int i, n = *nn, p1 = *pp1, p2 = *pp2, one = 1;
393
394 /* (Pointers to) Arrays - to be allocated */
395 double *t1, *t2, *y_tilde, *y_work, done = 1., dmone = -1.;
396 double *x1, *x2, *ot1, *oT2, *ptr;
397
398 if(*trace_lev > 0) Rprintf(
399 "lmrob_M_S(n = %d, nRes = %d, (p1,p2)=(%d,%d), (orth,subs,desc)=(%d,%d,%d))\n",
400 n, *nRes, p1, p2,
401 *orthogonalize, *subsample, *descent);
402
403 t1 = (double *) R_alloc(n, sizeof(double)); /* size n needed for rllarsbi */
404 t2 = (double *) R_alloc(p2, sizeof(double));
405 ot1 = (double *) R_alloc(p1, sizeof(double));
406 oT2 = (double *) R_alloc(p2*p1, sizeof(double));
407 y_work = (double *) R_alloc(n, sizeof(double));
408 COPY(y, y_work, n);
409 y_tilde = (double *) R_alloc(n, sizeof(double));
410 x1 = (double *) R_alloc(n*p1, sizeof(double));
411 x2 = (double *) R_alloc(n*p2, sizeof(double));
412 COPY(X2, x2, n*p2);
413
414 /* Variables required for rllarsbi
415 * (l1 / least absolut residuals - estimate) */
416 int NIT=0, K=0, KODE=0;
417 double SIGMA = 0.,
418 *SC1 = (double *) R_alloc(n, sizeof(double)),
419 *SC2 = (double *) R_alloc(p1, sizeof(double)),
420 *SC3 = (double *) R_alloc(p1, sizeof(double)),
421 *SC4 = (double *) R_alloc(p1, sizeof(double));
422 double BET0 = 0.773372647623; /* = pnorm(0.75) */
423
424 /* STEP 1: Orthgonalize X2 and y from X1 */
425 if (*orthogonalize) {
426 COPY(X1, x1, n*p1);
427 F77_CALL(rllarsbi)(x1, y_work, &n, &p1, &n, &n, rel_tol,
428 &NIT, &K, &KODE, &SIGMA, t1, y_tilde, SC1, SC2,
429 SC3, SC4, &BET0);
430 COPY(t1, ot1, p1);
431 for (i=0; i < p2; i++) {
432 COPY(X1, x1, n*p1);
433 ptr = X2+i*n; COPY(ptr, y_work, n);
434 F77_CALL(rllarsbi)(x1, y_work, &n, &p1, &n, &n, rel_tol,
435 &NIT, &K, &KODE, &SIGMA, t1, x2+i*n, SC1, SC2,
436 SC3, SC4, &BET0);
437 ptr = oT2+i*p1; COPY(t1, ptr, p1);
438 }
439 COPY(y_tilde, y_work, n);
440 /* compare with Maronna & Yohai 2000:
441 * y_work and y_tilde now contain \tilde y, ot1 -> t_1,
442 * x2 -> \tilde x2, oT2 -> T_2 */
443 }
444
445 /* STEP 2: Subsample */
446 if (*subsample) {
447 m_s_subsample(X1, y_work, n, p1, p2, *nRes, *max_it_scale,
448 *rel_tol, *inv_tol, *scale_tol, bb,
449 rho_c, *ipsi, scale, *trace_lev,
450 b1, b2, t1, t2, y_tilde, res, x1, x2,
451 &NIT, &K, &KODE, &SIGMA, &BET0,
452 SC1, SC2, SC3, SC4, *mts, *ss);
453
454 if (*scale < 0)
455 error(_("m_s_subsample() stopped prematurely (scale < 0)."));
456 }
457
458 /* STEP 3: Transform back */
459 if (*orthogonalize) {
460 /* t1 = ot1 + b1 - oT2 %*% b2 */
461 for(int i=0; i < p1; i++) t1[i] = ot1[i] + b1[i];
462 F77_CALL(dgemv)("N", &p1, &p2, &dmone, oT2, &p1, b2, &one, &done, t1, &one FCONE);
463 COPY(t1, b1, p1);
464 /* restore x2 */
465 COPY(X2, x2, n*p2);
466 }
467
468 /* update / calculate residuals */
469 COPY(y, res, n);
470 F77_CALL(dgemv)("N", &n, &p1, &dmone, X1, &n, b1, &one, &done, res, &one FCONE);
471 F77_CALL(dgemv)("N", &n, &p2, &dmone, X2, &n, b2, &one, &done, res, &one FCONE);
472
473 /* STEP 4: Descent procedure */
474 if (*descent) {
475 *converged = m_s_descent(
476 X1, X2, y, n, p1, p2, *K_m_s, *max_k, *max_it_scale,
477 *rel_tol, *scale_tol, bb, rho_c, *ipsi, scale, *trace_lev,
478 b1, b2, t1, t2, y_tilde, res, y_work, x1, x2,
479 &NIT, &K, &KODE, &SIGMA, &BET0, SC1, SC2, SC3, SC4);
480 }
481 }
482
483 /* This function performs RWLS iterations starting from
484 * an S-regression estimator (and associated residual scale).
485 * So, in itself, this is ``just'' an M-estimator -- called from R's
486 * lmrob..M..fit() [ ../R/lmrob.MM.R ]
487 * ~~~~~~~~~~~~~~~
488 * NOTE: rel_tol now controls the *relative* changes in beta,
489 * instead of being hard-wired to EPS = 1e-7 and bounding the
490 * absolute || beta_1 - beta_2 ||
491 */
R_lmrob_MM(double * X,double * y,int * n,int * P,double * beta_initial,double * scale,double * beta_m,double * resid,int * max_it,double * rho_c,int * ipsi,double * loss,double * rel_tol,int * converged,int * trace_lev,int * mts,int * ss)492 void R_lmrob_MM(double *X, double *y, int *n, int *P,
493 double *beta_initial, double *scale,
494 double *beta_m, double *resid,
495 int *max_it, double *rho_c, int *ipsi, double *loss,
496 double *rel_tol, int *converged, int *trace_lev, int *mts, int *ss)
497 {
498 /* starting from the S-estimate (beta_initial), use
499 * irwls to compute the MM-estimate (beta_m) */
500
501 if(*trace_lev > 0)
502 Rprintf("lmrob_MM(): rwls():\n");
503
504 *converged = (int)rwls(X,y,*n,*P,beta_m, beta_initial, resid, loss,
505 *scale, *rel_tol,
506 max_it, rho_c, *ipsi, *trace_lev);
507 if (!converged)
508 COPY(beta_initial, beta_m, *P);
509 }
510
511
512 /* Call subsample() from R, for testing purposes only */
R_subsample(const double x[],const double y[],int * n,int * m,double * beta,int * ind_space,int * idc,int * idr,double * lu,double * v,int * p,double * _Dr,double * _Dc,int * _rowequ,int * _colequ,int * status,int * sample,int * mts,int * ss,double * tol_inv,int * solve)513 void R_subsample(const double x[], const double y[], int *n, int *m,
514 double *beta, int *ind_space, int *idc, int *idr,
515 double *lu, double *v, int *p,
516 double *_Dr, double *_Dc, int *_rowequ, int *_colequ,
517 int *status, int *sample, int *mts, int *ss, double *tol_inv,
518 int *solve)
519 {
520 int info;
521
522 /* set the seed */
523 GetRNGstate();
524
525 SETUP_EQUILIBRATION(*n, *m, x, 0);
526
527 *status = subsample(Xe, y, *n, *m, beta, ind_space, idc, idr, lu, v, p,
528 Dr, Dc, rowequ, colequ, (Rboolean)*sample, *mts, (Rboolean)*ss,
529 *tol_inv, (Rboolean)*solve);
530
531 COPY(Dr, _Dr, *n);
532 COPY(Dc, _Dc, *m);
533 *_rowequ = rowequ;
534 *_colequ = colequ;
535
536 CLEANUP_EQUILIBRATION;
537
538 PutRNGstate();
539 }
540
541 //---- Psi(), Rho(), Functions-----------------------------------------------------------
542
R_psifun(SEXP x_,SEXP c_,SEXP ipsi_,SEXP deriv_)543 SEXP R_psifun(SEXP x_, SEXP c_, SEXP ipsi_, SEXP deriv_) {
544 /*
545 * Calculate psi for vectorized x, scaled to get psi'(0) = 1
546 * deriv -1: rho(x) {*not* normalized}
547 * deriv 0: psi(x) = rho'(x)
548 * deriv 1: psi'(x) = rho''(x) {we always have psip(0) == 1}
549 * deriv 2: psi''(x)= rho'''(x)
550 */
551 int nprot = 1, ipsi = asInteger(ipsi_), deriv = asInteger(deriv_);
552 if (isInteger(x_)) {
553 x_ = PROTECT(coerceVector(x_, REALSXP)); nprot++;
554 }
555 if (!isReal(x_)) error(_("Argument '%s' must be numeric or integer"), "x");
556 if (!isReal(c_)) error(_("Argument '%s' must be numeric or integer"), "cc");
557 R_xlen_t i, n = XLENGTH(x_);
558 SEXP res = PROTECT(allocVector(REALSXP, n)); // the result
559 double *x = REAL(x_), *r = REAL(res), *cc = REAL(c_);
560
561 // put the for() loop *inside* the switch (--> speed for llength >> 1) :
562 #define for_i_n_NA for(i = 0; i < n; i++) r[i] = ISNAN(x[i]) ? x[i] :
563
564 switch(deriv) { // our rho() is rho~(), i.e., scaled to max = 1
565 case -1:
566 if(is_redescender(ipsi)) {
567 double rho_Inf = rho_inf(cc, ipsi);
568 for_i_n_NA rho(x[i], cc, ipsi) * rho_Inf;
569 } else { // huber, ..
570 for_i_n_NA rho(x[i], cc, ipsi);
571 }
572 break;
573
574 case 0: for_i_n_NA psi (x[i], cc, ipsi); break;
575 case 1: for_i_n_NA psip(x[i], cc, ipsi); break;
576 case 2: for_i_n_NA psi2(x[i], cc, ipsi); break;
577 default: error(_("'deriv'=%d is invalid"), deriv);
578 }
579 UNPROTECT(nprot);
580 return res;
581 }
582
R_chifun(SEXP x_,SEXP c_,SEXP ipsi_,SEXP deriv_)583 SEXP R_chifun(SEXP x_, SEXP c_, SEXP ipsi_, SEXP deriv_) {
584 /*
585 * Calculate chi for vectorized x, i.e. rho~(.) with rho~(inf) = 1:
586 * deriv 0: chi (x) = \rho(x) / \rho(Inf) =: \rho(x) * nc == our rho() C-function
587 * deriv 1: chi'(x) = psi(x) * nc
588 * deriv 2: chi''(x) = psi'(x) * nc
589 */
590 int nprot = 1, ipsi = asInteger(ipsi_), deriv = asInteger(deriv_);
591 if (isInteger(x_)) {
592 x_ = PROTECT(coerceVector(x_, REALSXP)); nprot++;
593 }
594 if (!isReal(x_)) error(_("Argument '%s' must be numeric or integer"), "x");
595 if (!isReal(c_)) error(_("Argument '%s' must be numeric or integer"), "cc");
596 R_xlen_t i, n = XLENGTH(x_);
597 SEXP res = PROTECT(allocVector(REALSXP, n)); // the result
598 double *x = REAL(x_), *r = REAL(res), *cc = REAL(c_);
599
600 // our rho() is rho~() == chi(), i.e., scaled to max = 1
601 double rI = (deriv > 0) ? rho_inf(cc, ipsi) : 0./* -Wall */;
602
603 switch(deriv) {
604 case 0: for_i_n_NA rho(x[i], cc, ipsi); break;
605 case 1: for_i_n_NA psi (x[i], cc, ipsi) / rI; break;
606 case 2: for_i_n_NA psip(x[i], cc, ipsi) / rI; break;
607 case 3: for_i_n_NA psi2(x[i], cc, ipsi) / rI; break;
608 default: error(_("'deriv'=%d is invalid"), deriv);
609 }
610 UNPROTECT(nprot);
611 return res;
612 }
613
R_wgtfun(SEXP x_,SEXP c_,SEXP ipsi_)614 SEXP R_wgtfun(SEXP x_, SEXP c_, SEXP ipsi_) {
615 /*
616 * Calculate wgt(x) = psi(x)/x for vectorized x
617 */
618 int nprot = 1, ipsi = asInteger(ipsi_);
619 if (isInteger(x_)) {
620 x_ = PROTECT(coerceVector(x_, REALSXP)); nprot++;
621 }
622 if (!isReal(x_)) error(_("Argument '%s' must be numeric or integer"), "x");
623 if (!isReal(c_)) error(_("Argument '%s' must be numeric or integer"), "cc");
624 R_xlen_t i, n = XLENGTH(x_);
625 SEXP res = PROTECT(allocVector(REALSXP, n)); // the result
626 double *x = REAL(x_), *r = REAL(res), *cc = REAL(c_);
627
628 for_i_n_NA wgt(x[i], cc, ipsi);
629
630 UNPROTECT(nprot);
631 return res;
632 }
633
634 #undef for_i_n_NA
635
R_rho_inf(SEXP cc,SEXP ipsi)636 SEXP R_rho_inf(SEXP cc, SEXP ipsi) {
637 if (!isReal(cc)) error(_("Argument 'cc' must be numeric"));
638 if (!isInteger(ipsi)) error(_("Argument 'ipsi' must be integer"));
639 return ScalarReal(rho_inf(REAL(cc), INTEGER(ipsi)[0]));
640 }
641
rho_inf(const double k[],int ipsi)642 double rho_inf(const double k[], int ipsi) {
643 /*
644 * Compute \rho(\infty) for psi functions
645 * (Note that our C function rho() is "rho~" and has rho(Inf) = 1)
646 */
647 double c = k[0];
648
649 switch(ipsi) {
650 default: error(_("rho_inf(): ipsi=%d not implemented."), ipsi);
651 case 0: return(R_PosInf); // huber
652 case 1: return(c*c/6.); // biweight
653 case 2: return(c*c); // GaussWeight / "Welsh"
654 case 3: return(3.25*c*c); // Optimal
655 case 4: return(0.5*k[0]*(k[1]+k[2]-k[0])); // Hampel
656 case 5: // GGW (Generalized Gauss Weight)
657 switch((int)c) {
658 default:
659 case 0: return(k[4]); break; // k[4] == cc[5] in R -- must be correct!
660 case 1: return(5.309853); break;
661 case 2: return(2.804693); break;
662 case 3: return(0.3748076); break;
663 case 4: return(4.779906); break;
664 case 5: return(2.446574); break;
665 case 6: return(0.4007054); break;
666 };
667 case 6: // LQQ aka 'lin psip'
668 return (k[2]*k[1]*(3*k[1]+2*k[0]) + (k[0]+k[1])*(k[0]+k[1])) / (6.*(k[2]-1.));
669 }
670 } // rho_inf()
671
normcnst(const double k[],int ipsi)672 double normcnst(const double k[], int ipsi) {
673 /*
674 * return normalizing constant for psi functions := 1 / \rho(\infty)
675 */
676
677 double c = k[0];
678
679 switch(ipsi) {
680 default: error(_("normcnst(): ipsi=%d not implemented."), ipsi);
681 case 0: return(0.); // huber {normcnst() should never be used for that!}
682 case 1: return(6./(c*c)); // biweight
683 case 2: return(1./(c*c)); // GaussWeight / "Welsh"
684 case 3: return(1./3.25/(c*c)); // Optimal
685 case 4: return(2./(k[0]*(k[1]+k[2]-k[0]))); // Hampel
686 case 5: // GGW
687 switch((int)c) {
688 default:
689 case 0: return(1./ k[4]); break; // k[4] == cc[5] in R -- must be correct!
690 case 1: return(1./5.309853); break;
691 case 2: return(1./2.804693); break;
692 case 3: return(1./0.3748076); break;
693 case 4: return(1./4.779906); break;
694 case 5: return(1./2.446574); break;
695 case 6: return(1./0.4007054); break;
696 };
697 case 6: // LQQ aka 'lin psip'
698 return((6*(k[2]-1))/(k[2]*k[1]*(3*k[1]+2*k[0])+(k[0]+k[1])*(k[0]+k[1])));
699 }
700 } // normcnst()
701
702
rho(double x,const double c[],int ipsi)703 double rho(double x, const double c[], int ipsi)
704 {
705 /*
706 * return the correct rho according to ipsi
707 * This rho() is normalized to 1, called rho~() or chi() in other contexts
708 */
709 switch(ipsi) {
710 default: error(_("rho(): ipsi=%d not implemented."), ipsi);
711 case 0: return(rho_huber(x, c)); // huber
712 case 1: return(rho_biwgt(x, c)); // biweight
713 case 2: return(rho_gwgt(x, c)); // GaussWeight / "Welsh"
714 case 3: return(rho_opt(x, c)); // Optimal
715 case 4: return(rho_hmpl(x, c)); // Hampel
716 case 5: return(rho_ggw(x, c)); // GGW (Generalized Gauss Weight)
717 case 6: return(rho_lqq(x, c)); // LQQ := Linear-Quadratic-Quadratic
718 // was LGW := "lin psip" := piecewise linear psi'()
719 }
720 }
721
psi(double x,const double c[],int ipsi)722 double psi(double x, const double c[], int ipsi)
723 {
724 /*
725 * return the correct psi according to ipsi
726 * this is actually rho' and not psi
727 */
728 switch(ipsi) {
729 default: error(_("psi(): ipsi=%d not implemented."), ipsi);
730 case 0: return(psi_huber(x, c)); // huber
731 case 1: return(psi_biwgt(x, c)); // biweight
732 case 2: return(psi_gwgt(x, c)); // GaussWeight / "Welsh"
733 case 3: return(psi_opt(x, c)); // Optimal
734 case 4: return(psi_hmpl(x, c)); // Hampel
735 case 5: return(psi_ggw(x, c)); // GGW
736 case 6: return(psi_lqq(x, c)); // LQQ (piecewise linear psi')
737 }
738 }
739
psip(double x,const double c[],int ipsi)740 double psip(double x, const double c[], int ipsi)
741 {
742 /*
743 * return the correct ppsi according to ipsi
744 * this is actually rho'' and not psip
745 */
746 switch(ipsi) {
747 default: error(_("psip(): ipsi=%d not implemented."), ipsi);
748 case 0: return(psip_huber(x, c)); // huber
749 case 1: return(psip_biwgt(x, c)); // biweight
750 case 2: return(psip_gwgt(x, c)); // GaussWeight / "Welsh"
751 case 3: return(psip_opt(x, c)); // Optimal
752 case 4: return(psip_hmpl(x, c)); // Hampel
753 case 5: return(psip_ggw(x, c)); // GGW
754 case 6: return(psip_lqq(x, c)); // LQQ (piecewise linear psi')
755 }
756 }
757
psi2(double x,const double c[],int ipsi)758 double psi2(double x, const double c[], int ipsi)
759 {
760 /* Compute psi''(x) == rho'''(x)
761 */
762 switch(ipsi) {
763 // default: error(_("psi2: ipsi=%d not implemented."), ipsi);
764 case 0: return(psi2_huber(x, c)); // huber
765 case 1: return(psi2_biwgt(x, c)); // biweight
766 case 4: return(psi2_hmpl(x, c)); // Hampel
767 case 6: return(psi2_lqq(x, c)); // LQQ (piecewise linear psi')
768
769 default: error(_("psi2(): ipsi=%d not implemented."), ipsi);
770 /*
771 case 2: return(psi2_gwgt(x, c)); // GaussWeight / "Welsh"
772 case 3: return(psi2_opt(x, c)); // Optimal
773 case 5: return(psi2_ggw(x, c)); // GGW
774 */
775 }
776 }
777
wgt(double x,const double c[],int ipsi)778 double wgt(double x, const double c[], int ipsi)
779 {
780 /*
781 * return the correct wgt according to ipsi
782 * wgt: rho'(x) / x
783 */
784 switch(ipsi) {
785 default:
786 case 0: return(wgt_huber(x, c)); // huber
787 case 1: return(wgt_biwgt(x, c)); // biweight
788 case 2: return(wgt_gwgt(x, c)); // GaussWeight / "Welsh"
789 case 3: return(wgt_opt(x, c)); // Optimal
790 case 4: return(wgt_hmpl(x, c)); // Hampel
791 case 5: return(wgt_ggw(x, c)); // GGW
792 case 6: return(wgt_lqq(x, c)); // LQQ (piecewise linear psi')
793 }
794 }
795
796 //--- Huber's rho / psi / ...
797 //--- -------
798
799 /* Huber's rho(): contrary to all the redescenders below,
800 this can NOT be scaled to rho(Inf)=1 : */
rho_huber(double x,const double c[])801 double rho_huber(double x, const double c[])
802 {
803 return (fabs(x) <= c[0]) ? x*x*0.5 : c[0]*(fabs(x) - c[0]/2);
804 }
805
psi_huber(double x,const double c[])806 double psi_huber(double x, const double c[])
807 {
808 // Huber's psi = rho'()
809 return (x <= -c[0]) ? -c[0] : ((x < c[0]) ? x : c[0]);
810 }
811
psip_huber(double x,const double c[])812 double psip_huber(double x, const double c[])
813 {
814 // psi' = rho'' : Second derivative of Huber's loss function
815 return (fabs(x) >= c[0]) ? 0. : 1.;
816 }
817
psi2_huber(double x,const double c[])818 double psi2_huber(double x, const double c[])
819 {
820 // psi'' = rho''' : Third derivative of Huber's loss function
821 return 0;
822 // FIXME? return NaN when |x| == c ?? -- then also for psi2_hmpl()
823 }
824
wgt_huber(double x,const double c[])825 double wgt_huber(double x, const double c[])
826 {
827 /*
828 * Weights for Huber's loss function w(x) = psi(x)/x
829 */
830 return (fabs(x) >= c[0]) ? c[0]/fabs(x) : 1.;
831 }
832
833 //--- Biweight = Bisquare = Tukey's Biweight ...
834 //--- --------------------------------------
835
rho_biwgt(double x,const double c[])836 double rho_biwgt(double x, const double c[])
837 {
838 /*
839 * Tukey's bisquare loss function == R's tukeyChi()
840 */
841 if (fabs(x) > (*c))
842 return(1.);
843 else {
844 double t = x / (*c);
845 t *= t; /* = t^2 */
846 return( t*(3. + t*(-3. + t)) );
847 }
848 }
849
psi_biwgt(double x,const double c[])850 double psi_biwgt(double x, const double c[])
851 {
852 /*
853 * First derivative of Tukey's bisquare loss function
854 */
855 if (fabs(x) > (*c))
856 return(0.);
857 else {
858 double a = x / (*c),
859 u = 1. - a*a;
860 return( x * u * u );
861 }
862 }
863
psip_biwgt(double x,const double c[])864 double psip_biwgt(double x, const double c[])
865 {
866 /*
867 * Second derivative of Tukey's bisquare loss function
868 */
869 if (fabs(x) > (*c))
870 return(0.);
871 else {
872 x /= *c;
873 double x2 = x*x;
874 return( (1. - x2) * (1 - 5 * x2));
875 }
876 }
877
psi2_biwgt(double x,const double c[])878 double psi2_biwgt(double x, const double c[])
879 {
880 /** 3rd derivative of Tukey's bisquare loss function rho()
881 *= 2nd derivative of psi() :
882 */
883 if (fabs(x) >= c[0]) // psi''() *is* discontinuous at x = c[0]: use "middle" value there:
884 return (fabs(x) == c[0]) ? 4*x/c[0] : 0.;
885 else {
886 x /= c[0];
887 double x2 = x*x;
888 return 4*x/c[0] * (5 * x2 - 3.);
889 }
890 }
891
wgt_biwgt(double x,const double c[])892 double wgt_biwgt(double x, const double c[])
893 {
894 /*
895 * Weights for Tukey's bisquare loss function
896 */
897 if( fabs(x) > *c )
898 return(0.);
899 else {
900 double a = x / (*c);
901 a = (1. - a)*(1. + a);
902 return( a * a );
903 }
904 }
905
906 //---------- gwgt == Gauss Weight Loss function =: "Welsh" --------------------
907
rho_gwgt(double x,const double c[])908 double rho_gwgt(double x, const double c[])
909 {
910 /*
911 * Gauss Weight Loss function
912 */
913 double ac = x / (*c);
914 return(-expm1(-(ac*ac)/2));
915 }
916
917 // Largest x such that exp(-x) does not underflow :
918 static double MIN_Exp = -708.4; // ~ = M_LN2 * DBL_MIN_EXP = -log(2) * 1022 = -708.3964 */
919
920 // Largest x such that exp(-x^2/2) does not underflow :
921 static double MAX_Ex2 = 37.7; // ~ = sqrt(- 2. * M_LN2 * DBL_MIN_EXP);
922 /* max {x | exp(-x^2/2) < .Machine$double.xmin } =
923 * min {x | x^2 > -2*log(2)* .Machine$double.min.exp } =
924 * = sqrt(-2*log(2)* .Machine$double.min.exp) = {IEEE double}
925 * = sqrt(log(2) * 2044) = 37.64031 */
926
psi_gwgt(double x,const double c[])927 double psi_gwgt(double x, const double c[])
928 {
929 /*
930 * Gauss Weight Psi()
931 */
932 double a = x / (*c);
933 if(fabs(a) > MAX_Ex2)
934 return 0.;
935 else
936 return x*exp(-(a*a)/2);
937 }
938
psip_gwgt(double x,const double c[])939 double psip_gwgt(double x, const double c[])
940 {
941 /*
942 * Gauss Weight Psi'()
943 */
944 x /= (*c);
945 if(fabs(x) > MAX_Ex2)
946 return 0.;
947 else {
948 double ac = x*x;
949 return exp(-ac/2) * (1. - ac);
950 }
951 }
952
wgt_gwgt(double x,const double c[])953 double wgt_gwgt(double x, const double c[])
954 {
955 /*
956 * Gauss Weight Loss function
957 */
958 double a = x / (*c);
959 return(exp(-(a*a)/2));
960 }
961
rho_opt(double x,const double c[])962 double rho_opt(double x, const double c[])
963 {
964 /*
965 * Optimal psi Function, thank you robust package
966 */
967 double ac = x / (*c), // AX=S/XK
968 ax = fabs(ac); // AX=ABST/XK
969 if (ax > 3) // IF (AX .GT. 3.D0) THEN
970 return(1); // rlRHOm2=3.25D0*XK*XK
971 else if (ax > 2.) {
972 const double R1 = -1.944/2., R2 = 1.728/4., R3 = -0.312/6., R4 = 0.016/8.;
973 ax *= ax; // = |x/c| ^ 2
974 return (ax*(R1+ ax*(R2+ ax*(R3+ ax*R4))) +1.792)/3.25;
975 // rlRHOm2=XK*XK*(R1*AX**2+R2*AX**4+R3*AX**6+R4*AX**8+1.792D0)
976 }
977 else
978 return(ac*ac/6.5); // rlRHOm2=S2/2.D0
979 }
980
psi_opt(double x,const double c[])981 double psi_opt(double x, const double c[])
982 {
983 /*
984 * Optimal psi Function, thank you robust package
985 */
986 double R1 = -1.944, R2 = 1.728, R3 = -0.312, R4 = 0.016;
987 double ax, ac;
988 ac = x / (*c); // AX=S/XK
989 ax = fabs(ac); // AX=ABST/XK
990 if (ax > 3.) // IF (AX .GT. 3.D0) THEN
991 return(0); // rlPSIm2=0.D0
992 else if (ax > 2.) { // ELSE IF(AX .GT. 2.D0) THEN
993 double a2 = ac*ac;
994 if (ac > 0.) // IF (AX .GT. 0.D0) THEN
995 return fmax2(0., (*c)*((((R4*a2 +R3)*a2 +R2)*a2 +R1)*ac));
996 // rlPSIm2=DMAX1(0.D0,XK*(R4*AX**7+R3*AX**5+R2*AX**3+R1*AX))
997 else
998 return -fabs((*c)*((((R4*a2 +R3)*a2 +R2)*a2 +R1)*ac));
999 // rlPSIm2=-DABS(XK*(R4*AX**7+R3*AX**5+R2*AX**3+R1*AX))
1000 }
1001 else
1002 return x;
1003 }
1004
psip_opt(double x,const double c[])1005 double psip_opt(double x, const double c[])
1006 {
1007 /*
1008 * psi'() for Optimal psi Function, thank you robust package
1009 */
1010 double ac = x / (*c),
1011 ax = fabs(ac);
1012 if (ax > 3.)
1013 return 0;
1014 else if (ax > 2.) {
1015 const double R1 = -1.944, R2 = 1.728, R3 = -0.312, R4 = 0.016;
1016 ax *= ax; // = |x/c| ^ 2
1017 return R1 + ax*(3*R2 + ax*(5*R3 + ax * 7*R4));
1018 } else
1019 return 1;
1020 }
1021
wgt_opt(double x,const double c[])1022 double wgt_opt(double x, const double c[])
1023 {
1024 /*
1025 * w(.) for optimal psi Function, thank you robust package
1026 */
1027 double ac = x / (*c),
1028 ax = fabs(ac);
1029 if (ax > 3.)
1030 return 0.;
1031 else if (ax > 2.) {
1032 const double R1 = -1.944, R2 = 1.728, R3 = -0.312, R4 = 0.016;
1033 ax *= ax; // = |x/c| ^ 2
1034 return fmax2(0., R1+ ax*(R2 + ax*(R3 + ax*R4)));
1035 }
1036 else
1037 return 1.;
1038 }
1039
1040
rho_hmpl(double x,const double k[])1041 double rho_hmpl(double x, const double k[])
1042 {
1043 /*
1044 * rho() for Hampel's redescending psi function
1045 * constants (a, b, r) == k[0:2] s.t. slope of psi is 1 in the center
1046 *
1047 * This function is normalized s.t. rho(inf) = 1
1048 */
1049 double u = fabs(x),
1050 nc = k[0]*(k[1]+k[2]-k[0])/2;
1051
1052 if (u <= k[0])
1053 return( x*x/2 / nc );
1054 else if (u <= k[1])
1055 return( ( u - k[0]/2 ) * k[0] / nc );
1056 else if (u <= k[2])
1057 return( ( k[1] - k[0]/2 + (u - k[1]) * (1 - ( u - k[1] ) / ( k[2] - k[1] ) / 2 )) * k[0] / nc);
1058 else
1059 return( 1 );
1060 }
1061
psi_hmpl(double x,const double k[])1062 double psi_hmpl(double x, const double k[])
1063 {
1064 /*
1065 * psi() for Hampel's redescending psi function
1066 * constants (a, b, r) == k[0:2] s.t. slope of psi is 1 in the center
1067 */
1068 // double sx = sign(x), u = fabs(x); :
1069 double sx, u;
1070 if (x < 0) { sx = -1; u = -x; } else { sx = +1; u = x; }
1071
1072 if (u <= k[0])
1073 return( x );
1074 else if (u <= k[1])
1075 return sx * k[0];
1076 else if (u <= k[2])
1077 return sx * k[0] * (k[2] - u) / (k[2] - k[1]);
1078 else
1079 return 0.;
1080 }
1081
psip_hmpl(double x,const double k[])1082 double psip_hmpl(double x, const double k[])
1083 {
1084 /*
1085 * psi'() for Hampel's redescending psi function
1086 * constants (a, b, r) == k[0:2] s.t. slope of psi is 1 in the center
1087 */
1088 double u = fabs(x);
1089
1090 if (u <= k[0])
1091 return( 1 );
1092 else if (u <= k[1])
1093 return( 0 );
1094 else if (u <= k[2])
1095 return( k[0] / ( k[1] - k[2] ) );
1096 else
1097 return( 0 );
1098 }
1099
psi2_hmpl(double x,const double k[])1100 double psi2_hmpl(double x, const double k[])
1101 {
1102 /*
1103 * psi''() for Hampel's redescending psi function
1104 * constants (a, b, r) == k[0:2] s.t. slope of psi is 1 in the center
1105 */
1106 return 0.; // even though psi'() is already discontinuous at k[j]
1107 }
1108
wgt_hmpl(double x,const double k[])1109 double wgt_hmpl(double x, const double k[])
1110 {
1111 /*
1112 * w(x) = psi(x)/x for Hampel's redescending psi function
1113 * Hampel redescending psi function
1114 * constants (a, b, r) == k[0:2] s.t. slope of psi is 1 in the center
1115 */
1116 double u = fabs(x);
1117
1118 if (u <= k[0])
1119 return( 1 );
1120 else if (u <= k[1])
1121 return( k[0] / u );
1122 else if (u <= k[2])
1123 return( k[0] * ( k[2] - u ) / ( k[2] - k[1] ) / u );
1124 else
1125 return( 0 );
1126 }
1127
1128
1129 //--- GGW := Generalized Gauss-Weight Koller and Stahel (2011)
1130 //--- ---
1131 // rho() & chi() need to be calculated by numerical integration -- apart from 6 pre-stored cases
rho_ggw(double x,const double k[])1132 double rho_ggw(double x, const double k[])
1133 {
1134 /*
1135 * Gauss Weight with constant center
1136 */
1137 if (k[0] > 0) { // for hard-coded constants --- use a *polynomial* approximation
1138 const double C[6][20] = { // 0: b = 1, 95% efficiency
1139 {0.094164571656733, -0.168937372816728, 0.00427612218326869,
1140 0.336876420549802, -0.166472338873754, 0.0436904383670537,
1141 -0.00732077121233756, 0.000792550423837942, -5.08385693557726e-05,
1142 1.46908724988936e-06, -0.837547853001024, 0.876392734183528,
1143 -0.184600387321924, 0.0219985685280105, -0.00156403138825785,
1144 6.16243137719362e-05, -7.478979895101e-07, -3.99563057938975e-08,
1145 1.78125589532002e-09, -2.22317669250326e-11},
1146 // 1: b = 1, 85% efficiency
1147 {0.174505224068561, -0.168853188892986, 0.00579250806463694,
1148 0.624193375180937, -0.419882092234336, 0.150011303015251,
1149 -0.0342185249354937, 0.00504325944243195, -0.0004404209084091,
1150 1.73268448820236e-05, -0.842160072154898, 1.19912623576069,
1151 -0.345595777445623, 0.0566407000764478, -0.00560501531439071,
1152 0.000319084704541442, -7.4279004383686e-06, -2.02063746721802e-07,
1153 1.65716101809839e-08, -2.97536178313245e-10},
1154 // 2: b = 1, bp 0.5
1155 {1.41117142330711, -0.168853741371095, 0.0164713906344165,
1156 5.04767833986545, -9.65574752971554, 9.80999125035463,
1157 -6.36344090274658, 2.667031271863, -0.662324374141645,
1158 0.0740982983873332, -0.84794906554363, 3.4315790970352,
1159 -2.82958670601597, 1.33442885893807, -0.384812004961396,
1160 0.0661359078129487, -0.00557221619221031, -5.42574872792348e-05,
1161 4.92564168111658e-05, -2.80432020951381e-06},
1162 // 3: b = 1.5, 95% efficiency
1163 {0.104604570079252, 0.0626649856211545, -0.220058184826331,
1164 0.403388189975896, -0.213020713708997, 0.102623342948069,
1165 -0.0392618698058543, 0.00937878752829234, -0.00122303709506374,
1166 6.70669880352453e-05, 0.632651530179424, -1.14744323908043,
1167 0.981941598165897, -0.341211275272191, 0.0671272892644464,
1168 -0.00826237596187364, 0.0006529134641922, -3.23468516804340e-05,
1169 9.17904701930209e-07, -1.14119059405971e-08},
1170 // 4: b = 1.5, 85% efficiency
1171 {0.205026436642222, 0.0627464477520301, -0.308483319391091,
1172 0.791480474953874, -0.585521414631968, 0.394979618040607,
1173 -0.211512515412973, 0.0707208739858416, -0.0129092527174621,
1174 0.000990938134086886, 0.629919019245325, -1.60049136444912,
1175 1.91903069049618, -0.933285960363159, 0.256861783311473,
1176 -0.0442133943831343, 0.00488402902512139, -0.000338084604725483,
1177 1.33974565571893e-05, -2.32450916247553e-07},
1178 // 5: b = 1.5, bp 0.5
1179 {1.35010856132000, 0.0627465630782482, -0.791613168488525,
1180 5.21196700244212, -9.89433796586115, 17.1277266427962,
1181 -23.5364159883776, 20.1943966645350, -9.4593988142692,
1182 1.86332355622445, 0.62986381140768, -4.10676399816156,
1183 12.6361433997327, -15.7697199271455, 11.1373468568838,
1184 -4.91933095295458, 1.39443093325178, -0.247689078940725,
1185 0.0251861553415515, -0.00112130382664914}};
1186
1187 double end[6] = {18.5527638190955, 13.7587939698492, 4.89447236180905,
1188 11.4974874371859, 8.15075376884422, 3.17587939698492};
1189 int j = ((int)k[0]) - 1;
1190 double c;
1191 switch(j) {
1192 // c : identical numbers to those in SET_ABC_GGW below
1193 case 0: c = 1.694; break;
1194 case 1: c = 1.2442567; break;
1195 case 2: c = 0.4375470; break;
1196 case 3: c = 1.063; break;
1197 case 4: c = 0.7593544; break;
1198 case 5: c = 0.2959132; break;
1199 default: error(_("rho_ggw(): case (%i) not implemented."), j+1);
1200 }
1201 x = fabs(x);
1202 if (x <= c)
1203 return(C[j][0]*x*x);
1204 else if (x <= 3*c)
1205 return(C[j][1] +
1206 x*(C[j][2] +
1207 x*(C[j][3] +
1208 x*(C[j][4] +
1209 x*(C[j][5] +
1210 x*(C[j][6] +
1211 x*(C[j][7] +
1212 x*(C[j][8] +
1213 x*(C[j][9])))))))));
1214 else if (x <= end[j])
1215 return(C[j][10] +
1216 x*(C[j][11] +
1217 x*(C[j][12] +
1218 x*(C[j][13] +
1219 x*(C[j][14] +
1220 x*(C[j][15] +
1221 x*(C[j][16] +
1222 x*(C[j][17] +
1223 x*(C[j][18]+
1224 x*(C[j][19]))))))))));
1225 else return(1.);
1226 }
1227 else { // k[0] == 0; k[1:4] = (a, b, c, rho(Inf)) = "general parameters"
1228 x = fabs(x);
1229 double a = 0., epsabs = R_pow(DOUBLE_EPS, 0.25), result, abserr;
1230 int neval, ier, last, limit = 100, lenw = 4 * limit;
1231 int *iwork = (int *) R_alloc(limit, sizeof(int));
1232 double *work = (double *) R_alloc(lenw, sizeof(double));
1233
1234 // --> calculate integral of psi(.); Rdqags() is from R's official API ("Writing R Extensions")
1235 Rdqags(psi_ggw_vec, (void *)k, &a, &x, &epsabs, &epsabs,
1236 &result, &abserr, &neval, &ier,
1237 &limit, &lenw, &last,
1238 iwork, work);
1239 if (ier >= 1)
1240 error(_("Error from Rdqags(psi_ggw*, k, ...): ier = %i"), ier);
1241 return(result/k[4]);
1242 }
1243 }
1244
psi_ggw_vec(double * x,int n,void * k)1245 void psi_ggw_vec(double *x, int n, void *k)
1246 {
1247 for (int i = 0; i<n; i++) x[i] = psi_ggw(x[i], k);
1248 }
1249
psi_ggw(double x,const double k[])1250 double psi_ggw(double x, const double k[])
1251 {
1252 /*
1253 * Gauss Weight with constant center
1254 */
1255 #define SET_ABC_GGW(NAME) \
1256 /* set a,b,c */ \
1257 double a, b, c; \
1258 switch((int)k[0]) { \
1259 /* user specified: */ \
1260 case 0: a = k[1]; b = k[2]; c = k[3]; break; \
1261 /* Set of predefined cases: */ \
1262 case 1: a = 0.648; b = 1.; c = 1.694; break; \
1263 case 2: a = 0.4760508; b = 1.; c = 1.2442567; break; \
1264 case 3: a = 0.1674046; b = 1.; c = 0.4375470; break; \
1265 case 4: a = 1.387; b = 1.5; c = 1.063; break; \
1266 case 5: a = 0.8372485; b = 1.5; c = 0.7593544; break; \
1267 case 6: a = 0.2036741; b = 1.5; c = 0.2959132; break; \
1268 default: error(#NAME "_ggw: Case not implemented."); \
1269 } \
1270 double ax = fabs(x);
1271
1272 SET_ABC_GGW(psi);
1273 if (ax < c) return x;
1274 else {
1275 a = -R_pow(ax-c,b)/2/a;
1276 return (a < MIN_Exp) ? 0. : x*exp(a);
1277 }
1278 }
1279
psip_ggw(double x,const double k[])1280 double psip_ggw(double x, const double k[])
1281 {
1282 /*
1283 * Gauss Weight with constant center
1284 */
1285 SET_ABC_GGW(psip);
1286 if (ax < c) return 1.;
1287 else {
1288 double ea;
1289 a *= 2.;
1290 ea = -R_pow(ax-c,b)/a;
1291 return (ea < MIN_Exp) ? 0. : exp(ea) * (1 - b/a*ax*R_pow(ax-c,b-1));
1292 }
1293 }
1294
wgt_ggw(double x,const double k[])1295 double wgt_ggw(double x, const double k[])
1296 {
1297 /*
1298 * Gauss Weight with constant center
1299 */
1300 SET_ABC_GGW(wgt);
1301 return (ax < c) ? 1. : exp(-R_pow(ax-c,b)/2/a);
1302 }
1303 #undef SET_ABC_GGW
1304
1305
1306 //--- LQQ := Linear-Quadratic-Quadratic ("lqq") --------------------------------
1307 //--- --- was LGW := "lin psip" := piecewise linear psi'() ------------------
1308
1309 // k[0:2] == (b, c, s) :
1310 // k[0]= b = bend adjustment
1311 // k[1]= c = cutoff of central linear part
1312 // k[2]= s : "slope of descending": 1 - s = min_x psi'(x) =: ms
1313
1314 // "lin psip" := piecewise linear psi'() :
psip_lqq(double x,const double k[])1315 double psip_lqq (double x, const double k[])
1316 {
1317 double ax = fabs(x);
1318 if (ax <= k[1])
1319 return(1.);
1320 else {
1321 double k01 = k[0] + k[1];// = b+c
1322 if (/*k[1] < ax && */ ax <= k01)
1323 return 1. - k[2]/k[0] * (ax - k[1]);
1324 else {
1325 double
1326 s5 = 1. - k[2], // = (1-s)
1327 a = (k[0] * k[2] - 2 * k01)/ s5;
1328 if (/* k01 < ax && */ ax < k01 + a)
1329 return -s5*((ax - k01)/a -1.);
1330 else
1331 return 0.;
1332 }
1333 }
1334 }
1335
1336 // piecewise linear psi'() ==> piecewise constant psi''():
psi2_lqq(double x,const double k[])1337 double psi2_lqq (double x, const double k[])
1338 {
1339 // double sx = sign(x), ax = fabs(x); :
1340 double sx, ax;
1341 if (x < 0) { sx = -1; ax = -x; } else { sx = +1; ax = x; }
1342
1343 // k[0:2] == (b, c, s) :
1344 if (ax <= k[1])
1345 return 0.;
1346 else {
1347 double k01 = k[0] + k[1];
1348 if (/*k[1] < ax && */ ax <= k01)
1349 return sx * (- k[2]/k[0]);
1350 else {
1351 double
1352 s5 = 1. - k[2], // = (1-s)
1353 a = (k[0] * k[2] - 2 * k01)/ s5;
1354 if (/* k01 < ax && */ ax < k01 + a)
1355 return sx * (- s5 / a);
1356 else
1357 return 0.;
1358 }
1359 }
1360 }
1361
psi_lqq(double x,const double k[])1362 double psi_lqq (double x, const double k[])
1363 {
1364 double ax = fabs(x);
1365 if (ax <= k[1])
1366 return(x);
1367 else {
1368 // k[0:2] == (b, c, s) :
1369 double k01 = k[0] + k[1];
1370 if (ax <= k01)
1371 return((double) (x>0 ? 1 : (x<0 ? -1 : 0)) *
1372 (ax - k[2] * pow(ax - k[1], 2.) / k[0] / 2.));
1373 else {
1374 double
1375 s5 = k[2] - 1., // s - 1
1376 s6 = -2 * k01 + k[0] * k[2]; // numerator( -a ) ==> s6/s5 = -a
1377 if (/* k01 < ax && */ ax < k01 - s6 / s5)
1378 return((double) (x>0 ? 1 : -1) *
1379 (-s6/2. - pow(s5, 2.) / s6 * (pow(ax - k01, 2.) / 2. + s6 / s5 * (ax - k01))));
1380 else
1381 return 0.;
1382 }
1383 }
1384 }
1385
rho_lqq(double x,const double k[])1386 double rho_lqq (double x, const double k[])
1387 {
1388 double ax = fabs(x), k01 = k[0] + k[1];
1389 if (ax <= k[1])
1390 return((3. * k[2] - 3.) /
1391 (k[2] * k[1] * (3. * k[1] + 2. * k[0]) +
1392 pow(k01, 2.)) * x * x);
1393 else if (/* k[1] < ax && */ ax <= k01) {
1394 double s0 = ax - k[1];
1395 return((6. * k[2] - 6.) /
1396 (k[2] * k[1] * (3. * k[1] + 2. * k[0]) + pow(k01, 2.)) *
1397 (x * x / 2. - k[2] / k[0] * pow(s0, 3.) / 6.));
1398 }
1399 else {
1400 double
1401 s5 = k[2] - 1.,
1402 s6 = -2 * k01 + k[0] * k[2];
1403 if (/* k01 < ax && */ ax < k01 - s6 / s5) {
1404 double s7 = ax - k01, k01_2 = pow(k01, 2.);
1405 return((6. * s5) /
1406 (k[2] * k[1] * (3. * k[1] + 2. * k[0]) + k01_2) *
1407 (k01_2 / 2. - k[2] * k[0] * k[0] / 6. -
1408 s7/2. * (s6 + s7 * (s5 + s7 * s5 * s5 / 3. / s6))));
1409 }
1410 else
1411 return 1.;
1412 }
1413 }
1414
wgt_lqq(double x,const double k[])1415 double wgt_lqq (double x, const double k[])
1416 {
1417 double ax = fabs(x);
1418 if (ax <= k[1])
1419 return(1.);
1420 else {
1421 double k01 = k[0] + k[1];
1422 if (ax <= k01) {
1423 double s0 = ax - k[1];
1424 return(1. - k[2] * s0 * s0 / (2 * ax * k[0]));
1425 }
1426 else {
1427 double
1428 s5 = k[2] - 1.,
1429 s6 = -2 * k01 + k[0] * k[2];
1430 if (ax < k01 - s6 / s5) {
1431 double s7 = ax - k01;
1432 return(-(s6/2. + s5 * s5 / s6 * s7 * (s7/2. + s6 / s5)) / ax);
1433 }
1434 else
1435 return(0.);
1436 }
1437 }
1438 }
1439 /*============================================================================*/
1440
1441
1442 /* this function finds the k-th place in the
1443 * vector a, in the process it permutes the
1444 * elements of a
1445 */
kthplace(double * a,int n,int k)1446 double kthplace(double *a, int n, int k)
1447 {
1448 int jnc,j;
1449 int l,lr;
1450 double ax,w;
1451 k--;
1452 l=0;
1453 lr=n-1;
1454 while (l < lr) {
1455 ax=a[k];
1456 jnc=l;
1457 j=lr;
1458 while (jnc <= j) {
1459 while (a[jnc] < ax) jnc++;
1460 while (a[j] > ax) j--;
1461 if (jnc <= j) {
1462 w=a[jnc];
1463 a[jnc]=a[j];
1464 a[j]=w;
1465 jnc++;
1466 j--;
1467 }
1468 }
1469 if (j < k) l=jnc;
1470 if (k < jnc) lr=j;
1471 }
1472 return(a[k]);
1473 }
1474
1475 /* This is from VR's bundle, MASS package VR/MASS/src/lqs.c : */
1476 /*
1477 Sampling k from 0:n-1 without replacement.
1478 */
sample_noreplace(int * x,int n,int k,int * ind_space)1479 static void sample_noreplace(int *x, int n, int k, int *ind_space)
1480 {
1481 int i, j, nn=n;
1482 #define II ind_space
1483
1484 for (i = 0; i < n; i++) II[i] = i;
1485 for (i = 0; i < k; i++) {
1486 j = nn * unif_rand();
1487 x[i] = II[j];
1488 II[j] = II[--nn];
1489 }
1490 #undef II
1491 }
1492
1493 /* RWLS iterations starting from i_estimate,
1494 * ---- the workhorse of the "lmrob_MM" algorithm, called only from R_lmrob_MM(),
1495 * which itself is called only from R's lmrob..M..fit().
1496 * In itself, ``just'' an M-estimator :
1497 */
rwls(const double X[],const double y[],int n,int p,double * estimate,double * i_estimate,double * resid,double * loss,double scale,double epsilon,int * max_it,double * rho_c,const int ipsi,int trace_lev)1498 Rboolean rwls(const double X[], const double y[], int n, int p,
1499 double *estimate, double *i_estimate,
1500 double *resid, double* loss,
1501 double scale, double epsilon,
1502 int *max_it, /* on Input: maximal number of iterations;
1503 on Output: number of iterations */
1504 double *rho_c, const int ipsi, int trace_lev)
1505 {
1506 int lwork = -1, one = 1, info = 1;
1507 double work0, *work, wtmp, *weights;
1508 double done = 1., dmone = -1., d_beta = 0.;
1509 int j, k, iterations = 0;
1510 Rboolean converged = FALSE;
1511
1512 double
1513 *wx = (double *) R_alloc(n*p, sizeof(double)),
1514 *wy = (double *) R_alloc(n, sizeof(double)),
1515 *beta0 = (double *) R_alloc(p, sizeof(double));
1516
1517 INIT_WLS(wx, wy, n, p);
1518
1519 COPY(i_estimate, beta0, p);
1520 /* calculate residuals */
1521 COPY(y, resid, n);
1522 F77_CALL(dgemv)("N", &n, &p, &dmone, X, &n, beta0, &one, &done, resid, &one FCONE);
1523
1524 /* main loop */
1525 while(!converged && ++iterations < *max_it) {
1526 R_CheckUserInterrupt();
1527 /* compute weights */
1528 get_weights_rhop(resid, scale, n, rho_c, ipsi, weights);
1529 /* solve weighted least squares problem */
1530 COPY(y, wy, n);
1531 FIT_WLS(X, wx, wy, n, p);
1532 COPY(wy, estimate, p);
1533 /* calculate residuals */
1534 COPY(y, resid, n);
1535 F77_CALL(dgemv)("N", &n, &p, &dmone, X, &n, estimate, &one, &done, resid, &one FCONE);
1536 d_beta = norm1_diff(beta0,estimate, p);
1537 if(trace_lev >= 3) {
1538 /* get the loss for the new estimate */
1539 *loss = sum_rho_sc(resid,scale,n,0,rho_c,ipsi);
1540 Rprintf(" it %4d: L(b1) = %#.12g ", iterations, *loss);
1541 if(trace_lev >= 4) {
1542 Rprintf("\n b1 = (");
1543 for(j=0; j < p; j++)
1544 Rprintf("%s%.11g", (j > 0)? ", " : "", estimate[j]);
1545 Rprintf(");");
1546 }
1547 Rprintf(" ||b0 - b1||_1 = %g\n", d_beta);
1548 }
1549 /* check for convergence */
1550 converged = d_beta <= epsilon * fmax2(epsilon, norm1(estimate, p));
1551 COPY(estimate, beta0, p);
1552 } /* end while(!converged & iter <=...) */
1553
1554 if(0 < trace_lev) {
1555 if(trace_lev < 3) *loss = sum_rho_sc(resid,scale,n,0,rho_c,ipsi);
1556 Rprintf(" rwls() used %2d it.; last ||b0 - b1||_1 = %#g, L(b1) = %.12g; %sconvergence\n",
1557 iterations, d_beta, *loss, (converged ? "" : "NON-"));
1558 }
1559
1560 *max_it = iterations;
1561
1562 CLEANUP_WLS;
1563
1564 return converged;
1565 } /* rwls() */
1566
1567 /* sets the entries of a matrix to zero */
zero_mat(double ** a,int n,int m)1568 void zero_mat(double **a, int n, int m)
1569 {
1570 int i,j;
1571 for(i=0; i < n; i++)
1572 for(j=0; j < m; j++)
1573 a[i][j] = 0.;
1574 }
1575
1576 /*
1577 *
1578 * 2004 / 5 -- Matias Salibian-Barrera & Victor Yohai
1579 * Department of Statistics, University of British Columbia
1580 * matias@stat.ubc.ca
1581 * Department of Mathematics, University of Buenos Aires
1582 * vyohai@uolsinectis.com.ar
1583 *
1584 *
1585 * Reference: A fast algorithm for S-regression estimates,
1586 * 2005, Salibian-Barrera and Yohai.
1587 */
1588
1589 /* This function implements the "large n" strategy
1590 */
fast_s_large_n(double * X,double * y,int * nn,int * pp,int * nRes,int * max_it_scale,double * res,int * ggroups,int * nn_group,int * K,int * max_k,double rel_tol,double inv_tol,double scale_tol,int * converged,int * best_r,double * bb,double * rrhoc,int * iipsi,double * bbeta,double * sscale,int trace_lev,int mts,Rboolean ss)1591 void fast_s_large_n(double *X, double *y,
1592 int *nn, int *pp, int *nRes, int *max_it_scale, double *res,
1593 int *ggroups, int *nn_group,
1594 int *K, int *max_k, double rel_tol, double inv_tol, double scale_tol, int *converged,
1595 int *best_r, double *bb, double *rrhoc, int *iipsi,
1596 double *bbeta, double *sscale,
1597 int trace_lev, int mts, Rboolean ss)
1598 {
1599 /* *X = the n x p design matrix (incl. intercept if appropriate),
1600 * in column order as sent by R)
1601 * *y = the ( n ) response vector
1602 * *nn =: n = the length of y
1603 * *pp =: p = the number of columns in X
1604 * *nRes = number of re-sampling candidates to be used in each partition
1605
1606 * *ggroups = number of groups in which to split the
1607 * random subsample
1608 * *nn_group = size of each of the (*ggroups) groups
1609 * to use in the random subsample
1610 * *K = number of refining steps for each candidate (typically 1 or 2)
1611 * *max_k = number of refining steps for each candidate (typically 1 or 2)
1612 [used to be hard coded to MAX_ITER_REFINE_S = 50 ]
1613 * *rel_tol= convergence tolerance for iterative refinement iterations
1614 [used to be hard coded to EPS = 1e-7 ]
1615 * *converged: will become 0(FALSE) iff at least one of the best_r iterations
1616 * did not converge (in max_k steps to rel_tol precision)
1617 * *best_r = no. of best candidates to be iterated further ("refined")
1618 * *bb = right-hand side of S-equation (typically 1/2)
1619 * *rrhoc = tuning constant for loss function
1620 * (this should be associated with *bb)
1621 * *iipsi = indicator for type of psi function to be used
1622 * *bbeta = final estimator
1623 * *sscale = associated scale estimator (or -1 when problem)
1624 */
1625 int i,j,k, ij, freedsamp = 0, initwls = 0;
1626 int n = *nn, p = *pp, kk = *K, ipsi = *iipsi;
1627 int groups = *ggroups, n_group = *nn_group, sg = groups * n_group;
1628 double b = *bb, sc, best_sc, worst_sc;
1629 /* (Pointers to) Arrays - to be allocated */
1630 int *indices, *ind_space;
1631 double **best_betas, *best_scales;
1632 double *xsamp, *ysamp, *beta_ref;
1633 double **final_best_betas, *final_best_scales;
1634
1635 #define CALLOC_MAT(_M_, _n_, _d_) \
1636 _M_ = (double **) Calloc(_n_, double *); \
1637 for(int i=0; i < _n_; i++) \
1638 _M_[i] = (double *) Calloc(_d_, double)
1639
1640 beta_ref = (double *) Calloc(p, double);
1641 CALLOC_MAT(final_best_betas, *best_r, p);
1642 final_best_scales = (double *) Calloc(*best_r, double);
1643 k = *best_r * groups;
1644 best_scales = (double *) Calloc(k, double );
1645 CALLOC_MAT(best_betas, k, p);
1646 indices = (int *) Calloc(sg, int);
1647 ind_space = (int *) Calloc(n, int);
1648 xsamp = (double *) Calloc(n_group*p, double);
1649 ysamp = (double *) Calloc(n_group, double);
1650
1651 /* assume that n > 2000 */
1652
1653 /* set the seed */
1654 GetRNGstate();
1655
1656 /* get a sample of k indices */
1657 sample_noreplace(indices, n, sg, ind_space);
1658 /* FIXME: define groups using nonsingular subsampling? */
1659 /* would also need to allow observations to be part */
1660 /* of multiple groups at the same time */
1661 Free(ind_space);
1662 /* FIXME: Also look at lqs_setup(),
1663 * ----- and xr[.,.] "fortran-like" matrix can be used from there!*/
1664
1665 /* For each (of 'groups') group : get the *best_r best betas : */
1666
1667 #define X(_k_, _j_) X[_j_*n+_k_]
1668 #define xsamp(_k_, _j_) xsamp[_j_*n_group+_k_]
1669
1670 for(i=0; i < groups; i++) {
1671 /* populate matrix */
1672 for(j = 0; j < n_group; j++) {
1673 ij = i*n_group + j;
1674 for (k = 0; k < p; k++)
1675 xsamp(j, k) = X(indices[ij], k);
1676 ysamp[j] = y[indices[ij]];
1677 }
1678 if (trace_lev)
1679 Rprintf(" Subsampling to find candidate betas in group %d:\n", i);
1680 if(fast_s_with_memory(xsamp, ysamp, res,
1681 &n_group, pp, nRes, max_it_scale, K, max_k,
1682 rel_tol, inv_tol, scale_tol,
1683 trace_lev, best_r, bb, rrhoc,
1684 iipsi, best_betas + i* *best_r,
1685 best_scales+ i* *best_r, mts, ss)) {
1686 *sscale = -1.; /* problem */
1687 goto cleanup_and_return;
1688 }
1689 }
1690
1691 Free(xsamp); Free(ysamp); freedsamp = 1;
1692 #undef xsamp
1693
1694 /* now iterate (refine) these "best_r * groups"
1695 * best betas in the (xsamp,ysamp) sample
1696 * with kk C-steps and keep only the "best_r" best ones
1697 */
1698 /* initialize new work matrices */
1699 double *wx, *wy;
1700 wx = (double *) R_alloc(n*p, sizeof(double)); // need only k here,
1701 wy = (double *) R_alloc(n, sizeof(double)); // but n in the last step
1702 xsamp = (double *) Calloc(sg*p, double);
1703 ysamp = (double *) Calloc(sg, double);
1704 freedsamp = 0;
1705
1706 #define xsamp(_k_,_j_) xsamp[_j_*sg+_k_]
1707
1708 for (ij = 0; ij < sg; ij++) {
1709 for (k = 0; k < p; k++)
1710 xsamp(ij, k) = X(indices[ij],k);
1711 ysamp[ij] = y[indices[ij]];
1712 }
1713
1714 int lwork = -1, one = 1, info = 1;
1715 double work0, *work, *weights;
1716 INIT_WLS(wx, wy, n, p); initwls = 1;
1717
1718 Rboolean conv = FALSE;
1719 int pos_worst_scale = 0;
1720 for(i=0; i < *best_r; i++)
1721 final_best_scales[i] = INFI;
1722 worst_sc = INFI;
1723 /* set the matrix to zero */
1724 zero_mat(final_best_betas, *best_r, p);
1725 for(i=0; i < (*best_r * groups); i++) {
1726 if(trace_lev >= 3) {
1727 Rprintf(" Sample[%3d]: before refine_(*, conv=FALSE):\n", i);
1728 if(i > 0) {
1729 Rprintf(" beta_cand : "); disp_vec(best_betas[i],p);
1730 Rprintf(" with scale %.15g\n", best_scales[i]);
1731 }
1732 }
1733 refine_fast_s(xsamp, wx, ysamp, wy, weights, sg, p, res,
1734 work, lwork, best_betas[i],
1735 kk, &conv/* = FALSE*/, *max_k, rel_tol, trace_lev,
1736 b, rrhoc, ipsi, best_scales[i], /* -> */ beta_ref, &sc);
1737 if(trace_lev >= 3) {
1738 Rprintf(" after refine: beta_ref : "); disp_vec(beta_ref,p);
1739 Rprintf(" with scale %.15g\n", sc);
1740 }
1741 if ( sum_rho_sc(res, worst_sc, sg, p, rrhoc, ipsi) < b ) {
1742 int scale_iter = *max_it_scale;
1743 /* scale will be better */
1744 sc = find_scale(res, b, rrhoc, ipsi, sc, sg, p, &scale_iter, scale_tol, trace_lev >= 3);
1745 int k2 = pos_worst_scale;
1746 final_best_scales[ k2 ] = sc;
1747 COPY(beta_ref, final_best_betas[k2], p);
1748 pos_worst_scale = find_max(final_best_scales, *best_r);
1749 worst_sc = final_best_scales[pos_worst_scale];
1750 }
1751 }
1752
1753 Free(xsamp); Free(ysamp); freedsamp = 1;
1754
1755 /* now iterate the best "best_r"
1756 * betas in the whole sample until convergence (max_k, rel_tol)
1757 */
1758 best_sc = INFI; *converged = 1; k = 0;
1759 if(trace_lev)
1760 Rprintf(" Now refine() to convergence for %d very best ones:\n",
1761 *best_r);
1762
1763 for(i=0; i < *best_r; i++) {
1764 conv = TRUE;
1765 int it_k = refine_fast_s(X, wx, y, wy, weights, n, p, res,
1766 work, lwork, final_best_betas[i], kk,
1767 &conv/* = TRUE */, *max_k, rel_tol, trace_lev,
1768 b, rrhoc, ipsi, final_best_scales[i],
1769 /* -> */ beta_ref, &sc);
1770 if(trace_lev)
1771 Rprintf(" Best[%d]: %sconvergence (%d iter.)",
1772 i, conv ? "" : "NON ", it_k);
1773 if(best_sc > sc) {
1774 if(trace_lev)
1775 Rprintf(": -> improved scale to %.15g", sc);
1776 best_sc = sc;
1777 COPY(beta_ref, bbeta, p);
1778 }
1779 if (trace_lev) Rprintf("\n");
1780 if (!conv && *converged) *converged = 0;
1781 if (k < it_k) k = it_k;
1782 }
1783 *sscale = best_sc;
1784 *max_k = k;
1785
1786 /* Done. Now clean-up. */
1787
1788 cleanup_and_return:
1789 PutRNGstate();
1790
1791 Free(best_scales);
1792 k = *best_r * groups;
1793 for(i=0; i < k; i++) Free( best_betas[i] );
1794 Free(best_betas); Free(indices);
1795 for(i=0; i < *best_r; i++)
1796 Free(final_best_betas[i]);
1797 Free(final_best_betas);
1798 Free(final_best_scales);
1799 Free(beta_ref);
1800
1801 if (freedsamp == 0) {
1802 Free(xsamp); Free(ysamp);
1803 }
1804
1805 if (initwls) {
1806 CLEANUP_WLS;
1807 }
1808
1809 #undef X
1810 #undef xsamp
1811
1812 } /* fast_s_large_n() */
1813
fast_s_with_memory(double * X,double * y,double * res,int * nn,int * pp,int * nRes,int * max_it_scale,int * K,int * max_k,double rel_tol,double inv_tol,double scale_tol,int trace_lev,int * best_r,double * bb,double * rrhoc,int * iipsi,double ** best_betas,double * best_scales,int mts,Rboolean ss)1814 Rboolean fast_s_with_memory(double *X, double *y, double *res,
1815 int *nn, int *pp, int *nRes, int *max_it_scale,
1816 int *K, int *max_k, double rel_tol, double inv_tol, double scale_tol,
1817 int trace_lev, int *best_r, double *bb, double *rrhoc,
1818 int *iipsi, double **best_betas, double *best_scales,
1819 int mts, Rboolean ss)
1820 {
1821 /*
1822 * Called from fast_s_large_n(), the adjustment for large "n",
1823 * same as fast_s, but it returns the best_r best betas,
1824 * and their associated scales.
1825 *
1826 * x : an n x p design matrix (including intercept if appropriate)
1827 * y : an n vector
1828 * res : an n vector of residuals
1829 * *nn = n, *pp = p
1830 * *nRes = number of re-sampling candidates to be taken
1831 * *K = number of refining steps for each candidate
1832 * *best_r = number of (refined) to be retained for full iteration
1833 * *bb = right-hand side of the S-equation (typically 1/2)
1834 * *rrhoc = tuning constant for loss function
1835 * (this should be associated with *bb)
1836 * *iipsi = indicator for type of loss function to be used
1837 * *best_betas = returning the best ... coefficient vectors
1838 * *best_scales = returning their associated residual scales
1839 */
1840
1841 int i,j,k;
1842 int n = *nn, p = *pp, nResample = *nRes;
1843 Rboolean conv = FALSE,
1844 sing = FALSE; // sing = TRUE|FALSE the final result
1845 int ipsi = *iipsi;
1846 double b = *bb, sc, worst_sc = INFI;
1847 double work0, *weights, *work;
1848 int lwork = -1, one = 1, info = 1;
1849
1850 SETUP_SUBSAMPLE(n, p, X, 1);
1851 INIT_WLS(X, y, n, p);
1852
1853 double
1854 *wx = (double *) Calloc(n*p, double),
1855 *wy = (double *) Calloc(n, double),
1856 *beta_cand = (double *) Calloc(p, double),
1857 *beta_ref = (double *) Calloc(p, double);
1858
1859 for(i=0; i < *best_r; i++)
1860 best_scales[i] = INFI;
1861 int pos_worst_scale = 0;
1862
1863 /* resampling approximation */
1864
1865 for(i=0; i < nResample; i++) {
1866 R_CheckUserInterrupt();
1867 /* find a candidate */
1868 sing = (Rboolean) // 0 |--> FALSE (= success); {1,2} |-> TRUE
1869 subsample(Xe, y, n, p, beta_cand, ind_space, idc, idr, lu, v, pivot,
1870 Dr, Dc, rowequ, colequ, 1, mts, ss, inv_tol, 1);
1871 if (sing) {
1872 for (k=0; k< *best_r; k++) best_scales[i] = -1.;
1873 goto cleanup_and_return;
1874 }
1875 /* FIXME: is_ok ?? */
1876
1877 /* improve the re-sampling candidate */
1878
1879 /* conv = FALSE : do *K refining steps */
1880 refine_fast_s(X, wx, y, wy, weights, n, p, res,
1881 work, lwork, beta_cand, *K, &conv/* = FALSE*/,
1882 *max_k, rel_tol, trace_lev, b, rrhoc, ipsi, -1.,
1883 /* -> */ beta_ref, &sc);
1884
1885 /* FIXME: if sc ~ 0 ---> return beta_cand and be done */
1886
1887 if ( sum_rho_sc(res, worst_sc, n, p, rrhoc, ipsi) < b ) {
1888 int scale_iter = *max_it_scale;
1889 /* scale will be better */
1890 sc = find_scale(res, b, rrhoc, ipsi, sc, n, p,
1891 &scale_iter, scale_tol, trace_lev >= 3);
1892 k = pos_worst_scale;
1893 best_scales[ k ] = sc;
1894 for(j=0; j < p; j++)
1895 best_betas[k][j] = beta_ref[j];
1896 pos_worst_scale = find_max(best_scales, *best_r);
1897 worst_sc = best_scales[pos_worst_scale];
1898 if (trace_lev >= 2) {
1899 Rprintf(" Sample[%3d]: found new candidate with scale %.7g in %d iter.\n",
1900 i, sc, scale_iter);
1901 Rprintf(" worst scale is now %.7g\n", worst_sc);
1902 }
1903 }
1904 } /* for(i ) */
1905
1906 cleanup_and_return:
1907
1908 CLEANUP_SUBSAMPLE;
1909 CLEANUP_WLS;
1910
1911 Free(wx); Free(wy);
1912 Free(beta_cand); Free(beta_ref);
1913
1914 return sing;
1915 } /* fast_s_with_memory() */
1916
fast_s(double * X,double * y,int * nn,int * pp,int * nRes,int * max_it_scale,double * res,int * K,int * max_k,double rel_tol,double inv_tol,double scale_tol,int * converged,int * best_r,double * bb,double * rrhoc,int * iipsi,double * bbeta,double * sscale,int trace_lev,int mts,Rboolean ss)1917 void fast_s(double *X, double *y,
1918 int *nn, int *pp, int *nRes, int *max_it_scale, double *res,
1919 int *K, int *max_k, double rel_tol, double inv_tol, double scale_tol, int *converged,
1920 int *best_r, double *bb, double *rrhoc, int *iipsi,
1921 double *bbeta, double *sscale, int trace_lev, int mts, Rboolean ss)
1922 {
1923 /* *X = the n x p design matrix (incl. intercept if appropriate),
1924 * in column order as sent by R)
1925 * *y = the ( n ) response vector
1926 * *nn =: n = the length of y
1927 * *pp =: p = the number of columns in X
1928 * *nRes = number of re-sampling candidates to be taken
1929 * *K = number of refining steps for each candidate
1930 * *best_r = number of (refined) to be retained for full iteration
1931 * *converged: will become FALSE iff at least one of the best_r iterations
1932 * did not converge (in max_k steps to rel_tol precision)
1933 * *bb = right-hand side of the S-equation (typically 1/2)
1934 * *rrhoc = tuning constant for loss function
1935 * (this should be associated with *bb)
1936 * *iipsi = indicator for type of loss function to be used
1937 * *bbeta = final estimator
1938 * *sscale = associated scale estimator (or -1 when problem)
1939 */
1940 int i,k;
1941 int n = *nn, p = *pp, nResample = *nRes, ipsi = *iipsi;
1942 double b = *bb;
1943 double sc, best_sc, aux;
1944 int lwork = -1, one = 1, info = 1;
1945 double work0, *work, *weights;
1946
1947 /* Rprintf("fast_s %d\n", ipsi); */
1948
1949 SETUP_SUBSAMPLE(n, p, X, 0);
1950
1951 // More arrays, allocated:
1952 double
1953 *wx = (double *) R_alloc(n*p, sizeof(double)),
1954 *wy = (double *) R_alloc(n, sizeof(double)),
1955 *beta_cand = (double *) Calloc(p, double),
1956 *beta_ref = (double *) Calloc(p, double),
1957 *best_scales = (double *) Calloc(*best_r, double),
1958 // matrix:
1959 **best_betas = (double **) Calloc(*best_r, double *);
1960 for(i=0; i < *best_r; i++) {
1961 best_betas[i] = (double*) Calloc(p, double);
1962 best_scales[i] = INFI;
1963 }
1964
1965 INIT_WLS(wx, wy, n, p);
1966
1967 /* disp_mat(x, n, p); */
1968
1969 int pos_worst_scale = 0;
1970 Rboolean conv = FALSE;
1971 double worst_sc = INFI;
1972
1973 /* srand((long)*seed_rand); */
1974 GetRNGstate();
1975
1976 /* resampling approximation */
1977
1978 if (trace_lev)
1979 Rprintf(" Subsampling to find candidate betas:\n", i);
1980
1981 for(i=0; i < nResample; i++) {
1982
1983 R_CheckUserInterrupt();
1984 /* find a candidate */
1985 Rboolean sing = (Rboolean) // 0 |--> FALSE (= success); {1,2} |-> TRUE
1986 subsample(Xe, y, n, p, beta_cand, ind_space, idc, idr, lu, v, pivot,
1987 Dr, Dc, rowequ, colequ, 1, mts, ss, inv_tol, 1);
1988 if (sing) {
1989 *sscale = -1.;
1990 goto cleanup_and_return;
1991 }
1992 if (trace_lev >= 5) {
1993 Rprintf(" Sample[%3d]: idc = ", i); disp_veci(idc, p);
1994 }
1995
1996 /* disp_vec(beta_cand,p); */
1997
1998 /* improve the re-sampling candidate */
1999
2000 /* conv = FALSE : do *k refining steps */
2001 refine_fast_s(X, wx, y, wy, weights, n, p, res,
2002 work, lwork, beta_cand, *K, &conv/* = FALSE*/,
2003 *max_k, rel_tol, trace_lev, b, rrhoc, ipsi, -1.,
2004 /* -> */ beta_ref, &sc);
2005 if(trace_lev >= 3) {
2006 double del = norm_diff(beta_cand, beta_ref, p);
2007 Rprintf(" Sample[%3d]: after refine_(*, conv=FALSE):\n", i);
2008 Rprintf(" beta_ref : "); disp_vec(beta_ref,p);
2009 Rprintf(" with ||beta_ref - beta_cand|| = %.12g, --> sc = %.15g\n",
2010 del, sc);
2011 }
2012 if(fabs(sc) == 0.) { /* exact zero set by refine_*() */
2013 if(trace_lev >= 1)
2014 Rprintf(" Too many exact zeroes -> leaving refinement!\n");
2015 *sscale = sc;
2016 COPY(beta_cand, bbeta, p);
2017 goto cleanup_and_return;
2018 }
2019 if ( sum_rho_sc(res, worst_sc, n, p, rrhoc, ipsi) < b ) {
2020 int scale_iter = *max_it_scale;
2021 /* scale will be better */
2022 sc = find_scale(res, b, rrhoc, ipsi, sc, n, p,
2023 &scale_iter, scale_tol, trace_lev >= 3);
2024 k = pos_worst_scale;
2025 best_scales[ k ] = sc;
2026 COPY(beta_ref, best_betas[k], p);
2027 pos_worst_scale = find_max(best_scales, *best_r);
2028 worst_sc = best_scales[pos_worst_scale];
2029 if (trace_lev >= 2) {
2030 Rprintf(" Sample[%3d]: found new candidate with scale %.7g in %d iter.\n",
2031 i, sc, scale_iter);
2032 Rprintf(" worst scale is now %.7g\n", worst_sc);
2033 }
2034 }
2035
2036 } /* for(i ) */
2037
2038 /* now look for the very best */
2039 if(trace_lev)
2040 Rprintf(" Now refine() to convergence for %d very best ones:\n",
2041 *best_r);
2042
2043 best_sc = INFI; *converged = 1; k = 0;
2044 for(i=0; i < *best_r; i++) {
2045 conv = TRUE;
2046 if(trace_lev >= 4) Rprintf(" i=%d:\n", i);
2047 int it_k = refine_fast_s(X, wx, y, wy, weights, n, p, res, work, lwork,
2048 best_betas[i], *K, &conv /* = TRUE */, *max_k,
2049 rel_tol, trace_lev, b, rrhoc, ipsi,
2050 best_scales[i], /* -> */ beta_ref, &aux);
2051 if(trace_lev)
2052 Rprintf(" Best[%d]: %sconvergence (%d iter.)",
2053 i, (conv) ? "" : "NON ", it_k);
2054 if(aux < best_sc) {
2055 if(trace_lev)
2056 Rprintf(": -> improved scale to %.15g", aux);
2057 best_sc = aux;
2058 COPY(beta_ref, bbeta, p);
2059 }
2060 if(trace_lev) Rprintf("\n");
2061 if (!conv && *converged) *converged = 0;
2062 if (k < it_k) k = it_k;
2063 }
2064 *sscale = best_sc;
2065 *max_k = k;
2066
2067 cleanup_and_return:
2068
2069 PutRNGstate();
2070
2071 CLEANUP_SUBSAMPLE;
2072 CLEANUP_WLS;
2073
2074 Free(best_scales);
2075 Free(beta_cand);
2076 Free(beta_ref);
2077 for(i=0; i < *best_r; i++)
2078 Free(best_betas[i]);
2079 Free(best_betas);
2080
2081 return;
2082 } /* fast_s() */
2083
refine_fast_s(const double X[],double * wx,const double y[],double * wy,double * weights,int n,int p,double * res,double * work,int lwork,double * beta_cand,int kk,Rboolean * conv,int max_k,double rel_tol,int trace_lev,double b,double * rrhoc,int ipsi,double initial_scale,double * beta_ref,double * scale)2084 int refine_fast_s(const double X[], double *wx, const double y[], double *wy,
2085 double *weights,
2086 int n, int p, double *res,
2087 double *work, int lwork, double *beta_cand,
2088 int kk, Rboolean *conv, int max_k, double rel_tol,
2089 int trace_lev,
2090 double b, double *rrhoc, int ipsi, double initial_scale,
2091 double *beta_ref, double *scale)
2092 {
2093 /*
2094 * X = matrix (n x p) of explanatory variables
2095 * y = vector ( n ) of responses
2096 * weights = robustness weights wt[] * y[] (of length n)
2097 * res = residuals y[] - x[,] * beta (of length n)
2098 * conv: FALSE means do kk refining steps (and conv stays FALSE)
2099 * TRUE means refine until convergence(rel_tol, max_k)
2100 * and in this case, 'conv' *returns* TRUE if refinements converged
2101 * beta_cand= candidate beta[] (of length p) Input *and* Output
2102 * is = initial scale input
2103
2104 * beta_ref = resulting beta[] (of length p) Output
2105 * scale = final scale Output
2106
2107 * for FIT_WLS, DGELS:
2108 * wx = matrix (n x p)
2109 * wy = vector of length n
2110 * work = vector of length lwork
2111 * lwork = length of vector work
2112 */
2113
2114 int i,j,k, zeroes=0, one = 1, info = 1;
2115 Rboolean converged = FALSE;/* Wall */
2116 double s0, done = 1., dmone = -1., wtmp;
2117
2118 if (trace_lev >= 4) {
2119 Rprintf(" beta_cand before refinement : "); disp_vec(beta_cand,p);
2120 }
2121
2122 /* calculate residuals */
2123 COPY(y, res, n);
2124 F77_CALL(dgemv)("N", &n, &p, &dmone, X, &n, beta_cand, &one, &done, res, &one FCONE);
2125 for(j=0; j < n; j++) {
2126 if( fabs(res[j]) < EPS_SCALE )
2127 zeroes++;
2128 }
2129 /* if "perfect fit", return it with a 0 assoc. scale */
2130 if( zeroes > (((double)n + (double)p)/2.) ) /* <<- FIXME: depends on 'b' ! */
2131 {
2132 COPY(beta_cand, beta_ref, p);
2133 *scale = 0.;
2134 return 0;
2135 }
2136
2137 if( initial_scale < 0. )
2138 initial_scale = MAD(res, n, 0., wy, weights);// wy and weights used as work arrays
2139 s0 = initial_scale;
2140 if(*conv)
2141 kk = max_k;
2142 for(i=0; i < kk; i++) {
2143 /* one step for the scale */
2144 s0 = s0 * sqrt( sum_rho_sc(res, s0, n, p, rrhoc, ipsi) / b );
2145 /* compute weights for IRWLS */
2146 get_weights_rhop(res, s0, n, rrhoc, ipsi, weights);
2147 /* solve weighted least squares problem */
2148 COPY(y, wy, n);
2149 FIT_WLS(X, wx, wy, n, p);
2150 COPY(wy, beta_ref, p);
2151 if(*conv) { /* check for convergence */
2152 double del = norm_diff(beta_cand, beta_ref, p);
2153 double nrmB= norm(beta_cand, p);
2154 if(trace_lev >= 4)
2155 Rprintf(" it %4d, ||b[i]||= %#.12g, ||b[i] - b[i-1]|| = %#.15g\n",
2156 i, nrmB, del);
2157 converged = (del <= rel_tol * fmax2(rel_tol, nrmB));
2158 if(converged)
2159 break;
2160 }
2161 /* calculate residuals */
2162 COPY(y, res, n);
2163 F77_CALL(dgemv)("N", &n, &p, &dmone, X, &n, beta_ref, &one, &done, res, &one FCONE);
2164 COPY(beta_ref, beta_cand, p);
2165 } /* for(i = 0; i < kk ) */
2166
2167 if(*conv) {
2168 /* was "if(0)", since default lead to 'NOT converged' */
2169 if(!converged) {
2170 *conv = FALSE;
2171 warning(_("S refinements did not converge (to refine.tol=%g) in %d (= k.max) steps"),
2172 rel_tol, i);
2173 }
2174 }
2175 *scale = s0;
2176 return i; /* number of refinement steps */
2177 } /* refine_fast_s() */
2178
2179 /* Subsampling part for M-S algorithm */
2180 /* Recreates RLFRSTML function found in src/lmrobml.f */
2181 /* of the robust package */
m_s_subsample(double * X1,double * y,int n,int p1,int p2,int nResample,int max_it_scale,double rel_tol,double inv_tol,double scale_tol,double * bb,double * rrhoc,int ipsi,double * sscale,int trace_lev,double * b1,double * b2,double * t1,double * t2,double * y_tilde,double * res,double * x1,double * x2,int * NIT,int * K,int * KODE,double * SIGMA,double * BET0,double * SC1,double * SC2,double * SC3,double * SC4,int mts,Rboolean ss)2182 void m_s_subsample(double *X1, double *y, int n, int p1, int p2,
2183 int nResample, int max_it_scale,
2184 double rel_tol, double inv_tol, double scale_tol, double *bb,
2185 double *rrhoc, int ipsi, double *sscale, int trace_lev,
2186 double *b1, double *b2, double *t1, double *t2,
2187 double *y_tilde, double *res, double *x1, double *x2,
2188 int *NIT, int *K, int *KODE, double *SIGMA, double *BET0,
2189 double *SC1, double *SC2, double *SC3, double *SC4, int mts, Rboolean ss)
2190 {
2191 int i, one = 1, p = p1 + p2, info;
2192 double b = *bb, sc = INFI, done = 1., dmone = -1.;
2193 *sscale = INFI;
2194
2195 if (trace_lev >= 2)
2196 Rprintf(" Starting subsampling procedure.. ");
2197
2198 SETUP_SUBSAMPLE(n, p2, x2, 0);
2199
2200 /* set the seed */
2201 GetRNGstate();
2202
2203 if (trace_lev >= 2) Rprintf(" [setup Ok]\n");
2204
2205 for(i=0; i < nResample; i++) {
2206 R_CheckUserInterrupt();
2207 /* STEP 1: Draw a subsample of size p2 from (X2, y) */
2208 Rboolean sing = (Rboolean) // 0 |--> FALSE (= success); {1,2} |-> TRUE
2209 subsample(Xe, y, n, p2, t2, ind_space, idc, idr, lu, v, pivot,
2210 Dr, Dc, rowequ, colequ, /* sample= */ TRUE, mts,
2211 ss, inv_tol, /*solve = */ TRUE);
2212 if (sing) {
2213 *sscale = -1.;
2214 goto cleanup_and_return;
2215 }
2216 /* calculate partial residuals */
2217 COPY(y, y_tilde, n);
2218 F77_CALL(dgemv)("N", &n, &p2, &dmone, x2, &n, t2, &one, &done, y_tilde, &one FCONE);
2219 /* STEP 3: Obtain L1-estimate of b1 */
2220 COPY(X1, x1, n*p1);
2221 F77_CALL(rllarsbi)(x1, y_tilde, &n, &p1, &n, &n, &rel_tol,
2222 NIT, K, KODE, SIGMA, t1, res, SC1, SC2,
2223 SC3, SC4, BET0);
2224 if (*KODE > 1) {
2225 REprintf("m_s_subsample(): Problem in RLLARSBI (RILARS). KODE=%d. Exiting.\n",
2226 *KODE);
2227 *sscale = -1.;
2228 goto cleanup_and_return;
2229 }
2230 /* STEP 4: Check if candidate looks promising */
2231 if (sum_rho_sc(res, *sscale, n, p, rrhoc, ipsi) < b) {
2232 int scale_iter = max_it_scale;
2233 /* scale will be better */
2234 /* STEP 5: Solve for sc */
2235 sc = find_scale(res, b, rrhoc, ipsi, sc, n, p, &scale_iter,
2236 scale_tol, trace_lev >= 4);
2237 if(trace_lev >= 2)
2238 Rprintf(" Sample[%3d]: new candidate with sc = %#10.5g in %d iter\n",
2239 i, sc, scale_iter);
2240 /* STEP 6: Update best fit */
2241 *sscale = sc;
2242 COPY(t1, b1, p1);
2243 COPY(t2, b2, p2);
2244 if (sc < EPS_SCALE) {
2245 REprintf("\nScale too small\n",
2246 "Aborting m_s_subsample()\n\n");
2247 *sscale = -1.;
2248 goto cleanup_and_return;
2249 }
2250 }
2251 } /* for(i ) */
2252
2253 /* STEP 7: Clean up and return */
2254 if (trace_lev >= 1) {
2255 Rprintf(" Finished M-S subsampling with scale = %.5f\n",*sscale);
2256 #define maybe_SHOW_b1_b2 \
2257 if (trace_lev >= 3) { \
2258 Rprintf(" b1: "); disp_vec(b1,p1);\
2259 Rprintf(" b2: "); disp_vec(b2,p2);\
2260 }
2261 maybe_SHOW_b1_b2;
2262 }
2263
2264 cleanup_and_return:
2265 CLEANUP_SUBSAMPLE;
2266 PutRNGstate();
2267 } /* m_s_subsample() */
2268
2269 /* Descent step for M-S algorithm
2270 * Return value: convergence; note that convergence is *not* guaranteed
2271 */
m_s_descent(double * X1,double * X2,double * y,int n,int p1,int p2,int K_m_s,int max_k,int max_it_scale,double rel_tol,double scale_tol,double * bb,double * rrhoc,int ipsi,double * sscale,int trace_lev,double * b1,double * b2,double * t1,double * t2,double * y_tilde,double * res,double * res2,double * x1,double * x2,int * NIT,int * K,int * KODE,double * SIGMA,double * BET0,double * SC1,double * SC2,double * SC3,double * SC4)2272 Rboolean m_s_descent(double *X1, double *X2, double *y,
2273 int n, int p1, int p2, int K_m_s, int max_k, int max_it_scale,
2274 double rel_tol, double scale_tol, double *bb, double *rrhoc, int ipsi,
2275 double *sscale, int trace_lev,
2276 double *b1, double *b2, double *t1, double *t2,
2277 double *y_tilde, double *res, double *res2, double *x1, double *x2,
2278 int *NIT, int *K, int *KODE, double *SIGMA, double *BET0,
2279 double *SC1, double *SC2, double *SC3, double *SC4)
2280 {
2281 int j, k, nnoimpr = 0, nref = 0;
2282 int p = p1 + p2;
2283 Rboolean converged = FALSE;
2284 double b = *bb;
2285 double sc = *sscale, done = 1., dmone = -1.;
2286 int lwork = -1, one = 1, info = 1;
2287 double work0, *work, wtmp, *weights;
2288
2289 COPY(b1, t1, p1);
2290 COPY(b2, t2, p2);
2291 COPY(res, res2, n);
2292
2293 if (trace_lev >= 2)
2294 Rprintf(" Starting descent procedure...\n");
2295
2296 INIT_WLS(x2, y, n, p2);
2297
2298 if (trace_lev >= 3) {
2299 Rprintf(" Scale: %.5f\n", *sscale);
2300 if (trace_lev >= 5) {
2301 Rprintf(" res2: "); disp_vec(res2,n);
2302 }
2303 }
2304
2305 /* Do descent steps until there is no improvement for */
2306 /* K_m_s steps or we are converged */
2307 /* (convergence is not guaranteed) */
2308 while ( (nref++ < max_k) & (!converged) & (nnoimpr < K_m_s) ) {
2309 R_CheckUserInterrupt();
2310 /* STEP 1: update b2 (save it to t2) */
2311 /* y_tilde = y - x1 %*% t1 */
2312 COPY(y, y_tilde, n);
2313 COPY(X1, x1, n*p1);
2314 F77_CALL(dgemv)("N", &n, &p1, &dmone, x1, &n, t1, &one, &done, y_tilde, &one FCONE);
2315 /* compute weights */
2316 get_weights_rhop(res2, sc, n, rrhoc, ipsi, weights);
2317 /* solve weighted least squares problem */
2318 FIT_WLS(X2, x2, y_tilde, n, p2);
2319 COPY(y_tilde, t2, p2);
2320 /* get (intermediate) residuals */
2321 COPY(y, res2, n);
2322 F77_CALL(dgemv)("N", &n, &p2, &dmone, X2, &n, t2, &one, &done, res2, &one FCONE);
2323 /* STEP 2: Obtain L1-estimate of b1 */
2324 COPY(res2, y_tilde, n);
2325 F77_CALL(rllarsbi)(x1, y_tilde, &n, &p1, &n, &n, &rel_tol,
2326 NIT, K, KODE, SIGMA, t1, res2,
2327 SC1, SC2, SC3, SC4, BET0);
2328 if (*KODE > 1) {
2329 CLEANUP_WLS;
2330 error(_("m_s_descent(): Problem in RLLARSBI (RILARS). KODE=%d. Exiting."),
2331 *KODE);
2332 }
2333 /* STEP 3: Compute the scale estimate */
2334 int scale_iter = max_it_scale;
2335 sc = find_scale(res2, b, rrhoc, ipsi, sc, n, p, &scale_iter,
2336 scale_tol, trace_lev >= 4); // <- here only if higher trace_lev
2337 /* STEP 4: Check for convergence */
2338 /* FIXME: check convergence using scale ? */
2339 double del = sqrt(norm_diff2(b1, t1, p1) + norm_diff2(b2, t2, p2));
2340 double nrmB = sqrt(norm2(t1, p1) + norm2(t2, p2));
2341 converged = (del < rel_tol * fmax2(rel_tol, nrmB));
2342 if (trace_lev >= 3) {
2343 if(converged) Rprintf(" -->> converged\n");
2344 if (trace_lev >= 4) {
2345 Rprintf(" Ref.step %3d: #{no-improvements}=%3d; (del,dB)=(%12.7g,%12.7g)\n",
2346 nref, nnoimpr, del, rel_tol * fmax2(rel_tol, nrmB));
2347 if (trace_lev >= 5) {
2348 Rprintf(" weights: "); disp_vec(weights,n);
2349 Rprintf(" t2: "); disp_vec(t2,p2);
2350 Rprintf(" t1: "); disp_vec(t1,p1);
2351 Rprintf(" res2: "); disp_vec(res2,n);
2352 }
2353 }
2354 }
2355 /* STEP 5: Update best fit */
2356 if (sc < *sscale) {
2357 COPY(t1, b1, p1);
2358 COPY(t2, b2, p2);
2359 COPY(res2, res, n);
2360 *sscale = sc;
2361 if (trace_lev >= 2)
2362 Rprintf(" Refinement step %3d: better fit, scale: %#10.5g\n",
2363 nref, sc);
2364 nnoimpr = 0;
2365 } else {
2366 if (trace_lev >= 3)
2367 Rprintf(" Refinement step %3d: no improvement, scale: %#10.5g\n",
2368 nref, sc);
2369 nnoimpr++;
2370 }
2371 } // while(.)
2372
2373 if ( (!converged) & (nref == max_k) )
2374 warning(_(" M-S estimate: maximum number of refinement steps reached."));
2375
2376 if (trace_lev >= 1) {
2377 Rprintf(" Descent procedure: %sconverged (best scale: %.5g, last step: %.5g)\n",
2378 converged ? "" : "not ", *sscale, sc);
2379 if (nnoimpr == K_m_s)
2380 Rprintf(" The procedure stopped after %d steps because there was no improvement in the last %d steps.\n To increase this number, use the control parameter 'k.m_s'.\n", nref, nnoimpr);
2381 else if (trace_lev >= 2)
2382 Rprintf(" No improvements in %d out of %d steps.\n", nnoimpr, nref);
2383 maybe_SHOW_b1_b2;
2384 }
2385
2386 CLEANUP_WLS;
2387
2388 return converged;
2389 } /* m_s_descent() */
2390
2391 /* draw a subsample of observations and calculate a candidate *
2392 * starting value for S estimates *
2393 * uses a custom LU decomposition, which acts on the transposed design *
2394 * matrix. In case of a singular subsample, the subsample is modified *
2395 * until it is non-singular (for ss == TRUE (== 1)). *
2396 * *
2397 * Parts of the algorithm are based on the Gaxpy version of the LU *
2398 * decomposition with partial pivoting by *
2399 * Golub G. H., Van Loan C. F. (1996) - MATRIX Computations */
subsample(const double x[],const double y[],int n,int m,double * beta,int * ind_space,int * idc,int * idr,double * lu,double * v,int * pivot,double * Dr,double * Dc,int rowequ,int colequ,Rboolean sample,int mts,Rboolean ss,double tol_inv,Rboolean solve)2400 int subsample(const double x[], const double y[], int n, int m,
2401 double *beta, int *ind_space, int *idc, int *idr,
2402 double *lu, double *v, int *pivot,
2403 double *Dr, double *Dc, int rowequ, int colequ,
2404 Rboolean sample, int mts, Rboolean ss, double tol_inv, Rboolean solve)
2405 {
2406 /* x: design matrix (n x m)
2407 y: response vector
2408 n: length of y, nrow of x
2409 m: ncol of x ( == p )
2410 beta: [out] candidate parameters (length m)
2411 ind_space: (required in sample_noreplace, length n)
2412 holds the index permutation
2413 idc: (required in sample_noreplace, !! length n !!)
2414 [out] index of observations used in subsample
2415 idr: work array of length m
2416 lu: [out] LU decomposition of subsample of xt (m x m)
2417 Note: U has is not rescaled by 1 / *cf, as it should,
2418 this is done R_subsample().
2419 v: work array of length m
2420 pivot: [out] pivoting table of LU decomposition (length m-1)
2421 Dr: row equilibration (as calculated in SETUP_EQUILIBRATION)
2422 Dc: column equilibration
2423 rowequ: whether rows were equilibrated
2424 coleq: whether cols were equilibrated
2425 sample: whether to sample or not
2426 mts: the number of singular samples allowed before
2427 giving up (Max Try Samples)
2428 ss: type of subsampling to be used:
2429 0 (FALSE): simple subsampling
2430 1 (TRUE): nonsingular subsampling
2431 tol_inv: tolerance for declaring a matrix singular
2432 solve: solve the least squares problem on the subsample?
2433 (0: no, 1: yes)
2434
2435 return value ('condition'):
2436 0: success
2437 1: singular (matrix xt does not contain a m dim. full rank submatrix)
2438 2: too many singular resamples (simple subsampling case)
2439 */
2440 int j, k, l, one = 1, mu = 0, tmpi, i = 0, attempt = 0;
2441 double tmpd;
2442 Rboolean sing;
2443
2444 #define xt(_k_, _j_) x[idr[_k_]*n+idc[_j_]]
2445 #define U(_k_, _j_) lu[_j_*m+_k_]
2446 #define u(_k_, _j_) lu + (_j_*m+_k_)
2447 #define L(_k_, _j_) lu[_j_*m+_k_]
2448 #define l(_k_, _j_) lu + (_j_*m+_k_)
2449
2450 Start:
2451 /* STEP 1: Calculate permutation of 1:n */
2452 if (sample) {
2453 sample_noreplace(ind_space, n, n, idc);
2454 } else for(k=0;k<n;k++) ind_space[k] = k;
2455 for(k=0;k<m;k++) idr[k] = k;
2456
2457 /* STEP 2: Calculate LU decomposition of the first m cols of xt *
2458 * using the order in ind_space */
2459 for(j = 0; j < m; j++) {
2460 sing=TRUE;
2461 do {
2462 if (i+j == n) {
2463 warning(_("subsample(): could not find non-singular subsample."));
2464 return(1);
2465 }
2466 idc[j] = ind_space[i+j];
2467 if (j == 0) {
2468 for(k=j;k<m;k++) v[k] = xt(k, j);
2469 } else {
2470 for(k=0;k<j;k++) U(k,j) = xt(k, j);
2471 /* z = solve(lu[0:(j-1), 0:(j-1)], xt[0:(j-1), j]) */
2472 F77_CALL(dtrsv)("L", "N", "U", &j, lu, &m, u(0, j), &one FCONE FCONE FCONE);
2473 /* Rprintf("Step %d: z = ", j); */
2474 /* for(i=0; i < j; i++) Rprintf("%lf ",U(i, j)); */
2475 /* Rprintf("\n"); */
2476 /* v[j:(m-1)] = xt[j:(m-1), j] - L[j:(m-1), 0:(j-1)] %*% z */
2477 for(k=j;k<m;k++) {
2478 v[k] = xt(k, j);
2479 for(l=0;l<j;l++) v[k] -= L(k, l) * U(l, j);
2480 }
2481 /* Rprintf("v = "); disp_vec(v, m); */
2482 }
2483 if (j < m-1) {
2484 /* find pivot */
2485 tmpd=fabs(v[j]); mu = j;
2486 for(k=j+1;k<m;k++) if (tmpd < fabs(v[k])) { mu = k; tmpd = fabs(v[k]); }
2487 /* debug possumDiv example, see tests/subsample.R */
2488 /* if (j == 36) { */
2489 /* Rprintf("Step %d: ", j+1); */
2490 /* for(k=j;k<m;k++) Rprintf("%lf ", fabs(v[k])); */
2491 /* Rprintf("\n %d %lf\n", mu+1, v[mu]); */
2492 /* Rprintf("47 > 51: %x\n", fabs(v[46]) > fabs(v[50])); */
2493 /* Rprintf("47 < 51: %x\n", fabs(v[46]) < fabs(v[50])); */
2494 /* } */
2495 /* continue only if pivot is large enough */
2496 if (tmpd >= tol_inv) {
2497 pivot[j] = mu;
2498 tmpd = v[j]; v[j] = v[mu]; v[mu] = tmpd;
2499 tmpi = idr[j]; idr[j] = idr[mu]; idr[mu] = tmpi;
2500 for(k=j+1;k<m;k++) L(k, j) = v[k] / v[j];
2501 if (j > 0) {
2502 for(k=0;k<j;k++) {
2503 tmpd = L(j, k); L(j, k) = L(mu, k); L(mu, k) = tmpd;
2504 }
2505 }
2506 }
2507 }
2508 if (fabs(v[j]) < tol_inv) {
2509 if (ss == 0) {
2510 attempt++;
2511 if (attempt >= mts) {
2512 warning(_("Too many singular resamples. Aborting subsample().\n See parameter 'subsampling; in help of lmrob.config()."));
2513 return(2);
2514 }
2515 goto Start;
2516 }
2517 /* drop observation and try next one */
2518 i++;
2519 } else {
2520 sing = FALSE;
2521 U(j, j) = v[j];
2522 }
2523 } while(sing);
2524 } /* end for loop */
2525
2526 /* Rprintf("lu:"); disp_vec(lu, m*m); */
2527 /* Rprintf("pivot:"); disp_veci(pivot, m-1); */
2528 /* Rprintf("idc:"); disp_veci(idc, m); */
2529
2530 /* STEP 3: Solve for candidate parameters if requested */
2531 if (solve == 0) {
2532 for(k=0;k<m;k++) beta[k] = NA_REAL;
2533 } else {
2534 for(k=0;k<m;k++) beta[k] = y[idc[k]];
2535 /* scale y ( = beta ) */
2536 if (rowequ) for(k=0;k<m;k++) beta[k] *= Dr[idc[k]];
2537 /* solve U\tr L\tr \beta = y[subsample] */
2538 F77_CALL(dtrsv)("U", "T", "N", &m, lu, &m, beta, &one FCONE FCONE FCONE);
2539 F77_CALL(dtrsv)("L", "T", "U", &m, lu, &m, beta, &one FCONE FCONE FCONE);
2540 /* scale the solution */
2541 if (colequ) for(k=0;k<m;k++) beta[k] *= Dc[idr[k]];
2542 /* undo pivoting */
2543 for(k=m-2;k>=0;k--) {
2544 tmpd = beta[k]; beta[k] = beta[pivot[k]]; beta[pivot[k]] = tmpd;
2545 }
2546 }
2547
2548 return(0);
2549
2550 #undef Xt
2551 #undef U
2552 #undef u
2553 #undef L
2554 #undef l
2555 }
2556
get_weights_rhop(const double r[],double s,int n,const double rrhoc[],int ipsi,double * w)2557 void get_weights_rhop(const double r[], double s, int n, const double rrhoc[], int ipsi,
2558 double *w)
2559 {
2560 for(int i=0; i < n; i++)
2561 w[i] = wgt(r[i] / s, rrhoc, ipsi);
2562 }
2563
find_scale(const double r[],double b,const double rrhoc[],int ipsi,double initial_scale,int n,int p,int * iter,double scale_tol,Rboolean trace)2564 double find_scale(const double r[], double b, const double rrhoc[], int ipsi,
2565 double initial_scale, int n, int p,
2566 int* iter, // input: max_iter, output: #{iterations used}
2567 double scale_tol,
2568 Rboolean trace)
2569 {
2570 if(initial_scale <= 0.) {
2571 warning(_("find_scale(*, initial_scale = %g) -> final scale = 0"), initial_scale);
2572 return 0.;
2573 }
2574 // else
2575 double scale = initial_scale;
2576 if(trace) Rprintf("find_scale(*, ini.scale =%#15.11g):\nit | new scale\n", scale);
2577 for(int it = 0; it < iter[0]; it++) {
2578 scale = initial_scale *
2579 sqrt( sum_rho_sc(r, initial_scale, n, p, rrhoc, ipsi) / b ) ;
2580 if(trace) Rprintf("%2d | %#13.10g\n", it, scale);
2581 if(fabs(scale - initial_scale) <= scale_tol*initial_scale) { // converged:
2582 *iter = it; return(scale);
2583 }
2584 initial_scale = scale;
2585 }
2586 warning(_("find_scale() did not converge in '%s' (= %d) iterations with tol=%g, last rel.diff=%g"),
2587 "maxit.scale", /* <- name from lmrob.control() */ *iter, scale_tol,
2588 (scale - initial_scale) / initial_scale);
2589
2590 return(scale);
2591 }
2592
2593
2594 // As R's which.max(a), return()ing zero-based k in {0,1,...,n-1}
find_max(double * a,int n)2595 int find_max(double *a, int n)
2596 {
2597 int k = 0;
2598 if(n > 1) {
2599 double tt = a[0];
2600 for(int i=1; i < n; i++)
2601 if(tt < a[i]) {
2602 tt = a[i];
2603 k = i;
2604 }
2605 }
2606 return k;
2607 }
2608
sum_rho_sc(const double r[],double scale,int n,int p,const double c[],int ipsi)2609 double sum_rho_sc(const double r[], double scale, int n, int p, const double c[], int ipsi)
2610 {
2611 double s = 0;
2612 for(int i=0; i < n; i++)
2613 s += rho(r[i]/scale, c, ipsi);
2614 return(s / ((double) n - p));
2615 }
2616
2617 /* ||x||_2^2 */
norm2(double * x,int n)2618 double norm2(double *x, int n)
2619 {
2620 double s = 0.;
2621 int one = 1;
2622 s = F77_CALL(dnrm2)(&n, x, &one);
2623 return( s*s );
2624 }
2625
2626 /* ||x||_2 */
norm(double * x,int n)2627 double norm(double *x, int n)
2628 {
2629 int one = 1;
2630 return(F77_CALL(dnrm2)(&n, x, &one));
2631 }
2632
norm1(double * x,int n)2633 double norm1(double *x, int n)
2634 {
2635 int one = 1;
2636 return(F77_CALL(dasum)(&n, x, &one));
2637 }
2638
2639 /* ||x-y||_2^2 */
norm_diff2(double * x,double * y,int n)2640 double norm_diff2(double *x, double *y, int n)
2641 {
2642 int i;
2643 double s = 0;
2644 for(i=0; i < n; i++)
2645 s += (x[i]-y[i])*(x[i]-y[i]);
2646 return( s );
2647 }
2648
2649 /* ||x-y||_2 */
norm_diff(double * x,double * y,int n)2650 double norm_diff(double *x, double *y, int n)
2651 {
2652 int i;
2653 double s = 0;
2654 for(i=0; i < n; i++)
2655 s += (x[i]-y[i])*(x[i]-y[i]);
2656 return( sqrt(s) );
2657 }
2658
2659
2660 /* ||x-y||_1 */
norm1_diff(double * x,double * y,int n)2661 double norm1_diff(double *x, double *y, int n)
2662 {
2663 int i;
2664 double s = 0;
2665 for(i=0; i < n; i++)
2666 s += fabs(x[i]-y[i]);
2667 return(s);
2668 }
2669
2670
MAD(double * a,int n,double center,double * b,double * tmp)2671 double MAD(double *a, int n, double center, double *b,
2672 double *tmp)
2673 {
2674 /* if center == 0 then do not center */
2675 int i;
2676 /* if( fabs(center) > 0.) { */
2677 for(i=0; i < n; i++)
2678 b[i] = a[i] - center;
2679 /* } */
2680 return( median_abs(b,n,tmp) * 1.4826 );
2681 }
2682
median(double * x,int n,double * aux)2683 double median(double *x, int n, double *aux)
2684 {
2685 double t;
2686 for(int i=0; i < n; i++) aux[i]=x[i];
2687 if ( (n/2) == (double) n / 2 )
2688 t = ( kthplace(aux,n,n/2) + kthplace(aux,n,n/2+1) ) / 2.0 ;
2689 else t = kthplace(aux,n, n/2+1 ) ;
2690 return(t);
2691 }
2692
median_abs(double * x,int n,double * aux)2693 double median_abs(double *x, int n, double *aux)
2694 {
2695 double t;
2696 for(int i=0; i < n; i++) aux[i]=fabs(x[i]);
2697 if ( (n/2) == (double) n / 2 )
2698 t = ( kthplace(aux,n,n/2) + kthplace(aux,n,n/2+1) ) / 2.0 ;
2699 else t = kthplace(aux,n, n/2+1 ) ;
2700 return(t);
2701 }
2702
disp_vec(double * a,int n)2703 void disp_vec(double *a, int n)
2704 {
2705 for(int i=0; i < n; i++) Rprintf("%lf ",a[i]);
2706 Rprintf("\n");
2707 }
2708
disp_veci(int * a,int n)2709 void disp_veci(int *a, int n)
2710 {
2711 for(int i=0; i < n; i++) Rprintf("%d ",a[i]);
2712 Rprintf("\n");
2713 }
2714
disp_mat(double ** a,int n,int m)2715 void disp_mat(double **a, int n, int m)
2716 {
2717 for(int i=0; i < n; i++) {
2718 Rprintf("\n");
2719 for(int j=0; j < m; j++) Rprintf("%10.8f ",a[i][j]);
2720 }
2721 Rprintf("\n");
2722 }
2723
R_find_D_scale(double * rr,double * kkappa,double * ttau,int * llength,double * sscale,double * cc,int * iipsi,int * ttype,double * rel_tol,int * max_k,int * converged)2724 void R_find_D_scale(double *rr, double *kkappa, double *ttau, int *llength,
2725 double *sscale, double *cc, int *iipsi, int *ttype, double *rel_tol,
2726 int *max_k, int *converged)
2727 {
2728 /* compute D_scale using iterative algorithm
2729 type: 1: d1
2730 2: d2
2731 3: dt1
2732 4: dt2
2733 */
2734 *converged = 0;
2735 for (int k=0; k < *max_k; k++) {
2736 double scale = *sscale, tsum1 = 0, tsum2 = 0;
2737 // calculate weights
2738 for(int i=0; i < *llength; i++) {
2739 double a, w = wgt(rr[i] / ttau[i] / scale, cc, *iipsi);
2740 switch(*ttype) {
2741 case 1: // d1
2742 a = rr[i]/ttau[i];
2743 tsum1 += a*a*w;
2744 tsum2 += w; break;
2745 case 2: // d2
2746 a = rr[i]/ttau[i]*w;
2747 tsum1 += a*a;
2748 tsum2 += w*w; break;
2749 default:
2750 case 3: // dt1
2751 tsum1 += rr[i]*rr[i]*w;
2752 tsum2 += w*ttau[i]*ttau[i]; break;
2753 case 4: // dt2
2754 a = rr[i]*w;
2755 tsum1 += a*a;
2756 a = ttau[i]*w;
2757 tsum2 += a*a; break;
2758 };
2759 }
2760
2761 *sscale = sqrt(tsum1 / tsum2 / *kkappa);
2762
2763 // Rprintf("\n type = %d, scale = %10.8f \n", *ttype, *sscale);
2764
2765 if (fabs(scale - *sscale) < *rel_tol * fmax2(*rel_tol, scale)) {
2766 *converged = 1;
2767 break;
2768 }
2769 }
2770 }
2771
2772 /* specialized function calc_fitted */
2773 /* calculates fitted values from simulation output array. */
2774 /* this is used to process simulation output in the */
2775 /* lmrob_simulation vignette */
R_calc_fitted(double * XX,double * bbeta,double * RR,int * nn,int * pp,int * nnrep,int * nnproc,int * nnerr)2776 void R_calc_fitted(double *XX, double *bbeta, double *RR, int *nn, int *pp, int *nnrep,
2777 int *nnproc, int *nnerr)
2778 {
2779 unsigned long A, B, C, D, E;
2780 A = (unsigned long)*nnerr; B = (unsigned long)*nnproc; C = (unsigned long)*nnrep;
2781 D = (unsigned long)*nn; E = (unsigned long)*pp;
2782 // calculate fitted values over errstr, procstr and replicates
2783 for(unsigned long a = 0; a < A; a++) { // errstr
2784 for(unsigned long b = 0; b < B; b++) { // procstr
2785 for(unsigned long c = 0; c < C; c++) { // replicates
2786 // check for NAs
2787 if (!ISNA(bbeta[c + /* 0*C + */ b*C*E + a*B*E*C])) {
2788 for(unsigned long d = 0; d < D; d++) { // observations
2789 RR[d + c*D + b*C*D + a*B*C*D] = 0; // initialize result
2790 for(unsigned long e = 0; e < E; e++) { // predictors
2791 RR[d + c*D + b*C*D + a*B*C*D] += bbeta[c + e*C + b*C*E + a*B*E*C] *
2792 XX[d + e*D + c*E*D + a*C*E*D];
2793 }
2794 }
2795 }
2796 }
2797 }
2798 }
2799 }
2800