Why does FindFit work and NonlinearModelFit does not?

Question

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]

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]

Answer 1

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.}}

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.