octave-bug-tracker
[Top][All Lists]

## [Octave-bug-tracker] [bug #52751] The gls script for generalized least s

 From: Dan Sebald Subject: [Octave-bug-tracker] [bug #52751] The gls script for generalized least squares fails when SIGMA is computed. Date: Thu, 28 Dec 2017 01:25:48 -0500 (EST) User-agent: Mozilla/5.0 (X11; Ubuntu; Linux x86_64; rv:55.0) Gecko/20100101 Firefox/55.0

```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
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

_______________________________________________________

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

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

```