[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[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
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/
[Prev in Thread] |
Current Thread |
[Next in Thread] |
- [Octave-bug-tracker] [bug #52751] The gls script for generalized least squares fails when SIGMA is computed.,
Dan Sebald <=