bugGNU Octave - Bugs: bug #60348, [octave forge] (statistics)...


bug #60348: [octave forge] (statistics) logistic_regression

Submitter:  Muhali <muhali>
Submitted:  Tue 06 Apr 2021 09:14:52 AM UTC
Category:  Octave Package Severity:  3 - Normal
Priority:  5 - Normal Item Group:  Incorrect Result
Status:  Fixed Assigned to:  None
Originator Name:  Open/Closed:  * Closed
Release:  * 7.2.0 Operating System:  * Any
Fixed Release:  None Planned Release:  None
* Mandatory Fields

Add a New Comment Rich Markup

Jump to the original submission

Tue 27 Sep 2022 02:15:49 PM UTC, comment #9: 

for reference additional discussion on this function has been moved to https://github.com/gnu-octave/statistics/issues/44

Nicholas Jankowski <nrjank>
Group Member
Mon 26 Sep 2022 02:53:38 AM UTC, comment #8: 

Excellent. Sorry I let it go dormant for a month, glad he was able to sort it out.  Closing as fixed.

Nicholas Jankowski <nrjank>
Group Member
Sun 25 Sep 2022 09:36:07 PM UTC, comment #7: 

Andy Penn has fixed the logistic_regression function with PR#43 found at https://github.com/gnu-octave/statistics/pull/43 .
BISTs have been added as well. This bug can be closed.

Andreas Bertsatos <pr0m1th3as>
Wed 17 Aug 2022 10:02:56 PM UTC, comment #6: 

ok. looking through the algorithm a bit more.  some notes:

1 - right at the start the function performs x = -x.  in the optimization routine that evaluates gamma, it expands g = [z,x]*[theta;beta], where z is a unit valued matrix with a column for each discrete response level in Y. z and z1 basically serve as a logical index arrays for low and high values.  so that equation is a row-select matrix creating "theta_i + (-x) * beta" for each level i and value of x. This of course does match the original help text. It hints as to why the algorithm solves for theta to have an opposite sign of beta, but removing that sign flip doesn't restore the standard definition, as x gets used elsewhere in the gradient calcs.  will continue to look and see if it's easy enough to restore to the standard definition.

2 - returned p is not gamma.  it's a level dependent array of the values used in the log loss calculation. essentially for a 2-level (0-1) fit it's an array of:

y = [0 ,   p = [1-px
     1]         px]

or, it's the likelihood of a correct prediction for that input value. 

as you said, the doc states 'P holds estimates for the conditional distribution of Y given X'. my understanding is that that indicates it should be just px, or 'gamma' as it's generally given.

e.g., with my trivial comment #5 example, the result is:

>> [theta beta dev d1 d2 p] = logistic_regression (y=[0 0 0 0 1 1 1 1 1 1]',x=[0:9]')
theta = 98.193
beta = 28.047
dev = 3.2503e-06
d1 =


d2 =

  -1.6251e-06   5.7115e-06
   5.7115e-06  -2.0479e-05

p =


>> gamma = (1-y).*(1-p) + y.*p
gamma =


Nicholas Jankowski <nrjank>
Group Member
Wed 10 Aug 2022 06:02:49 PM UTC, comment #5: 

closed the other logistic regression reports and pointed them here.  additional testing shows that despite what the help text, most reference material, and part of the source show, the output seems to fit -theta + x*beta instead of what was stated earlier.

the theta and beta are supposed to fit the data to a sigmoid function with linear argument. a simple example:

x = [1 2 3 4 5 6 7 8 9]';
y = [0 0 0 0 1 1 1 1 1]';

the above data should produce a sigmoid that is 0 to the left of 4 and rises to 1 to the right of 5, with a smooth transition between.


[theta, beta] = logistic_regression(y,x)
theta = 98.193
beta = 28.047.

