1
$\begingroup$

I have defined the following function in Mathematica:

$\int_0^t[\int_0^a e^{-\frac{r^2+r'^2}{4Diff(t-t')}}I_o (\frac{rr'}{2Diff(t-t')})r'^2 dr']\times \frac{e^{-\alpha z + \alpha^2 Diff(t-t')}}{Diff(t-t')} Erfc(\frac{2Diff \alpha(t-t')-z}{\sqrt{4Diff(t-t')}})dt'$

TempDist[r_, z_, t_, a_, Diff_, \[Alpha]_] := NIntegrate[(ro^2) Exp[-(ro^2 + r^2)/(4 Diff (t - ti))] BesselI[0, ro r/(2 Diff (t - ti))] Exp[-\[Alpha]*z + (\[Alpha]^2)*Diff (t - ti)]*Erfc[(2 Diff \[Alpha] (t - ti) - z)/(Sqrt[4 Diff (t - ti)])]/(Diff (t - ti)), {ro, 0, a}, {ti, 0, t}, PrecisionGoal -> 4, MaxRecursion -> 10];

I am trying to fit a dataset to this function (1000 points). The independent variables are $r$ and $t$. The values of $z$ and $\alpha$ are already known, so the only parameters I'm trying to get are $a$ and $Diff$.

The problem is that when I ask Mathematica to fit:

nlm = NonlinearModelFit[data, 4.15*TempDist[r, 2, t, a, Diff, 0.023], {{Diff, 0.125}, {a, 4.0}}, {r, t}]

The following messages appear:

NIntegrate::izero: Integral and error estimates are 0 on all integration subregions. Try increasing the value of the MinRecursion option. If value of integral may be 0, specify a finite value for the AccuracyGoal option. >>

General::stop: Further output of NIntegrate::izero will be suppressed during this calculation. >>

And the calculation never finishes. Am I missing something on the way the integral is being computed?

Here is the dataset I am trying to fit:

