2
$\begingroup$

I have a somewhat complicated function that I have defined.

myfunction[t_, qi_, di_, b_, dt_] :=
 Module[{din, dtn, tdt, qtdt},
  din = If[b == 0, -Log[1 - di], ((1 - di)^-b - 1)/b];
  If[b == 0, (qi - qi*Exp[-din*t/12])/din,
    dtn = -Log[1 - dt];
    tdt = (din/dtn - 1)*12/(b*din);
    If[t <= tdt,
     If[b == 1, qi/din*Log[(1 + b*din*t/12)^(1/b)], qi^b/((b - 1)*din)*((qi*(1 + b*din*t/12)^(-1/b))^(1 - b) - qi^(1 - b))],
     qtdt = qi*(1 + b*din*tdt/12)^(-1/b);
     If[b == 1, qi/din*Log[(1 + b*din*tdt/12)^(1/b)] + (qtdt - qtdt*Exp[-dtn*(t - tdt)/12])/dtn,
      qi^b/((b - 1)*din)*((qi*(1 + b*din*tdt/12)^(-1/b))^(1 - b) - qi^(1 - b)) + (qtdt - qtdt*Exp[-dtn*(t - tdt)/12])/ dtn]]]*12]

And here is sample data to fit to:

data= {{1, 1116}, {2, 5116}, {3, 8168}, {4, 10362}, {5, 13381}, {6,   15989}, {7, 17606}, {8, 19233}, {9, 20738}, {10, 22006}, {11,   23180}, {12, 24087}, {13, 25547}, {14, 29835}, {15, 33275}, {16,   35482}, {17, 37634}, {18, 39340}, {19, 41056}, {20, 43071}, {21,   44697}, {22, 46414}, {23, 47917}, {24, 49870}, {25, 51350}, {26,   53339}, {27, 55418}, {28, 57095}, {29, 58724}, {30, 60600}, {31,   61704}, {32, 63050}, {33, 64510}, {34, 66156}, {35, 67429}, {36,  68764}, {37, 69902}, {38, 71286}, {39, 72714}, {40, 73667}, {41,   74864}, {42, 75072}, {43, 75344}, {44, 75344}, {45, 76453}, {46,   76964}, {47, 77396}, {48, 78485}, {49, 80406}, {50, 81696}, {51,   82463}, {52, 84066}, {53, 85591}, {54, 86119}, {55, 86554}, {56,   87947}, {57, 88715}, {58, 89661}, {59, 90868}, {60, 91742}, {61,   92221}, {62, 93451}, {63, 95129}, {64, 98020}, {65, 119405}}

Now the function must have constraints and I have found that FindFit works a lot better if I give it starting values, so I use the following code and get the the following answer:

FindFit[data, {myfunction[t, qi, di, b, 0.1], {di >= 0, 0 <= b <= 3}}, {{qi, 4000}, {di, 0.5}, {b, 2}}, t]


{qi -> 2664.49, di -> 0.242422, b -> 0.770906}

I would like to have all of the functionality of NonlinearModelFit (residuals, confidence intervals, and all that good stuff), but I type in the same arguments and it spits out something unintelligible:

NonlinearModelFit[data, {myfunction[t, qi, di, b, 0.1], {di >= 0, 0 <= b <= 3}}, {{qi, 4000}, {di, 0.5}, {b, 2}}, t]

enter image description here

Why can I not get NonlinearModelFit to work? when FindFit works just fine?

UPDATE Based on somebody's comment, it was suggested to try to write my code with Piecewise. So I have created that below. However, with this function, neither fit function appears to work.

myfunctionpiecewise[t_, qi_, di_, b_, dt_] :=
Module[{din, dtn, tdt, qtdt},
      din = Piecewise[{{-Log[1 - di], b == 0}, {((1 - di)^-b - 1)/b, 
          b > 0}}];
      If[b != 0,
       dtn = -Log[1 - dt];
       tdt = (din/dtn - 1)*12/(b*din);
       qtdt = qi*(1 + b*din*tdt/12)^(-1/b)];
      Piecewise[{{(qi - qi*Exp[-din*t/12])/din, b == 0},
         {Piecewise[{{qi/din*Log[(1 + b*din*t/12)^(1/b)], 
             t <= tdt}, {qi/din*
               Log[(1 + b*din*tdt/12)^(1/b)] + (qtdt - 
                 qtdt*Exp[-dtn*(t - tdt)/12])/dtn, t > tdt}}], b == 1},
         {Piecewise[{{qi^
                b/((b - 1)*din)*((qi*(1 + b*din*t/12)^(-1/b))^(1 - b) - 
                qi^(1 - b)), 
             t <= tdt}, {qi^
                 b/((b - 1)*din)*((qi*(1 + b*din*tdt/12)^(-1/b))^(1 - b) -
                  qi^(1 - b)) + (qtdt - qtdt*Exp[-dtn*(t - tdt)/12])/dtn, 
             t > tdt}}], b != 1 || b > 0}}]*12]
