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

 
 

bug #56118: [octave forge] (statistics) logistic_regression gives incorrect result

Submitter:  Muhali <muhali>
Submitted:  Thu 11 Apr 2019 03:08:21 PM UTC
   
 
Category:  Octave Package Severity:  3 - Normal
Priority:  5 - Normal Item Group:  Incorrect Result
Status:  Duplicate 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

Wed 10 Aug 2022 05:55:24 PM UTC, comment #7: 

ok. after some testing it appears that despite the help text and some parts of the code, the output appears to fit -theta + beta * x.  I'm going to be close this report and shifting discussion and corrections over to bug #60348, since some discussion already over there and all of these logistic_regression reports are related.

Nicholas Jankowski <nrjank>
Group Member
Mon 08 Aug 2022 07:11:25 PM UTC, comment #6: 

thanks for the pointers.  looking at how they all seem to hinge on the help text definition of what's happening under the hood, we probably could have combined these. will take a look.

@muhali, could you also look at those as well in light of the recent changes. I started looking but am not sure that definition clarification alone is enough to fix the problems

Nicholas Jankowski <nrjank>
Group Member
Mon 08 Aug 2022 06:47:59 PM UTC, comment #5: 


comment #3:

> So, I'll recommend that change to the help text. Given that do you see any other issues that would prevent closing this report?


looks good to me now. Thx.

Muhali <muhali>
Mon 08 Aug 2022 05:27:06 PM UTC, comment #4: 

New patch pushed at https://github.com/gnu-octave/statistics/commit/ba48f6b7f85e91bdf257037b0b1a5759bee1aec2

There are also bugs #58651 and #60348 related to the logistic_regression function. The former is about documentation and the latter about incorrect results. It would be nice if you could see these through as well.

