[Octave-bug-tracker] [bug #54696] statistics' regress function produces division by zero and complex intervals with linear input

[Miguel V. S. Frasson]

URL:
  <https://savannah.gnu.org/bugs/?54696>

                 Summary: statistics' regress function produces division by
zero and complex intervals with linear input
                 Project: GNU Octave
            Submitted by: frasson
            Submitted on: Wed 19 Sep 2018 11:29:22 PM CEST
                Category: Libraries
                Severity: 3 - Normal
                Priority: 5 - Normal
              Item Group: None
                  Status: None
             Assigned to: None
         Originator Name: 
        Originator Email: 
             Open/Closed: Open
         Discussion Lock: Any
                 Release: 4.0.0
        Operating System: GNU/Linux

    _______________________________________________________

Details:

Hi

With an application using regress(), it happened that input data was a column
of ones (all percents = 100%), then a warning with division by zero appeared.

y = A * ones(n,1) ==> division by zero on line 167 of regress.m
y = A * X(:,2) + B ==> division by zero on line 169 of regress.m

for any A and B real numbers.

System:
=======

Xubuntu 16.04.

>> ver statistics
----------------------------------------------------------------------
GNU Octave Version: 4.0.0
GNU Octave License: GNU General Public License
Operating System: Linux 4.15.0-34-generic #37~16.04.1-Ubuntu SMP Tue Aug 28
10:44:06 UTC 2018 x86_64
----------------------------------------------------------------------
Package Name  | Version | Installation directory
--------------+---------+-----------------------
  statistics *|   1.2.4 | /usr/share/octave/packages/statistics-1.2.4 

Reproduzing:
============

1) Save file bug1.m with contents

1;
## bug with constant y

pkg load statistics

rand("seed",1)
A = rand()
n = 5;

X = [ones(n,1),(1:n)'];
y = A * ones(n,1);

[b, bint, r, rint, stats] = regress (y, X)

2) Save file bug2.m with contents

1;
## bug with linear y
pkg load statistics

rand("seed",1)
A = rand()
B = rand()
n = 5;

X = [ones(n,1),(1:n)'];

y = A * X(:,2) + B
[b, bint, r, rint, stats] = regress (y , X)

3) from terminal

octave --eval bug1

4) read from terminal

A =  0.86680
warning: division by zero
warning: called from
    regress at line 167 column 8
    bug at line 13 column 27
b1 =

   8.6680e-01
   5.5511e-17

bint1 =

   8.6680e-01   8.6680e-01
  -8.8732e-17   1.9975e-16

r1 =

   2.2204e-16
   1.1102e-16
   0.0000e+00
   0.0000e+00
   0.0000e+00

rint1 =

   2.2204e-16 - 3.5332e-16i   2.2204e-16 + 3.5332e-16i
  -2.8400e-16 - 0.0000e+00i   5.0605e-16 + 0.0000e+00i
  -4.9967e-16 - 0.0000e+00i   4.9967e-16 + 0.0000e+00i
  -4.6740e-16 - 0.0000e+00i   4.6740e-16 + 0.0000e+00i
  -3.5332e-16 - 0.0000e+00i   3.5332e-16 + 0.0000e+00i

stats1 =

         -Inf  -3.0000e+00   1.0000e+00   2.0543e-32

5) from terminal

octave --eval bug2

6) read from terminal

A =  0.86680
B =  0.91265
y =

   1.7795
   2.6463
   3.5131
   4.3799
   5.2467

warning: division by zero
warning: called from
    regress at line 169 column 7
    bug2 at line 13 column 27
b2 =

   0.91265
   0.86680

bint2 =

   0.91265   0.91265
   0.86680   0.86680

r2 =

   4.4409e-16
   0.0000e+00
  -4.4409e-16
  -8.8818e-16
  -8.8818e-16

rint2 =

  -1.2868e-15 - 0.0000e+00i   2.1750e-15 + 0.0000e+00i
  -2.6440e-15 - 0.0000e+00i   2.6440e-15 + 0.0000e+00i
  -3.0881e-15 - 0.0000e+00i   2.1999e-15 + 0.0000e+00i
  -2.6191e-15 - 0.0000e+00i   8.4274e-16 + 0.0000e+00i
  -8.8818e-16 - 2.9783e-23i  -8.8818e-16 + 2.9783e-23i

stats2 =

   1.00000       Inf   0.00000   0.00000

7) Observe that

* Notice warning of division by zero on lines 167 and 169
* Notice rint1 (intervals) have complex numbers! (probably sqrt of negative
small numbers)

Cause of problem:
=================

Lines from 165 to 174 from regress.m:
 
  if (nargout > 4)

    R2 = 1 - SSE / sum ((y - mean (y)) .^ 2);
#    F = (R2 / (p - 1)) / ((1 - R2) / dof);
    F = dof / (p - 1) / (1 / R2 - 1);
    pval = 1 - fcdf (F, p - 1, dof);

    stats = [R2 F pval v];

  endif

On line 167, denominator
  sum ((y - mean (y)) .^ 2)
is zero with y == A * ones(n,1), for all a real and n integer.

If (p == 1) or (R2 == 1) then it will happend division by zero on line 169.





    _______________________________________________________

File Attachments:


-------------------------------------------------------
Date: Wed 19 Sep 2018 11:29:22 PM CEST  Name: bug1.m  Size: 177B   By: frasson
Input files to reproduce bug
<http://savannah.gnu.org/bugs/download.php?file_id=45056>
-------------------------------------------------------
Date: Wed 19 Sep 2018 11:29:22 PM CEST  Name: bug2.m  Size: 185B   By: frasson
Input files to reproduce bug
<http://savannah.gnu.org/bugs/download.php?file_id=45057>

    _______________________________________________________

Reply to this item at:

  <https://savannah.gnu.org/bugs/?54696>

_______________________________________________
  Message sent via Savannah
  https://savannah.gnu.org/