(y = theta + beta*x should be a line passing through the 0-to-1 transition. with the data transition between 4 and 5, theta and beta above obviously don't do that.)

the sigmoid function would normally be:

p(x) = 1/(1+exp(-(theta+beta*x))

plotting the four different sign options:

xtest = -10:0.1:10;
p_p_p = 1/(1+exp(-(theta+beta*xtest));
p_p_n = 1/(1+exp(-(theta+beta*xtest));
p_n_p = 1/(1+exp(-(theta+beta*xtest));
p_n_n = 1/(1+exp(-(theta+beta*xtest));


plotting that you see that only the p_n_p function actually matches the expected sigmoid behavior.  So the output from octaves function is producing a constants for a p(x) = 1/(1+exp(-(-theta+beta*x)) model.

I've determined we should just ignore that gamma block at the end, as it's irrelevant to the code. the subfunction really does look like it should follow +theta + beta*x, so not sure yet why it doesn't.

Nicholas Jankowski <nrjank>
Group Member
Mon 08 Aug 2022 11:51:43 PM UTC, comment #4: 

After reading through the code a bit more, i think i can confirm that gamma as used is defined in the logistic_regression_likelihood subfunction to use theta + x*beta.  (E.g, 1/(1+exp(-(theta+x*beta)))

Nicholas Jankowski <nrjank>
Group Member
Mon 08 Aug 2022 08:31:05 PM UTC, comment #3: 

the opening parameter checks are concerning too.

the function definition is:

[theta, beta, dev, dl, d2l, p] = logistic_regression (y, x, print, theta, beta)

then we get to:

  if (nargin < 4)
    beta = zeros (nx, 1);
  if (nargin < 5)
    g = cumsum (sum (z))' ./ my;
    theta = log (g ./ (1 - g));
  tb = [theta; beta];

theta is the 4th argument, beta the 5th.  if I follow this right, if i supply 3 arguments, it will set beta and theta.  if I supply 4 arguments (meaning i've supplied theta), it'll... define theta for me, but not define beta.

verifying this, if you supply theta but not beta, beta is never defined and you get an error at line 119.

[THETA, BETA, DEV, DL, D2L, P] = logistic_regression (Y, X,0,THETA) ;
error: Invalid call to beta.  Correct usage is:

 -- beta (A, B)
error: called from
    print_usage at line 105 column 5
    beta at line 54 column 5
    logistic_regression at line 119 column 6
>> edit logistic_regression

I'll try to come back to this, but IMHO it seems like this function is in need of a complete overhaul by someone more familiar with the statistical process being used.

Nicholas Jankowski <nrjank>
Group Member
Mon 08 Aug 2022 08:21:08 PM UTC, comment #2: 

more concerning now that I look at this function - the help text gives the calling form as: 

[THETA, BETA, DEV, DL, D2L, P] = logistic_regression (Y, X, PRINT,
          THETA, BETA)

but farther down you see:

  The full form is
             = logistic_regression (Y, X, PRINT, THETA, BETA)

The function definition in the code itself indicates P not GAMMA.  But at the end of the code, there is:

  if (nargout == 6)
    if (nx > 0)
      e = ((x beta) ones (1, nz)) + ((y 0 + 1) theta');
      e = (y 0 + 1) theta';
    gamma = diff ([(y 0), (exp (e) ./ (1 + exp (e))), (y 0 + 1)]')';

which makes one think that the function INTENDED to return Gamma, not P as the last parameter.  But in fact gamma as defined in lines 181-188 is never used and the help text gives the definition of return variable P saying it "holds estimates for the conditional distribution of Y given X".

Hence, I'm now questioning everything said in the other reports about what gamma(x) = +/-Theta +/- x*beta should actually look like to match the code, since the gamma = ... lines I based it on are irrelevant to what's actually being returned.  I have no idea why it's even calculating Gamma. I notice more than one author, so there may have been some mixing of intent quite some time ago.

Nicholas Jankowski <nrjank>
Group Member
Mon 08 Aug 2022 08:18:52 PM UTC, comment #1: 

ok, rechecking this against what is thought to be the corrected definition in bug #56118 bug #58318, where it is suggested the formula should have been written:

logit(gamma_i(x)) = theta_i + x * beta

noting from bug #58651, you thought it must work out to be -THETA + X * BETA.  checking both below from the same database gives:

[THETA, BETA, DEV, DL, D2L, P] = logistic_regression (Y, X) ;
P2 = 1 ./ (1 + exp(-(THETA + X * BETA)));

corr(P2, Y)
ans = 0.7252

P2 = 1 ./ (1 + exp(-(-THETA + X * BETA)));
corr(P2, Y)
ans = 0.8790

(and the old definition)
P2 = 1 ./ (1 + exp(-(THETA - X * BETA)));
corr(P2, Y)
ans = -0.8790

(matlab 2022a still gives [-0.8770; 0.8770])

Nicholas Jankowski <nrjank>
Group Member
Tue 06 Apr 2021 09:14:52 AM UTC, original submission:  

The code below shows that the logistic_regression function outputs incorrect probabilities.

pkg load statistics

urlwrite('http://archive.ics.uci.edu/ml/machine-learning-databases/spambase/spambase.data', '/tmp/spambase.data') ;
load /tmp/spambase.data

X = spambase(:,1:end-1) ;
Y = spambase(:,end) ;

[THETA, BETA, DEV, DL, D2L, P] = logistic_regression (Y, X) ;
P2 = 1 ./ (1 + exp(-(THETA - X * BETA))) ;

corr(P, Y)
corr(P2, Y)

Note that according to the docs, "P holds estimates for the conditional distribution of Y given X."

The results are corr(P,Y)=-0.12 and corr(P2,Y)=-0.88. The latter coincides with the results obtained by Matlab, using

[B DEV STATS] = mnrfit(X,Y+1) ;
P = mnrval(B,X) ;
corr(P, Y+1)

Inspecting P it turns out that, in this particular example, one has

P1(1:1813) = 1 - P2(1:1813)


P1(1814:end) = P2(1814:end)

This may or may not be related to similar open bugs of the same function, such as #56118 and #58651.

Muhali <muhali>


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

Attach Files:

No files currently attached


Carbon-Copy List
  • -email is unavailable- added by pr0m1th3as (Posted a comment)
  • -email is unavailable- added by nrjank
  • -email is unavailable- added by nrjank (Posted a comment)
  • -email is unavailable- added by muhali (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.


    Follow 9 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2022-09-26 nrjank StatusConfirmed Fixed
        Open/ClosedOpen Closed
    2022-08-10 nrjank StatusNone Confirmed
        Release6.1.0 7.2.0
        Operating SystemGNU/Linux Any
    2022-08-10 nrjank Dependencies- bugs #52417 is dependent
    2022-08-10 nrjank Dependencies- bugs #56118 is dependent
    2022-08-10 nrjank Dependencies- bugs #58651 is dependent
    2022-08-08 nrjank Carbon-Copy- Added pr0m1th3as

    Back to the top

    Powered by Savane 3.13-04b1.
    Corresponding source code