[Octave-bug-tracker] [bug #52751] The gls script for generalized least squares fails when SIGMA is computed.

[Dan Sebald]

URL:
  <http://savannah.gnu.org/bugs/?52751>

                 Summary: The gls script for generalized least squares fails
when SIGMA is computed.
                 Project: GNU Octave
            Submitted by: sebald
            Submitted on: Thu 28 Dec 2017 06:25:47 AM UTC
                Category: Octave Function
                Severity: 3 - Normal
                Priority: 5 - Normal
              Item Group: Unexpected Error
                  Status: None
             Assigned to: None
         Originator Name: 
        Originator Email: 
             Open/Closed: Open
         Discussion Lock: Any
                 Release: dev
        Operating System: Any

    _______________________________________________________

Details:

The following illustrates how the routine gls.m fails when the second output
is desired:


s = [0.5 1.0 1.5];
B = [0.5 0 0.5; 0 0.5 0.5; 0 0 0.5];
X = randn(1000,3) * [1 0 0; 0 2 0; 0 0 3];
O = kron(diag(s)*diag(s), eye(size(X,1)));
E = randn(size(X)) * diag(s);
Y = X*B + E;
[BETA, SIGMA, R] = gls(Y,X,O);
plot(X,Y,'.');
BETA
SIGMA
sqrt(SIGMA)


resulting in


octave:14> [BETA, SIGMA, R] = gls(Y,X,O);
error: gls: operator /: nonconformant arguments (op1 is 1x1, op2 is 1000x3)
error: called from
    gls at line 120 column 9


However, the routine completes if there is only the first output computed,
i.e.,


octave:15> [BETA,~,~] = gls(Y,X,O);


The issue is that the gls.m routine is mistakenly using the variable 'r' for
two things, the rank and the residuals.  Here is a diff hunk to fix it, and I
will post a changeset after thinking about this a bit because I'm wondering
what to do in the case the observation length minus rank is zero ((rx*cy -
rnk) = 0 ... probably isn't very likely).


diff --git a/scripts/statistics/base/gls.m b/scripts/statistics/base/gls.m
--- a/scripts/statistics/base/gls.m
+++ b/scripts/statistics/base/gls.m
@@ -104,9 +104,9 @@ function [beta, v, r] = gls (y, x, o)
   z = o * z;
   y1 = o * reshape (y, ry*cy, 1);
   u = z' * z;
-  r = rank (u);
+  rnk = rank (u);
 
-  if (r == cx*cy)
+  if (rnk == cx*cy)
     b = inv (u) * z' * y1;
   else
     b = pinv (z) * y1;
@@ -117,7 +117,7 @@ function [beta, v, r] = gls (y, x, o)
   if (isargout (2) || isargout (3))
     r = y - x * beta;
     if (isargout (2))
-      v = (reshape (r, ry*cy, 1))' * (o^2) * reshape (r, ry*cy, 1) / (rx*cy -
r);
+      v = (reshape (r, ry*cy, 1))' * (o^2) * reshape (r, ry*cy, 1) / (rx*cy -
rnk);
     endif
   endif






    _______________________________________________________

Reply to this item at:

  <http://savannah.gnu.org/bugs/?52751>

_______________________________________________
  Message sent via/by Savannah
  http://savannah.gnu.org/