bugGNU Octave - Bugs: bug #52751, The gls script for generalized...

 
 

bug #52751: The gls script for generalized least squares fails when SIGMA is computed.

Submitter:  Dan Sebald <sebald>
Submitted:  Thu 28 Dec 2017 06:25:47 AM UTC
   
 
Category:  Octave Function Severity:  3 - Normal
Priority:  5 - Normal Item Group:  Unexpected Error or Warning
Status:  None Assigned to:  None
Originator Name:  Open/Closed:  * Open
Release:  * dev Operating System:  * Any
Fixed Release:  None Planned Release:  None
* Mandatory Fields

Add a New Comment Rich Markup
   

Thu 28 Dec 2017 06:25:47 AM UTC, original submission:  

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


Dan Sebald <sebald>

 

(Note: upload size limit is set to 16384 kB, after insertion of the required escape characters.)

Attach Files:
   
   
Comment:
   

No files currently attached

 

Depends on the following items: None found

Items that depend on this one: None found

 

Carbon-Copy List
  • -email is unavailable- added by sebald (Submitted the item)
  •  

    There are 0 votes so far. Votes easily highlight which items people would like to see resolved in priority, independently of the priority of the item set by tracker managers.

    Only group members can vote.

     

    No changes have been made to this item

    Back to the top

    Powered by Savane 3.13-4448.
    Corresponding source code