# Numerical Integration as Model for Nonlinear Fit

I'm trying to do a fit for parameters of a function within an integral but I'm getting errors when I try and run it. Essentially I want the following fit to work out:

model[a_, b_, x_] := Sum[NIntegrate[totalνLayer[p, m] (a (1 - b^-p)),
{p, 0.2, 10}], {m, x - .45, x + .45, .02}]

nlm = NonlinearModelFit[data, model[a, b, x], {a,b}, x, Weights -> 1/errors^2]

Where totalνLayer is a function I defined earlier. I just want the fit values for a and b and their covariance matrix, but I'm clearly doing something wrong. Any help is much appreciated!

For those who care, all of the necessary code is below:

amp[x_] := 47.5111 - 1.79354 Abs[x]
γ[x_] := 2.5587996836804026 Exp[-0.20766051151975093 Abs[x]] + 0.5696086481954964

crossSection[x_] := -600.089 + 1201.5577451256304 x - 1751.9101332336447 x^2 + 1899.988095515613 x^3 - 1556.709382526493 x^4 + 974.3641634962205 x^5 - 468.5129814014657 x^6 + 172.90071370757002 x^7 - 48.48328424613457 x^8 + 10.067733931487915 x^9 - 1.4558045159248463 x^10 + 0.1215353564558521 x^11 -0.00009554527467999401 x^12 - 0.001211849211458331 x^13 + 0.00009244381215772915 x^14 + 6.85533417357146*^-6 x^15 - 1.3215490607411126*^-6 x^16 - 1.213701585221995*^-8 x^17 + 1.4963340473233928*^-8 x^18 - 6.405376189932708*^-10 x^19 - 1.3047056106410797*^-10 x^20 + 2.0016864863803638*^-11 x^21 - 1.2643359066392302*^-12 x^22 + 4.036927989211094*^-14 x^23 - 5.368234695690585*^-16 x^24 + 214.8882519560799 x^-1 - 54.26505937368572 x^-2 + 9.480463574208194 x^-3 - 1.1135829819290948 x^-4 + 0.08354380635765371 x^-5 - 0.0036045838406917876 x^-6 + 0.0000678896567028218 x^-7

And the definition of totalνLayer:

crossSection[x] 10^-38*1*6300*6.022*10^23)

And the data:

data = {{0, 1185.2847069383793}, {4, 780.6026466248076}, {8, 493.5502399832769},
{12, 312.1838499130911}, {16, 193.16881196193143}, {20, 113.41366552793787},
{24, 64.20654682074144}};

errors = {132.12627600307064, 90.55918065608893, 60.884899393246535,
41.88064141772561, 29.076284534300132, 20.049728493077968,
13.950371635495452};

And last but not least the model and fit:

model[a_, b_, x_] := Sum[NIntegrate[
totalνLayer[p, m] (a (1 - b^-p))*10^-38*1.4*150*6.022*10^23,
{p, 0.2, 10}],
{m, x - .45, x + .45, .02}]

nlm = NonlinearModelFit[data, model[α, β, x], {α, β}, x, Weights -> 1/errors^2]

I apologize for all of the weird numbers and functions; this is neutrino data from my research.

-
Can't try it right now, but does it help when you change the model definition LHS to model[a_?NumericQ,b_?NumericQ,x_?NumericQ] ? –  Sjoerd C. de Vries Jul 10 '13 at 6:20
I just tried it and it still doesn't work. –  feynmansafineman Jul 10 '13 at 8:43

So, this isn't a very charming answer, and falls under the category of brute force and ignorance, (with a heavy dose of the latter).

The problem is that NonlinearModelFit appears to behave similarly to FindMinimum and FindRoot, where the function is initially evaluated with the parameters being symbolic. The Nintegrate within your model requires that the parameters are evaluated numerically. I'll let the more educated Mathematica users chime in here because the suggestion to replace model[a_,b_,x_] with model[a_?NumericQ,b_?NumericQ,x_] should have worked if this were the only problem.