{{0,0,0.788942},{0,6,0.},{0,12,0.},{0,18,0.788942},{0,24,1.05192},{0,30,1.05192},{0,36,0.920433},{0,42,2.62981},{0,48,2.10385},{0,54,2.36683},{0,60,1.70937},{0,66,3.55024},{0,72,2.89279},{0,78,2.89279},{0,84,3.28726},{0,90,3.68173},{0,96,3.15577},{0,102,3.28726},{0,108,3.41875},{0,114,2.89279},{0,120,2.49832},{0,126,3.41875},{0,132,3.15577},{0,138,2.89279},{0,144,3.41875},{0,150,3.02428},{0,156,3.41875},{0,162,3.94471},{0,168,3.02428},{0,174,3.15577},{0,180,3.81322},{0,186,3.68173},{0,192,3.81322},{0,198,3.02428},{0,204,3.41875},{0,210,3.15577},{0,216,3.28726},{0,222,3.68173},{0,228,3.55024},{0,234,4.20769},{0,240,3.15577},{0,246,4.20769},{0,252,3.94471},{0,258,4.0762},{0,264,4.60216},{0,270,3.28726},{0,276,3.94471},{0,282,3.94471},{0,288,4.33918},{0,294,3.68173},{0,300,3.41875},{0,306,3.81322},{0,312,4.73365},{0,318,3.41875},{0,324,3.41875},{0,330,3.81322},{0,336,4.47067},{0,342,4.20769},{0,348,3.68173},{0,354,3.94471},{0,360,3.81322},{0,366,4.73365},{0,372,4.33918},{0,378,4.99663},{0,384,4.60216},{0,390,4.73365},{0,396,4.33918},{0,402,4.47067},{0,408,4.0762},{0,414,5.25962},{0,420,4.20769},{0,426,3.28726},{0,432,4.86514},{0,438,4.60216},{0,444,4.33918},{0,450,4.33918},{0,456,3.94471},{0,462,3.94471},{0,468,4.20769},{0,474,3.68173},{0,480,4.86514},{0,486,3.28726},{0,492,5.12813},{0,498,4.0762},{0,504,5.78558},{0,510,5.65409},{0,516,4.47067},{0,522,5.12813},{0,528,5.5226},{0,534,4.60216},{0,540,4.33918},{0,546,4.73365},{0,552,5.12813},{0,558,5.12813},{0,564,4.73365},{0,570,4.60216},{0,576,4.20769},{0,582,4.33918},{0,588,4.60216},{0,594,4.0762},{0,600,4.73365},{0,606,6.18005},{0,612,5.25962},{0,618,4.73365},{0,624,5.39111},{0,630,4.73365},{0,636,5.25962},{0,642,4.33918},{0,648,5.39111},{0,654,4.47067},{0,660,3.68173},{0,666,3.94471},{2,0,1.16512},{2,6,0.404215},{2,12,0.020846},{2,18,0.824866},{2,24,1.14326},{2,30,1.48389},{2,36,1.25155},{2,42,2.55875},{2,48,1.82298},{2,54,2.03449},{2,60,1.3055},{2,66,2.92561},{2,72,2.79616},{2,78,2.98746},{2,84,3.009},{2,90,3.63993},{2,96,3.47814},{2,102,3.22542},{2,108,3.74531},{2,114,3.56434},{2,120,2.98114},{2,126,3.84612},{2,132,3.73338},{2,138,3.7685},{2,144,3.78359},{2,150,3.49877},{2,156,3.59027},{2,162,3.81559},{2,168,3.89999},{2,174,3.62181},{2,180,3.73346},{2,186,3.93535},{2,192,3.82107},{2,198,3.62158},{2,204,3.65875},{2,210,3.88814},{2,216,3.61604},{2,222,4.00039},{2,228,4.12578},{2,234,4.2247},{2,240,3.53187},{2,246,4.31957},{2,252,4.19584},{2,258,4.49837},{2,264,4.33971},{2,270,4.36204},{2,276,3.96769},{2,282,4.4017},{2,288,4.58119},{2,294,4.18808},{2,300,4.15935},{2,306,4.10752},{2,312,4.76225},{2,318,4.11924},{2,324,4.01998},{2,330,4.05084},{2,336,4.53771},{2,342,4.71143},{2,348,4.36024},{2,354,3.88906},{2,360,3.93614},{2,366,4.64211},{2,372,4.35517},{2,378,4.84314},{2,384,4.86164},{2,390,4.69317},{2,396,4.47494},{2,402,4.52927},{2,408,4.50602},{2,414,5.1862},{2,420,4.48845},{2,426,4.11992},{2,432,4.87873},{2,438,4.41536},{2,444,4.33737},{2,450,4.38906},{2,456,4.34917},{2,462,4.20792},{2,468,4.69261},{2,474,4.02874},{2,480,5.06588},{2,486,4.14234},{2,492,5.32931},{2,498,4.41691},{2,504,5.69617},{2,510,5.58655},{2,516,4.54162},{2,522,4.52879},{2,528,4.99459},{2,534,4.61873},{2,540,4.32658},{2,546,4.94736},{2,552,4.49702},{2,558,4.5138},{2,564,4.30028},{2,570,4.03123},{2,576,4.11642},{2,582,3.96611},{2,588,4.40981},{2,594,4.07012},{2,600,4.47845},{2,606,5.30463},{2,612,5.37677},{2,618,4.50874},{2,624,5.18174},{2,630,5.29879},{2,636,4.49508},{2,642,4.34152},{2,648,5.21204},{2,654,4.17029},{2,660,4.52553},{2,666,4.87164},{4,0,0.549135},{4,6,0.},{4,12,0.},{4,18,0.339802},{4,24,0.416546},{4,30,0.511533},{4,36,0.507435},{4,42,1.71517},{4,48,1.32578},{4,54,1.36941},{4,60,0.831852},{4,66,2.01597},{4,72,1.59388},{4,78,1.66633},{4,84,1.95171},{4,90,2.26955},{4,96,1.98952},{4,102,2.20216},{4,108,2.45161},{4,114,2.02402},{4,120,1.44727},{4,126,2.45963},{4,132,1.97976},{4,138,2.40216},{4,144,2.29553},{4,150,1.84934},{4,156,1.97598},{4,162,2.40894},{4,168,2.16071},{4,174,2.12593},{4,180,2.21221},{4,186,2.19838},{4,192,2.00368},{4,198,2.30499},{4,204,2.22338},{4,210,2.03393},{4,216,1.90005},{4,222,2.30767},{4,228,2.35241},{4,234,2.58658},{4,240,2.04671},{4,246,2.57541},{4,252,3.12554},{4,258,2.9447},{4,264,3.14002},{4,270,2.55241},{4,276,2.96393},{4,282,2.66818},{4,288,2.83905},{4,294,2.67855},{4,300,2.59335},{4,306,2.50766},{4,312,3.00962},{4,318,2.77949},{4,324,2.45352},{4,330,2.13979},{4,336,2.98304},{4,342,2.92944},{4,348,2.6658},{4,354,2.63274},{4,360,2.3105},{4,366,3.02869},{4,372,2.96175},{4,378,3.1147},{4,384,3.00837},{4,390,2.59277},{4,396,2.81954},{4,402,2.74473},{4,408,2.86459},{4,414,3.77413},{4,420,3.01357},{4,426,2.59038},{4,432,3.49491},{4,438,2.94494},{4,444,2.93704},{4,450,2.93437},{4,456,3.00539},{4,462,2.72976},{4,468,3.18475},{4,474,2.7241},{4,480,3.15029},{4,486,2.52815},{4,492,3.52055},{4,498,2.70347},{4,504,4.0036},{4,510,3.72909},{4,516,3.14209},{4,522,2.97925},{4,528,3.20306},{4,534,2.94306},{4,540,2.75863},{4,546,3.16478},{4,552,3.10525},{4,558,2.96339},{4,564,2.8054},{4,570,2.6101},{4,576,2.48377},{4,582,2.34236},{4,588,2.98774},{4,594,2.40894},{4,600,2.76505},{4,606,3.71554},{4,612,3.62671},{4,618,3.04395},{4,624,3.5903},{4,630,3.27014},{4,636,3.15057},{4,642,2.67107},{4,648,3.2207},{4,654,2.63072},{4,660,2.78677},{4,666,2.92504},{6,0,1.21744},{6,6,0.374256},{6,12,0.167872},{6,18,0.603895},{6,24,1.10018},{6,30,1.08806},{6,36,1.0809},{6,42,1.90261},{6,48,1.56218},{6,54,1.53414},{6,60,1.0403},{6,66,1.90379},{6,72,1.84383},{6,78,1.76297},{6,84,1.81313},{6,90,2.17598},{6,96,2.09796},{6,102,2.041},{6,108,2.26989},{6,114,2.09157},{6,120,1.70227},{6,126,2.30022},{6,132,2.20457},{6,138,2.19272},{6,144,2.19187},{6,150,1.90349},{6,156,2.19106},{6,162,2.13228},{6,168,1.97139},{6,174,2.11059},{6,180,2.13849},{6,186,2.29438},{6,192,1.97575},{6,198,2.24628},{6,204,2.30613},{6,210,2.04594},{6,216,1.98422},{6,222,2.39187},{6,228,2.32553},{6,234,2.52546},{6,240,2.16609},{6,246,2.56854},{6,252,2.55461},{6,258,2.65491},{6,264,2.56333},{6,270,2.40599},{6,276,2.33022},{6,282,2.27972},{6,288,2.73393},{6,294,2.42457},{6,300,2.374},{6,306,2.46517},{6,312,2.67852},{6,318,2.41419},{6,324,2.43558},{6,330,2.29213},{6,336,2.46355},{6,342,2.58413},{6,348,2.47835},{6,354,2.17946},{6,360,2.34481},{6,366,2.8032},{6,372,2.52823},{6,378,2.79229},{6,384,2.58217},{6,390,2.6091},{6,396,2.42302},{6,402,2.73182},{6,408,2.45298},{6,414,3.1597},{6,420,2.66403},{6,426,2.39254},{6,432,2.87258},{6,438,2.61748},{6,444,2.70841},{6,450,2.55228},{6,456,2.54762},{6,462,2.49003},{6,468,2.94169},{6,474,2.45882},{6,480,2.79414},{6,486,2.28356},{6,492,2.92389},{6,498,2.56894},{6,504,3.34326},{6,510,3.09164},{6,516,2.74151},{6,522,2.9421},{6,528,2.71061},{6,534,2.60536},{6,540,2.64271},{6,546,2.64826},{6,552,2.71457},{6,558,2.46986},{6,564,2.48719},{6,570,2.28519},{6,576,2.20632},{6,582,2.30377},{6,588,2.68933},{6,594,2.37031},{6,600,2.25064},{6,606,3.12559},{6,612,2.91535},{6,618,2.82004},{6,624,3.13842},{6,630,2.9301},{6,636,2.72388},{6,642,2.83937},{6,648,3.1243},{6,654,2.64058},{6,660,2.81251},{6,666,2.98226},{8,0,0.647083},{8,6,0.448819},{8,12,0.441372},{8,18,0.633842},{8,24,0.537127},{8,30,0.776977},{8,36,0.682159},{8,42,1.00228},{8,48,1.03336},{8,54,0.86845},{8,60,0.798648},{8,66,1.21947},{8,72,1.06527},{8,78,1.13433},{8,84,1.19649},{8,90,1.33958},{8,96,1.16118},{8,102,1.35805},{8,108,1.38032},{8,114,1.31299},{8,120,1.04441},{8,126,1.41363},{8,132,1.21928},{8,138,1.34441},{8,144,1.38456},{8,150,1.18587},{8,156,1.2505},{8,162,1.31973},{8,168,1.27488},{8,174,1.24806},{8,180,1.30616},{8,186,1.31515},{8,192,1.26549},{8,198,1.35743},{8,204,1.45514},{8,210,1.41813},{8,216,1.23963},{8,222,1.37289},{8,228,1.35309},{8,234,1.57017},{8,240,1.3751},{8,246,1.67303},{8,252,1.64209},{8,258,1.58623},{8,264,1.81903},{8,270,1.43628},{8,276,1.49471},{8,282,1.64683},{8,288,1.65975},{8,294,1.51223},{8,300,1.61269},{8,306,1.59499},{8,312,1.81595},{8,318,1.38265},{8,324,1.51008},{8,330,1.51451},{8,336,1.55474},{8,342,1.6666},{8,348,1.61281},{8,354,1.34824},{8,360,1.30458},{8,366,1.65495},{8,372,1.65538},{8,378,1.61235},{8,384,1.56087},{8,390,1.49175},{8,396,1.65757},{8,402,1.56799},{8,408,1.63076},{8,414,1.72441},{8,420,1.74952},{8,426,1.29295},{8,432,1.54347},{8,438,1.57697},{8,444,1.56348},{8,450,1.54073},{8,456,1.58157},{8,462,1.39144},{8,468,1.7468},{8,474,1.69158},{8,480,1.73185},{8,486,1.41786},{8,492,1.77357},{8,498,1.66453},{8,504,1.99042},{8,510,1.95037},{8,516,1.71123},{8,522,1.65104},{8,528,1.83421},{8,534,1.64569},{8,540,1.48086},{8,546,1.63882},{8,552,1.74687},{8,558,1.53387},{8,564,1.57706},{8,570,1.50582},{8,576,1.31741},{8,582,1.47397},{8,588,1.59465},{8,594,1.47101},{8,600,1.60332},{8,606,1.81642},{8,612,1.81595},{8,618,1.66391},{8,624,1.87916},{8,630,1.88918},{8,636,1.71072},{8,642,1.79382},{8,648,1.89856},{8,654,1.71103},{8,660,1.60331},{8,666,1.87195},{10,0,0.35654},{10,6,0.336505},{10,12,0.229219},{10,18,0.37277},{10,24,0.391889},{10,30,0.390797},{10,36,0.456228},{10,42,0.626833},{10,48,0.540682},{10,54,0.381638},{10,60,0.347661},{10,66,0.623156},{10,72,0.509822},{10,78,0.513299},{10,84,0.503252},{10,90,0.616285},{10,96,0.559117},{10,102,0.582412},{10,108,0.580618},{10,114,0.516496},{10,120,0.544474},{10,126,0.653425},{10,132,0.489784},{10,138,0.539393},{10,144,0.58631},{10,150,0.604227},{10,156,0.654631},{10,162,0.557624},{10,168,0.563593},{10,174,0.544871},{10,180,0.51509},{10,186,0.694853},{10,192,0.549551},{10,198,0.591988},{10,204,0.662093},{10,210,0.53143},{10,216,0.493782},{10,222,0.687478},{10,228,0.602734},{10,234,0.728726},{10,240,0.539286},{10,246,0.767974},{10,252,0.729111},{10,258,0.725236},{10,264,0.723239},{10,270,0.603335},{10,276,0.661284},{10,282,0.72075},{10,288,0.741072},{10,294,0.571773},{10,300,0.624261},{10,306,0.716177},{10,312,0.733817},{10,318,0.678443},{10,324,0.643768},{10,330,0.649156},{10,336,0.708197},{10,342,0.82293},{10,348,0.708308},{10,354,0.605526},{10,360,0.604033},{10,366,0.828415},{10,372,0.787685},{10,378,0.731335},{10,384,0.787792},{10,390,0.755715},{10,396,0.727333},{10,402,0.824938},{10,408,0.716367},{10,414,0.787467},{10,420,0.774437},{10,426,0.636991},{10,432,0.759897},{10,438,0.690772},{10,444,0.679318},{10,450,0.572061},{10,456,0.712393},{10,462,0.671355},{10,468,0.770058},{10,474,0.638785},{10,480,0.743968},{10,486,0.701233},{10,492,0.84526},{10,498,0.695755},{10,504,0.875228},{10,510,0.849722},{10,516,0.803206},{10,522,0.749636},{10,528,0.780403},{10,534,0.694566},{10,540,0.580113},{10,546,0.650045},{10,552,0.631019},{10,558,0.705208},{10,564,0.607808},{10,570,0.596751},{10,576,0.542486},{10,582,0.594083},{10,588,0.739288},{10,594,0.62833},{10,600,0.569983},{10,606,0.677631},{10,612,0.749733},{10,618,0.686289},{10,624,0.852324},{10,630,0.918744},{10,636,0.769948},{10,642,0.679418},{10,648,0.755305},{10,654,0.751323},{10,660,0.704026},{10,666,0.837176}}

