Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

I am working on a logistic regression problem which requires minimizing a cost function J[{theta0, theta1, theta2}, X, y] to find the optimal value for fitting parameter {theta0, theta1, theta2}.

X is a known M x 3 matrix of real numbers, and y an M-dimensional vector. M is the number of sample data points (approx. 300) in the training data set. The definition for J is

J[θ_, X_, y_] := -(1/Length[y]) (y.Log[h[θ, X]] + (1 - y).Log[1 - h[θ, X]])

where

h[θ_, X_] := Sigmoid[X.θ];
Sigmoid[z_] := 1/(1 + Exp[-z]);

When I try to minimize J with NMinimize and FindMinimum, both fail due to singularity:

NMinimize[J[{θ0, θ1, θ2}, X1, y1], {θ0, θ1, θ2}]
(* NMinimize::nnum: The function value Indeterminate is not a number at {θ0, θ1, θ2} = {0.673558,0.659492,0.0861047}. >> *)

and

FindMinimum[J[{θ0, θ1, θ2}, X1, y1], {θ0, 0}, {θ1, 0}, {θ2, 0}]
(* FindMinimum::nrnum: The function value Indeterminate is not a real number at {θ0, θ1, θ2} = {0.00303682,0.364698,0.342032}. >> *)

For comparison, in MATLAB, using fminunc, it suffers from NaN too

[theta, cost] = fminunc(@(t)(costFunction(t, X, y)), [0;0;0], optimset('MaxIter', 200));
%Warning: Gradient must be provided for trust-region algorithm; using line-search algorithm instead. 
% In fminunc at 365 
%Error using roots (line 28)
%Input to ROOTS must not contain NaN or Inf.
%<snip>

however, setting the 'GradObj' option makes it work:

[theta, cost] = fminunc(@(t)(costFunction(t, X, y)), [0;0;0], optimset('GradObj', 'on', 'MaxIter', 400));    
%Local minimum possible.
%fminunc stopped because the final change in function value relative to 
%its initial value is less than the default value of the function tolerance.    
%<stopping criteria details>

theta    
%theta =    
%  -24.9330
%    0.2044
%    0.1996

cost
% cost = 0.2035

Is there a way to tweak Mathematica to solve this? The notebook showing the details is here.

Edit (follow-up questions)

  1. Is Indeterminate always an unachievable $+\infty$?

  2. Is there a way to do automatic regularization (i.e. adding a regularization term $\lambda || \vec{\beta} ||^2 $ and optimizing $\lambda$) when fitting parameter vector $\vec{\beta}$ for a binary non-linear logistic regression by using some option/setting with LogitModelFit?

  3. The best fit found with NMinimize is J_min as 0.203498 and {\[Beta]0 -> -25.1613, \[Beta]1 -> 0.206232, \[Beta]2 -> 0.201472} as compared to J_min as 0.203506 and [-24.932998 0.204408 0.199618] found with Matlab's fminunc. It seems Mathematica wins. Right?

  4. Using LogitModelFit on the binary nonlinear logistic regression example, it found a decision boundary with a hole. Is there any way to control the complexity/topology of the decision boundary (thinking under-fitting and/or over-fitting)?

Appendix: complete code

The complete code and data is stored in Github. I intend to keep it for as long as Github allows.

Appendix: example data set

(as per comment)

