continuous.stan: deal with (nearly) singular case

stan-dev · Sep 20, 2019 · 5d91237 · 5d91237
1 parent f2d2d7f
commit 5d91237
Showing 1 changed file with 7 additions and 3 deletions.
diff --git a/src/stan_files/continuous.stan b/src/stan_files/continuous.stan
@@ -93,9 +93,13 @@ transformed data {
     matrix[N, has_intercept + K + K_smooth ] X_ = has_intercept ? append_col(rep_vector(1.0, N), 
                                                   (K_smooth > 0 ? append_col(X[1], S) : X[1])) : 
                                                   (K_smooth > 0 ? append_col(X[1], S) : X[1]);
-    XtX = crossprod(X_);
-    OLS = mdivide_left_spd(XtX, X_' * y);
-    SSR = dot_self(y - X_ * OLS);
+    matrix[cols(X_), cols(X_)] R = qr_thin_R(X_);
+    if (tail(diagonal(R), 1)[1] > 1e-16) {
+      matrix[N, cols(R)] Q = qr_thin_Q(X_);
+      XtX = crossprod(X_);
+      OLS = mdivide_right_tri_low(y' * Q, R')';
+      SSR = dot_self(y - X_ * OLS);
+    } else can_do_OLS = 0;
   }
   if (can_do_normalidglm) {
     XS = K_smooth > 0 ? append_col(X[1], S) : X[1];