$\endgroup$
  • $\begingroup$ It would be helpful if you shared at least a subset of your data. $\endgroup$ – JimB Apr 5 '17 at 15:37
  • $\begingroup$ @JimBaldwin Ok. I've just edited the question and put some points there. $\endgroup$ – J.Edwards Apr 5 '17 at 15:47
  • 1
    $\begingroup$ you should first verify TempDist works for reasonable values of the parameters before trying to model fit. $\endgroup$ – george2079 Apr 5 '17 at 15:58
  • $\begingroup$ @george2079 TempDistworks. The problem, I think, is the time it takes to evaluate it. $\endgroup$ – J.Edwards Apr 5 '17 at 17:10
5
$\begingroup$

Ok, so there are two parts to this problem. This first really boils down to how do I get NonlinearModelFit to stop throwing a tantrum? That can be handled pretty straightforwardly by requiring inputs to our model function TempDist to be numeric. The new function definition looks like:

TempDist[r_?NumericQ, z_?NumericQ, t_?NumericQ, a_?NumericQ, 
   Diff_?NumericQ, \[Alpha]_?NumericQ] := 
  NIntegrate[(ro^2) Exp[-(ro^2 + r^2)/(4 Diff (t - ti))] BesselI[0, 
     ro r/(2 Diff (t - ti))] Exp[-\[Alpha]*z + (\[Alpha]^2)*
       Diff (t - ti)]*
    Erfc[(2 Diff \[Alpha] (t - ti) - z)/(Sqrt[
         4 Diff (t - ti)])]/(Diff (t - ti)), {ro, 0, a}, {ti, 0, t}, 
   PrecisionGoal -> 4, MaxRecursion -> 10];