34.62365962451697,78.0246928153624,0
30.28671076822607,43.89499752400101,0
35.84740876993872,72.90219802708364,0
60.18259938620976,86.30855209546826,1
79.0327360507101,75.3443764369103,1
45.08327747668339,56.3163717815305,0
61.10666453684766,96.51142588489624,1
75.02474556738889,46.55401354116538,1
76.09878670226257,87.42056971926803,1
84.43281996120035,43.53339331072109,1
95.86155507093572,38.22527805795094,0
75.01365838958247,30.60326323428011,0
82.30705337399482,76.48196330235604,1
69.36458875970939,97.71869196188608,1
39.53833914367223,76.03681085115882,0
53.9710521485623,89.20735013750205,1
69.07014406283025,52.74046973016765,1
67.94685547711617,46.67857410673128,0
70.66150955499435,92.92713789364831,1
76.97878372747498,47.57596364975532,1
67.37202754570876,42.83843832029179,0
89.67677575072079,65.79936592745237,1
50.534788289883,48.85581152764205,0
34.21206097786789,44.20952859866288,0
77.9240914545704,68.9723599933059,1
62.27101367004632,69.95445795447587,1
80.1901807509566,44.82162893218353,1
93.114388797442,38.80067033713209,0
61.83020602312595,50.25610789244621,0
38.78580379679423,64.99568095539578,0
61.379289447425,72.80788731317097,1
85.40451939411645,57.05198397627122,1
52.10797973193984,63.12762376881715,0
52.04540476831827,69.43286012045222,1
40.23689373545111,71.16774802184875,0
54.63510555424817,52.21388588061123,0
33.91550010906887,98.86943574220611,0
64.17698887494485,80.90806058670817,1
74.78925295941542,41.57341522824434,0
34.1836400264419,75.2377203360134,0
83.90239366249155,56.30804621605327,1
51.54772026906181,46.85629026349976,0
94.44336776917852,65.56892160559052,1
82.36875375713919,40.61825515970618,0
51.04775177128865,45.82270145776001,0
62.22267576120188,52.06099194836679,0
77.19303492601364,70.45820000180959,1
97.77159928000232,86.7278223300282,1
62.07306379667647,96.76882412413983,1
91.56497449807442,88.69629254546599,1
79.94481794066932,74.16311935043758,1
99.2725269292572,60.99903099844988,1
90.54671411399852,43.39060180650027,1
34.52451385320009,60.39634245837173,0
50.2864961189907,49.80453881323059,0
49.58667721632031,59.80895099453265,0
97.64563396007767,68.86157272420604,1
32.57720016809309,95.59854761387875,0
74.24869136721598,69.82457122657193,1
71.79646205863379,78.45356224515052,1
75.3956114656803,85.75993667331619,1
35.28611281526193,47.02051394723416,0
56.25381749711624,39.26147251058019,0
30.05882244669796,49.59297386723685,0
44.66826172480893,66.45008614558913,0
66.56089447242954,41.09209807936973,0
40.45755098375164,97.53518548909936,1
49.07256321908844,51.88321182073966,0
80.27957401466998,92.11606081344084,1
66.74671856944039,60.99139402740988,1
32.72283304060323,43.30717306430063,0
64.0393204150601,78.03168802018232,1
72.34649422579923,96.22759296761404,1
60.45788573918959,73.09499809758037,1
58.84095621726802,75.85844831279042,1
99.82785779692128,72.36925193383885,1
47.26426910848174,88.47586499559782,1
50.45815980285988,75.80985952982456,1
60.45555629271532,42.50840943572217,0
82.22666157785568,42.71987853716458,0
88.9138964166533,69.80378889835472,1
94.83450672430196,45.69430680250754,1
67.31925746917527,66.58935317747915,1
57.23870631569862,59.51428198012956,1
80.36675600171273,90.96014789746954,1
68.46852178591112,85.59430710452014,1
42.0754545384731,78.84478600148043,0
75.47770200533905,90.42453899753964,1
78.63542434898018,96.64742716885644,1
52.34800398794107,60.76950525602592,0
94.09433112516793,77.15910509073893,1
90.44855097096364,87.50879176484702,1
55.48216114069585,35.57070347228866,0
74.49269241843041,84.84513684930135,1
89.84580670720979,45.35828361091658,1
83.48916274498238,48.38028579728175,1
42.2617008099817,87.10385094025457,1
99.31500880510394,68.77540947206617,1
55.34001756003703,64.9319380069486,1
74.77589300092767,89.52981289513276,1
share|improve this question
1  
can you add some example data? the code you gave doesn't really work without the data (it gives a different error and can be confusing). It's a bad idea to rely on externally-linked files as they'll eventually disappear, and then this question is impossible to understand. –  acl Oct 27 '12 at 23:05
    
@acl: Added data set. –  qazwsx Oct 29 '12 at 1:16
1  
You should probably reformulate your last follow up question as a new question. –  Andy Ross Oct 29 '12 at 2:31
    
this is a problem from Coursera, right? –  magma Nov 5 '12 at 1:51

1 Answer 1

up vote 10 down vote accepted

You could try something like the following.

bounded[Indeterminate] := $MaxMachineNumber;
bounded[x_?NumericQ] := x

This gives a very large number instead of Indeterminate so NMinimize keeps searching.

NMinimize[bounded[J[{θ0, θ1, θ2}, X1, y1]], {θ0, θ1, θ2}]    
(* {0.203498, {θ0 -> -25.1613, θ1 -> 0.206231, θ2 -> 0.201471}} *)

Incidentally, is there any reason you don't use LogitModelFit?

LogitModelFit[Transpose[Transpose[X1]~Join~{y1}],
    {θ0, θ1, θ2}, {θ0, θ1, θ2}, IncludeConstantBasis -> False]["BestFitParameters"]

(*{-25.1613, 0.206232, 0.201472}*)

Edit:

In response to the edits. The only way for Indeterminate to arise beyond the length of y being zero is for the sigmoid function to return exact zero or 1. This should never happen with real input.

Reduce[Sigmoid[z] == 0 || Sigmoid[z] == 1, z, Reals]

(*False *)

LogitModelFit cannot perform regularization as of version 8. You will have to do this using your NMinimize approach.

share|improve this answer
    
Could you explain a little that why when Sigmoid returns exactly 0 or 1, there will then be a Indeterminate? Is it that y_i - Sigmoid[...] is somehow used as a denominator in the computation? –  qazwsx Oct 29 '12 at 1:50
    
Because your cost function contains Log[h[...]] and Log[1-h[...]] If h[...] gives 0 or 1 one of the logs would return Indeterminate inside NMinimize since -Infinity is not a real number. –  Andy Ross Oct 29 '12 at 2:26
    
I'm still not exactly sure why converting Indeterminate to $MaxMachineNumber is safe that no $-\infty$ will not be mistakenly converted. –  qazwsx Oct 29 '12 at 3:14
    
Because over the machine reals that NMinmize will try there is no instance that gives -Infinity. –  Andy Ross Oct 29 '12 at 4:09

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.