$\endgroup$
  • $\begingroup$ Have you considered reformulating your model to use Piecewise[]? $\endgroup$ – J. M. is away Aug 17 '15 at 15:49
  • $\begingroup$ I could certainly do that. Would it help? $\endgroup$ – Jud Aug 17 '15 at 16:43
  • $\begingroup$ Jud, I can't reproduce your FindFit result. NonlinearModelFit and FindFit behave similarly on my system, in that neither one of them is able to converge on a good fit with your starting conditions. I am using MMA 10.2 on Win7-64. What version / OS are you using? $\endgroup$ – MarcoB Aug 17 '15 at 19:13
  • $\begingroup$ There are things to do after your question is answered. It's a good idea to stay vigilant for some time, better approaches may come later improving over previous replies. Experienced users may point alternatives, caveats or limitations. New users should test answers before voting and wait 24 hours before accepting the best one. Participation is essential for the site, please come back to do your part tomorrow $\endgroup$ – rhermans Feb 4 '16 at 17:03
5
$\begingroup$

By removing the restrictions in the parameters things work fine (and much faster than leaving them in). In this case the resulting solution satisfies the restrictions so there's not reason to include them. However, because of the way you've written your function there are some warning messages (which I'm sure that very experienced Mathematica users can fix).

Here is the code I used:

nlm = NonlinearModelFit[data, myfunction[t, qi, di, b, 0.1], {{qi, 4000}, {di, 0.5}, {b, 2}}, t];
sol = nlm["BestFitParameters"]
nlm["CorrelationMatrix"]
Show[ListPlot[data], Plot[myfunction[t, qi /. sol, di /. sol, b /. sol, 0.1], {t, 0, 70}], Frame -> True]

with the following results

{qi -> 3745.17, di -> 0.5, b -> 2.91276}
{{1., 0.99956, 0.944353}, {0.99956, 1., 0.951849}, {0.944353, 0.951849, 1.}}

Data and regression curve

Note that the estimates of the parameters are far away from your estimates. This is because your model is close to having really only two parameters rather than three as can be seen from the large correlations in the parameter estimates. In this case different starting values can result in wildly different parameter estimates. However, the predictions will be nearly identical.

The two warnings I get are both Set::write: Tag Slot in #1 is Protected. >>. Again, I suspect that this is because of the way you've written your function and how that function interacts with NonlinearModelFit. The indicators of this issue is that the following command

nlm[5]

results in

12 If[5 <= (
   4.11981 (-1 + 
      9.49122 If[b == 0, -Log[1 - di], ((1 - di)^-b - 1)/b]))/
   If[b == 0, -Log[1 - di], ((1 - di)^-b - 1)/b], 
  If[2.91276 == 1, (
   3745.17 Log[(1 + (2.91276 din$36788 5)/12)^(1/2.91276)])/
   din$36788, (
   3745.17^2.91276 ((3745.17 (1 + (2.91276 din$36788 5)/12)^(-1/
         2.91276))^(1 - 2.91276) - 3745.17^(
      1 - 2.91276)))/((2.91276 - 1) din$36788)], 
  qtdt$36788 = 
   3745.17 (1 + (2.91276 din$36788 tdt$36788)/12)^(-1/2.91276); 
  If[2.91276 == 1, (
    3745.17 Log[(1 + (2.91276 din$36788 tdt$36788)/12)^(1/2.91276)])/
    din$36788 + (
    qtdt$36788 - qtdt$36788 Exp[1/12 (-dtn$36788) (5 - tdt$36788)])/
    dtn$36788, (
    3745.17^2.91276 ((3745.17 (1 + (2.91276 din$36788 tdt$36788)/
           12)^(-1/2.91276))^(1 - 2.91276) - 3745.17^(
       1 - 2.91276)))/((2.91276 - 1) din$36788) + (
    qtdt$36788 - qtdt$36788 Exp[1/12 (-dtn$36788) (5 - tdt$36788)])/
    dtn$36788]]

rather than a number for when t = 5.

$\endgroup$
  • $\begingroup$ Note also that the value of the di parameter is essentially unchanged in the optimization. If you use a different starting value, the same value will be returned by the fit. @Jud should really rethink his model thoroughly. $\endgroup$ – MarcoB Aug 17 '15 at 19:50
  • $\begingroup$ @MarcoB. I agree. And the extreme serial correlation in the residuals should also cause some concern but I've noticed that such deviations in residual behavior (or local lack of fit) tends to be tolerated or ignored at this site. $\endgroup$ – JimB Aug 17 '15 at 19:55
  • 1
    $\begingroup$ Thank you all for the responses. A couple things to note: the constraints are important. This is an engineering equation and the theoretical physics fail if left to regress to numbers outside of those bounds. I also realize now, that there was not a problem all along with the regression with the output I got before, I still could get the parameters the way that @JimB has described. I had tried something similar to nlm[5], and because that failed, I thought there was a problem. But the regression still worked. I feel like an idiot now for not trying that. Thanks again. $\endgroup$ – Jud Aug 17 '15 at 22:33
  • $\begingroup$ @jud Change the definition of the function to myfunction[t_?NumericQ, qi_?NumericQ, di_?NumericQ, b_?NumericQ, dt_?NumericQ] will suppress the Warning. (Remember to Clear@myfunction first. ) This issue has been addressed for several times before, just search in this site for more details. $\endgroup$ – xzczd Oct 11 '15 at 2:52

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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