The second question requires a bit more consideration. Namely, does my model represent my data? If we take a look at the data you will notice that these data have some offset from the $\left(r,t\right)$ plane.

enter image description here

This is a problem because your function, as it is defined, can't have a nonzero value at $\left(r,t\right) = \left(0,0\right)$. Therefore, if you let the original model be defined as $f\left(r,t\right)$, construct a new model function $g\left(r,t\right)=f\left(r,t\right)+\Delta$, it is possible to fit an appropriate fit. Now, after some computation time, the call to NonlinearModelFit will yield a result:

nlm = NonlinearModelFit[data, 
  TempDist[r, 2, t, a, Diff, 0.023] + \[CapitalDelta], {{Diff, 
    0.5}, {a, .75}, {\[CapitalDelta], 34}}, {r, t}]
(*FittedModel[31.8606 +TempDist[r,2,t,2.90992,10.5851,0.023]]*)

Checking against the original data:

testdata = 
  Flatten[Table[{r, t, nlm[r, t]}, {r, Subdivide[0, 10, 10]}, {t, 
     Subdivide[0, 600, 10]}], 1];
Show[ListPlot3D[testdata], ListPointPlot3D[data]]

enter image description here

$\endgroup$
  • $\begingroup$ ! Thank you. Indeed, I've just realized that I forgot to subtract the baseline value (which you called $\Delta$) from the original data. I'm doing this right now and trying your suggestions. $\endgroup$ – J.Edwards Apr 5 '17 at 17:24
