Numerical Integration as Model for Nonlinear Fit

Question

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]
μ[mrad_] := 5.449517720787487 Exp[-0.0693211022591832 Abs[mrad]]
σ[mrad_] := 1.7679270770338407 Exp[-0.12510841305113626 Abs[mrad]]
γ[x_] := 2.5587996836804026 Exp[-0.20766051151975093 Abs[x]] + 0.5696086481954964
ω[mrad_] := 7.63444841081276 Exp[-0.08637416159312404 Abs[mrad]]
lin[mrad_] := 0.9648001654360766 - 0.10434481940699446 Abs[mrad] + 0.003062473453032527 Abs[mrad]^2

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:

totalνLayer[x_, mrad_] := .3*3.78*10^-5*810^2*5 (amp[mrad] Exp[-(x - μ[mrad])^2/(
                          2 σ[mrad]^2)] + lin[mrad] x + (amp[mrad] γ[mrad]^2)/(
                          4 (γ[mrad]^2/4 + (x - ω[mrad])^2)))/(
                          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.

Answer 1

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.

Answer 2

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)$