Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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.

share|improve this question
    
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
add comment

2 Answers

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:

first attempt

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]]

finding minimum

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]] 

Best fit so far

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.

share|improve this answer
    
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
add comment

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:

Error Plot

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

share|improve this answer
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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