Here's my attempt to help you out and find an answer. A little digging in your code makes me think that the values for the parameters a and b lie in the order of magnitude $10^{-1}$ and $10^0$, respectively. I'm going to choose a = 0.2 and b = 2 and I get:

Not too bad for a not-so-educated guess. We can calculate the sum of squares

Total[(data[[All, 2]] - trial1[[All, 2]])^2]
(* 4771.29 *)

and use this as a starting point to improve our fit. Evaulating the sum of squares at various points a la

out = Table[{a, b, Total[(data[[All, 2]] -
Evaluate[model[a, b, #] & /@ data[[All, 1]]])^2]},
{a, 0.18, 0.22, 0.01}, {b, 1.8, 2.2, 0.1}]
ListPlot3D[Partition[Flatten[out], 3]]

I then search for the minimum:

Position[out, Min[Partition[Flatten[latest], 3][[All, 3]]]]
(* {{2,5,3}} *)
out[[2,5]]
(* {0.19, 2.2, 644.488} *)

So the brute force and ignorance approach is working. A few more iterations and the best I've found is a = 0.185, b = 2.5 (sum squares = 95.7263), which gives us:

ListPlot[{data, {#, model[0.185, 2.5, #]} & /@ data[[All, 1]]},
PlotStyle -> PointSize[0.02]]

Parameter b needs to be optimized further since it is continually on the upper edge of the range that I search. I fully expect one of the 10k+ rep guys to come in here and say "all you need to do is Reap Sow HoldPattern" and save you all this work, but I hope this exercise at least shows you one type of work around to get to the answer you are looking for.

And yes, I do know that you wanted a weighted least squares, but I have no idea how to do an a priori graphical weighted nonlinear least squares analysis. Sorry about that.

-
Thanks for this! Do you have any ideas on how to use the errors and also somehow get a covariance matrix from this? –  feynmansafineman Jul 15 '13 at 7:17
@feynmansafineman That's well beyond my statistics knowledge, sorry. –  bobthechemist Jul 16 '13 at 12:05

To "use the errors" you need to do the following:

Needs["ErrorBarPlots"]

{r, t} = {8.31, 0.002};

data = {{1244, 0.90237}, {1289, 0.83072}, {1306, 0.76147}, {1360, 0.6094},
{1388, 0.40952}, {1428, 0.29419}, {1488, 0.15592}, {1494, 0.11718},
{1527, 0.09192}, {1537, 0.07757}, {1590, 0.04546},
{1594, 0.04563}, {1654, 0.03229}, {1737, 0.01976}, {1740, 0.01949}};

nlm = NonlinearModelFit[
data,
{Exp[-a*tt^n*Exp[-(ea/(r*tt))]*t], 10^12 < a < 10^14,
27*10^4 < ea < 32*10^4,
0 < n < 0.9},
{{a, 10^14}, {n, 0}, {ea, 30*10^4}},
tt, Method -> "NMinimize"]

ErrorListPlot[{{{0, 1185.2847069383793},
ErrorBar[132.12627600307064]}, {{4, 780.6026466248076},
ErrorBar[90.55918065608893]}, {{8, 493.5502399832769},
ErrorBar[60.884899393246535]}, {{12, 312.1838499130911},
ErrorBar[41.88064141772561]}, {{16, 193.16881196193143},
ErrorBar[29.076284534300132]}, {{20, 113.41366552793787},
ErrorBar[20.049728493077968]}, {{24, 64.20654682074144},
ErrorBar[13.950371635495452`]}},
Axes -> False]

Output:

nlm["CovarianceMatrix"] // MatrixForm

Output:

$\left( \begin{array}{ccc} \text{6.434561442695655$\grave{ }$*${}^{\wedge}$-11} & -\text{3.901357230169781$\grave{ }$*${}^{\wedge}$-6} & -0.325294 \\ -\text{3.901357230169781$\grave{ }$*${}^{\wedge}$-6} & 0.236838 & 19723. \\ -0.325294 & 19723. & 1.64449\times 10^9 \\ \end{array} \right)$

-