3
$\begingroup$

This is an extended comment. @Marchi has shown you how to fix your function so that NonlinearModelFit will provide a solution. What is recommended to do next is to examine the quality of the fit.

From the 3D plot above one sees that the assumption of constant variance assumed by NonlinearModelFit is just not true. Either the variability increases with increases predicted value or increases with increases of r. Another way to look at this is to construct estimates of the distribution of residuals for each value of r:

resids = nlm["FitResiduals"];
d0 = SmoothKernelDistribution[Pick[resids, data[[All, 1]], 0]];
d2 = SmoothKernelDistribution[Pick[resids, data[[All, 1]], 2]];
d4 = SmoothKernelDistribution[Pick[resids, data[[All, 1]], 4]];
d6 = SmoothKernelDistribution[Pick[resids, data[[All, 1]], 6]];
d8 = SmoothKernelDistribution[Pick[resids, data[[All, 1]], 8]];
d10 = SmoothKernelDistribution[Pick[resids, data[[All, 1]], 10]];

And now a plot of the estimated density functions for each value of r:

Plot[{PDF[d0, x], PDF[d2, x], PDF[d4, x], PDF[d6, x], PDF[d8, x], PDF[d10, x]},
 {x, -3, 3}, PlotRange -> {Automatic, {0, 2.6}},
 PlotStyle -> {{Thickness[0.007], Red}, {Thickness[0.007], Black}, {Thickness[0.007], Blue},
   {Thickness[0.007], Green}, {Thickness[0.007], Cyan}, {Thickness[0.007], Orange}},
 ImageSize -> Large, 
 PlotLegends -> {"r = 0", "r = 2", "r = 4", "r = 6", "r = 8", "r = 10"}]

Residual distributions for each value of r

One can see that estimate of constant variability is not justified and there is a considerable lack of fit: none of the residual distributions seem to be centered around zero.

$\endgroup$

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.