In the meantime I am trying to implement mnrfit (related to bug #58318, not a bug but rather a missing function) and mnrval, but they need quite some work yet.

Andreas Bertsatos <pr0m1th3as>
Mon 08 Aug 2022 05:06:33 PM UTC, comment #3: 

I was hoping you might be familiar enough with the correct/standard definition and could verify whether the problem is the output or the text.
From a bit of reading up on the topic, my suspicion is that the function is correct.  The help text says logit(gamma(x) = theta-beta x. The same part of the code mentioned in the other bug report that outputs gamma defines it as as theta+x*beta:


line 181:

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


using that, I tried just changing your exponential arguments in the numerator and denominator, but the output was still nonsensical. (noticed that Pfun of p2 was 1 then went to nan at higher values)

recentering around zero, and rerunning your code as:

x = (-10:10)' ; y = [zeros(11, 1) ; ones(10,1)] ;
Pfun = @(p, x) exp(p(2) * x + p(1)) ./ (1 + exp(p(2) * x + p(1))) ; Pin = [0 0]' ;

p1 = nonlin_curvefit (Pfun, Pin, x, y) ;
[p2(1) p2(2)] = logistic_regression(y, x) ;

subplot(2,1,1) ;
plot(x, y, "x", x, Pfun(p1, x), "-") ;
set(gca, "box", "off")
subplot(2,1,2) ;
plot(x, y, "x", x, Pfun(p2, x), "-") ;
set(gca, "box", "off")


it's not perfect but it seems to at least behave as expected.

So, I'll recommend that change to the help text. Given that do you see any other issues that would prevent closing this report?

(cc package maintainer and attaching a simple doc patch to make that help text change as well as add variable tags to the texinfo in the formulas)

(file #53544)

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

somehow, the function call has changed. Now one needs

Pfun = @(p, x) ... :


pkg load optim
x = (1:20)' ; y = [zeros(10, 1) ; ones(10,1)] ;
Pfun = @(p, x) exp(p(1) - p(2) * x) ./ (1 + exp(p(1) - p(2) * x)) ; Pin = [0 0]' ;

p1 = nonlin_curvefit (Pfun, Pin, x, y) ;
[p2(1) p2(2)] = logistic_regression(y, x) ;

subplot(2,1,1) ;
plot(x, y, "x", x, Pfun(p1, x), "-") ;
set(gca, "box", "off")
subplot(2,1,2) ;
plot(x, y, "x", x, Pfun(p2, x), "-") ;
set(gca, "box", "off")


Whether the function or the documentation is incorrect is just a question of how the regression equation is defined. They only need to be consistent.

Muhali <muhali>
Sat 06 Aug 2022 06:02:14 PM UTC, comment #1: 

trying to run comment #0 code produces an error:


>> pkg load statistics
>> x = (1:20)' ; y = [zeros(10, 1) ; ones(10,1)] ;
>> Pfun = @(p) exp(p(1) - p(2) * x) ./ (1 + exp(p(1) - p(2) * x)) ; Pin = [0 0]' ;
>> p1 = nonlin_curvefit (Pfun, Pin, x, y) ;
error: @<anonymous>: function called with too many inputs
error: called from
    @<anonymous>
    __maybe_limit_arg_count__>@<anonymous> at line 30 column 26
    nonlin_curvefit>@<anonymous> at line 87 column 24
    __nonlin_residmin__ at line 310 column 9
    nonlin_curvefit at line 86 column 22


preventing p1 from being defined. I don't know if that matters at all for the issue, as it does seem that the p2 output is inverted and the correct result is produced by Pfun(-p2) according to the p2 plot.

not clear enough from comment #0 if they think the problem is the  function description or the function itself.  Guessing from the description of the gamma function definition there and at wikipedia, where it is given as gamma = B0 + B1*x (where B0 is here called theta and B1 is beta), it seems that in addition to the ordering change in bug #52417, that perhaps a sign change is needed too.  again, the code shows that "x*beta + theta" is being performed, not "theta - beta * x".  Perhaps this accounts for the sign change expected by the user on the beta term?

if that's the case, this should be a fairly easy doc fix. if the code needs to be fixed, then we'll need more info.

Nicholas Jankowski <nrjank>
Group Member
Thu 11 Apr 2019 03:08:21 PM UTC, original submission:  

Applying a logistic regression in the following way


pkg load optim
x = (1:20)' ; y = [zeros(10, 1) ; ones(10,1)] ;
Pfun = @(p) exp(p(1) - p(2) * x) ./ (1 + exp(p(1) - p(2) * x)) ; Pin = [0 0]' ;

p1 = nonlin_curvefit (Pfun, Pin, x, y) ;
[p2(1) p2(2)] = logistic_regression(y, x) ;

subplot(2,1,1) ;
plot(x, y, "x", x, Pfun(p1), "-") ;
set(gca, "box", "off")
subplot(2,1,2) ;
plot(x, y, "x", x, Pfun(p2), "-") ;
set(gca, "box", "off")


one sees that the result cannot be correct. The reason is a mismatch between the function description in the help text and what is actually programmed.

Using Pfun(-p2) gives the correct result.

Muhali <muhali>

 

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

Attach Files:
   
   
Comment:
   

Attached Files

 

Digest:
   bug dependencies.

Items that depend on this one: None found

 

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 14 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2022-08-10 nrjank StatusPatch Submitted Duplicate
        Open/ClosedOpen Closed
        Dependencies- Depends on bugs #60348
    2022-08-08 nrjank StatusConfirmed Patch Submitted
    2022-08-08 nrjank Attached File- Added doc-logistic_regression-formula-fix-pt2-bug-56118.patch, #53544
    2022-08-06 nrjank Carbon-Copy- Added pr0m1th3as
    2022-08-06 nrjank Carbon-CopyRemoved pr0m1th3us -
    2022-08-06 nrjank Carbon-Copy- Added pr0m1th3us
    2022-08-06 nrjank Item GroupNone Incorrect Result
        StatusNone Confirmed
        Release5.1.0 7.2.0
        Operating SystemGNU/Linux Any
    2019-04-11 mtmiller Summary[octave forge] (optim) logistic_regression gives incorrect result [octave forge] (statistics) logistic_regression gives incorrect result
    2019-04-11 mtmiller Summary(optim) logistic_regression gives incorrect result [octave forge] (optim) logistic_regression gives incorrect result

    Back to the top

    Powered by Savane 3.13-4448.
    Corresponding source code