0
$\begingroup$

I'm trying to fit a model curve to some data by performing a chi square minimisation wrt three parameters $a,b$ and $NN$. The trouble I am having is that the variables with which I want to minimise the chi square with respect to appear in an integral. I attach the code I am working with.

Needs["ErrorBarPlots`"]

MasData5 = {{{44.8, 47.5}, ErrorBar[4.0]}, {{54.8, 50.1}, 
ErrorBar[4.2]}, {{64.8, 61.7}, ErrorBar[5.1]}, {{74.8, 64.8}, 
ErrorBar[5.5]}, {{84.9, 75}, ErrorBar[6.2]}, {{94.9, 81.2}, 
ErrorBar[6.7]}, {{104.9, 85.3}, ErrorBar[7.1]}, {{119.5, 94.5}, 
ErrorBar[7.5]}, {{144.1, 101.5}, ErrorBar[8.3]}, {{144.9, 101.9}, 
ErrorBar[10.9]}, {{162.5, 117.8}, 
ErrorBar[12.8]}, {{177.3, 130.2}, 
ErrorBar[13.4]}, {{194.8, 147.7}, 
ErrorBar[17.1]}, {{219.6, 137.4}, 
ErrorBar[20.1]}, {{244.8, 176.6}, 
ErrorBar[20.3]}, {{267.2, 178.7}, 
ErrorBar[21.1]}, {{292.3, 200.4}, ErrorBar[29.1]}, {{60, 55.8}, 
ErrorBar[4.838]}, {{80, 66.6}, ErrorBar[7.280]}, {{100, 73.4}, 
ErrorBar[6.426]}, {{120, 86.7}, ErrorBar[7.245]}, {{140, 104}, 
ErrorBar[12.083]}, {{160, 110}, ErrorBar[16.279]}, {{42.5, 43.8}, 
ErrorBar[3.482]}, {{55, 57.2}, ErrorBar[3.980]}, {{65, 62.5}, 
ErrorBar[4.614]}, {{75, 68.9}, ErrorBar[5.197]}, {{85, 72.1}, 
ErrorBar[5.523]}, {{100, 81.9}, ErrorBar[5.368]}, {{117.5, 95.7}, 
ErrorBar[6.277]}, {{132.5, 103.9}, ErrorBar[6.912]}, {{155, 115}, 
ErrorBar[7.920]}, {{185, 129.1}, ErrorBar[9.192]}, {{215, 141.7}, 
ErrorBar[10.666]}, {{245, 140.3}, ErrorBar[14.526]}, {{275, 189}, 
ErrorBar[24.274]}, {{49, 39.2}, ErrorBar[10]}, {{86, 75.7}, 
ErrorBar[14.414]}, {{167, 118}, ErrorBar[22.828]}, {{43.2, 50.7}, 
ErrorBar[1.5]}, {{50, 59.5}, ErrorBar[1.4]}, {{57.3, 61.8}, 
ErrorBar[1.9]}, {{65.3, 67.6}, ErrorBar[1.7]}, {{73.9, 72.4}, 
ErrorBar[1.9]}, {{83.2, 79.9}, ErrorBar[2.3]}, {{93.3, 84.4}, 
ErrorBar[2.1]}, {{104.3, 86.7}, ErrorBar[2.7]}, {{47.9, 55.4}, 
ErrorBar[2.1]}, {{68.4, 66.4}, ErrorBar[2.9]}}; (*this is the experimental data *)

gamma = 5.55*^-6;
MJpsi = 3.1;
alphaem = 1/137;
alphas[k_] = 4*Pi/9/Log[k/lambda]; 
lambda = 0.09;
Ca = 3; (*list of constants and k dependent function alphas*)

xg = NN*((4*qbar)/((4*qbar - MJpsi^2) + w^2))^(-a)*(qbar)^b*Exp[Sqrt[16*Ca/9*Log[(4*qbar - MJpsi^2 + w^2)/(4*qbar)]*
 Log[Log[qbar/lambda]/Log[qo/lambda]]]] (* def xg, appearing in below *)

F5[w_] = 3.89379*^5*1/4.5/16*4*Pi^3*MJpsi^3*
gamma/12/
 alphaem*(Nintegrate[
    alphas[k_]/qbar/(qbar + k_)*
     D[(2^(2*(D[Log[xg], 
                Log[1/(((4*qbar))/((4*qbar - MJpsi^2) + 
                w^2))]]) + 3)/Sqrt[Pi]*
         Gamma[D[Log[xg], 
             Log[1/(((4*qbar))/((4*qbar - MJpsi^2) + w^2))]] + 
            5/2]/ Gamma[
           D[Log[xg], 
             Log[1/(((4*qbar))/((4*qbar - MJpsi^2) + w^2))]] + 
            4])*NN* (((4*qbar))/((4*qbar - MJpsi^2) + 
            w^2))^(-a)*(k_)^b*
       Exp[Sqrt[
         16*Ca/9*Log[(4*qbar - MJpsi^2 + w^2)/(4*qbar)]*
          Log[Log[k_/lambda]/Log[qo/lambda]]]], k_], {k_, 
     qo, (w^2 - MJpsi^2)/4}] + 
   0.118/qbar/qo*Log[(qbar + qo)/qbar]*
    NN* ((4*qbar)/((4*qbar - MJpsi^2) + w^2))^(-a)*(qo)^
     b*(2^(2*(D[Log[xg], 
             Log[1/(((4*qbar))/((4*qbar - MJpsi^2) + w^2))]]) + 
          3)/Sqrt[Pi]*
      Gamma[D[Log[xg], 
          Log[1/(((4*qbar))/((4*qbar - MJpsi^2) + w^2))]] + 5/2]/ 
       Gamma[D[Log[xg], 
          Log[1/(((4*qbar))/((4*qbar - MJpsi^2) + w^2))]] + 
         4]))^2; /. { qbar -> 2.4025, qo -> 2} (*my model function - involves a numerical integration over variable k_ *)

  chisq5 = Sum[((MasData5[[k]][[1]][[2]] - F5[MasData5[[k]][[1]][[1]]])/
  MasData5[[k]][[2]][[1]])^2, {k, 1, Length[MasData5]}]; (*def of my chi square function *)

  rr = Minimize[chisq5, {a, b, NN}] (*minimise chi square wrt a,b,NN *)

This gives me an error that the output is not a number at {a,b,NN}. Can mathematica handle minimisation of a function wrt variables that are involved in an integration? I would have thought yes as the numerical integration involved is simply a nested operation but perhaps I have to adapt my syntax. Or perhaps the error is simply something else I overlooked? Thanks for any pointers.

Edit:

   Needs["ErrorBarPlots`"]

   MasData5 = {{{44.8, 47.5}, ErrorBar[4.0]}, {{54.8, 50.1}, 
ErrorBar[4.2]}, {{64.8, 61.7}, ErrorBar[5.1]}, {{74.8, 64.8}, 
ErrorBar[5.5]}, {{84.9, 75}, ErrorBar[6.2]}, {{94.9, 81.2}, 
ErrorBar[6.7]}, {{104.9, 85.3}, ErrorBar[7.1]}, {{119.5, 94.5}, 
ErrorBar[7.5]}, {{144.1, 101.5}, ErrorBar[8.3]}, {{144.9, 101.9}, 
ErrorBar[10.9]}, {{162.5, 117.8}, 
ErrorBar[12.8]}, {{177.3, 130.2}, 
ErrorBar[13.4]}, {{194.8, 147.7}, 
ErrorBar[17.1]}, {{219.6, 137.4}, 
ErrorBar[20.1]}, {{244.8, 176.6}, 
ErrorBar[20.3]}, {{267.2, 178.7}, 
ErrorBar[21.1]}, {{292.3, 200.4}, ErrorBar[29.1]}, {{60, 55.8}, 
ErrorBar[4.838]}, {{80, 66.6}, ErrorBar[7.280]}, {{100, 73.4}, 
ErrorBar[6.426]}, {{120, 86.7}, ErrorBar[7.245]}, {{140, 104}, 
ErrorBar[12.083]}, {{160, 110}, ErrorBar[16.279]}, {{42.5, 43.8}, 
ErrorBar[3.482]}, {{55, 57.2}, ErrorBar[3.980]}, {{65, 62.5}, 
ErrorBar[4.614]}, {{75, 68.9}, ErrorBar[5.197]}, {{85, 72.1}, 
ErrorBar[5.523]}, {{100, 81.9}, ErrorBar[5.368]}, {{117.5, 95.7}, 
ErrorBar[6.277]}, {{132.5, 103.9}, ErrorBar[6.912]}, {{155, 115}, 
ErrorBar[7.920]}, {{185, 129.1}, ErrorBar[9.192]}, {{215, 141.7}, 
ErrorBar[10.666]}, {{245, 140.3}, ErrorBar[14.526]}, {{275, 189}, 
ErrorBar[24.274]}, {{49, 39.2}, ErrorBar[10]}, {{86, 75.7}, 
ErrorBar[14.414]}, {{167, 118}, ErrorBar[22.828]}, {{43.2, 50.7}, 
ErrorBar[1.5]}, {{50, 59.5}, ErrorBar[1.4]}, {{57.3, 61.8}, 
ErrorBar[1.9]}, {{65.3, 67.6}, ErrorBar[1.7]}, {{73.9, 72.4}, 
ErrorBar[1.9]}, {{83.2, 79.9}, ErrorBar[2.3]}, {{93.3, 84.4}, 
ErrorBar[2.1]}, {{104.3, 86.7}, ErrorBar[2.7]}, {{47.9, 55.4}, 
ErrorBar[2.1]}, {{68.4, 66.4}, ErrorBar[2.9]}};
(*h1 2006 Q^2=0 data,zeus 2002,zeus 2004 and h1 2013 data for Q^2=0*)

gamma = 5.55*^-6;
MJpsi = 3.1;
alphaem = 1/137;
alphas[k_] = 4*Pi/9/Log[k/lambda] - 
16*16*Pi/9/9/9*Log[Log[k/lambda]]/(Log[k/lambda])^2;
lambda = 0.09;
Ca = 3;

xg = NN*((4*qbar)/((4*qbar - MJpsi^2) + w^2))^(-a)*(qbar)^b*Exp[Sqrt[
 16*Ca/9*Log[(4*qbar - MJpsi^2 + w^2)/(4*qbar)]*
  Log[Log[qbar/lambda]/Log[qo/lambda]]]];

F5[w_?NumericQ, a_?NumericQ, b_?NumericQ, NN_?NumericQ] := 
Module[{qbarr = 2.4025, qoo = 2, qbar, qo}, (3.89379*^5*1/(4.9 + 4*0.06*Log[w/90])/16*4*Pi^3*MJpsi^3*
  gamma/12/
   alphaem*(NIntegrate[
      alphas[k]/qbar/(qbar + k)*
         D[(2^(2*(D[Log[xg], 
                Log[1/(((4*qbar))/((4*qbar - MJpsi^2) + w^2))]]
                + 3))/Sqrt[Pi]*
             Gamma[D[Log[xg], 
                Log[1/(((4*qbar))/((4*qbar - MJpsi^2) + w^2))]]
                + 5/2]/
              Gamma[D[Log[xg], 
                Log[1/(((4*qbar))/((4*qbar - MJpsi^2) + w^2))]]
                + 4])*
           NN*(((4*qbar))/((4*qbar - MJpsi^2) + w^2))^(-a)*(k)^b*
           Exp[Sqrt[
             16*Ca/9*Log[(4*qbar - MJpsi^2 + w^2)/(4*qbar)]*
              Log[Log[k/lambda]/Log[qo/lambda]]]], k] /. 
        qbar -> qbarr /. qo -> qoo, {k, 
       qoo, (w^2 - MJpsi^2)/4}] + 
     0.118/qbar/qo*Log[(qbar + qo)/qbar]*
      NN*((4*qbar)/((4*qbar - MJpsi^2) + w^2))^(-a)*(qo)^
       b*(2^(2*(D[Log[xg], 
               Log[1/(((4*qbar))/((4*qbar - MJpsi^2) + w^2))]]
              + 3))/Sqrt[Pi]*
        Gamma[D[Log[xg], 
            Log[1/(((4*qbar))/((4*qbar - MJpsi^2) + w^2))]]
           + 5/2]/
         Gamma[D[Log[xg], 
            Log[1/(((4*qbar))/((4*qbar - MJpsi^2) + w^2))]]
           + 4]))^2) /. qbar -> qbarr /. qo -> qoo]

   chisq5[a_?NumericQ, b_?NumericQ, NN_?NumericQ] := Sum[((MasData5[[k, 1, 2]] - F5[MasData5[[k, 1, 1]], a, b, NN])/
  MasData5[[k, 2, 1]])^2, {k, 1, Length[MasData5]}];

   NMinimize[chisq5[a, b, NN], {a, b, NN}]

  (* General::ivar: 7.043210807822302` is not a valid variable.
     General::ivar: 7.043210807822302` is not a valid variable.
     General::ivar: 7.043210807822302` is not a valid variable.
     General::stop: Further output of General::ivar will be suppressed during this calculation.
     NIntegrate::inumr: The integrand (0.416233 (<<1>>/(k^<<18>> Gamma[4+<<1>><<1>>])-(<<18>> <<2>> <<1>>)/(k^<<18>> <<3>>)) ((4 \[Pi])/(9 Log[11.1111 k])-(256 \[Pi] Log[Log[11.1111 k]])/(729 Log[<<18>> k]^2)))/(2.4025 +k) has evaluated to non-numerical values for all sampling points in the region with boundaries {{2,2748.6}}.
     NIntegrate::inumr: The integrand (0.416233 (<<1>>/(k^<<18>> Gamma[4+<<1>><<1>>])-(<<18>> <<2>> <<1>>)/(k^<<18>> <<3>>)) ((4 \[Pi])/(9 Log[11.1111 k])-(256 \[Pi] Log[Log[11.1111 k]])/(729 Log[<<18>> k]^2)))/(2.4025 +k) has evaluated to non-numerical values for all sampling points in the region with boundaries {{2,2748.6}}.
    NIntegrate::inumr: The integrand (0.416233 (<<1>>/(k^<<18>> Gamma[4+<<1>><<1>>])-(<<18>> <<2>> <<1>>)/(k^<<18>> <<3>>)) ((4 \[Pi])/(9 Log[11.1111 k])-(256 \[Pi] Log[Log[11.1111 k]])/(729 Log[<<18>> k]^2)))/(2.4025 +k) has evaluated to non-numerical values for all sampling points in the region with boundaries {{2,2748.6}}.
    General::stop: Further output of NIntegrate::inumr will be suppressed during this calculation.
    NMinimize::nnum: The function value 413148. +0.0198373 (85.3 -1154.08 (-0.635717+<<1>>)^2)^2+0.0253802 (95.7 -<<18>> <<1>>)^2+<<21>> <<1>>+0.0190512 <<1>>^2+0.0209311 (103.9 -1141.12 (<<1>>)^2)^2+0.00684937 (104-1138.11 (-0.784293+NIntegrate[<<1>>,{<<3>>}])^2)^2 is not a number at {a,b,NN} = {0.363833,-3.72138,-5.43375}. *)

     (* {160.103, {a -> -0.0482274, b -> -1.9486, NN -> -2.66949}}*)
Improve this question
$\endgroup$
3
  • $\begingroup$ By default Minimize will try and and do symbolic manipulation before plugging numbers in, you need to make your function only be defined for numeric values. But that isn't the only issue here. You have Nintegrate instead of NIntegrate, you try and apply the substitutions qbar and qo at the end of the expression when it needs to be within the limits of NIntegrate, you are trying to use k_ (a pattern) as the variable in the integration. I'd also check whether the differentiation is doing what you expect, as you are differentiating with respect to an expression. $\endgroup$
    – SPPearce
    Commented Nov 15, 2017 at 12:56
  • $\begingroup$ @Krazug, Thanks! I've made your suggestions but can't seem to find a way to do the NIntegrate without first specifying what NN,a and b are. But of course they are to determined through minimisation so is there an easy way to proceed? $\endgroup$
    – CAF
    Commented Nov 15, 2017 at 16:16
  • $\begingroup$ Yes, you can't call NIntegrate without giving numerical values for all the arguments. $\endgroup$
    – SPPearce
    Commented Nov 16, 2017 at 10:53

1 Answer 1

Reset to default
1
$\begingroup$

Try something like this:

F5[w_?NumericQ, a_?NumericQ, b_?NumericQ, NN_?NumericQ] :=
 Module[{qbarr = 2.4025, qoo = 2, qbar, qo},
  (3.89379*^5*1/4.5/16*4*Pi^3*MJpsi^3*gamma/12/alphaem*
   (NIntegrate[alphas[k]/qbar/(qbar + k)*
             D[(2^(2*(D[Log[xg], Log[1/(((4*qbar))/((4*qbar - MJpsi^2) + w^2))]]) + 3)/Sqrt[Pi]*
                 Gamma[D[Log[xg], Log[1/(((4*qbar))/((4*qbar - MJpsi^2) + w^2))]] + 5/2]/Gamma[D[Log[xg], Log[1/(((4*qbar))/((4*qbar - MJpsi^2) + w^2))]] + 4])*NN*(((4*qbar))/((4*qbar - MJpsi^2) + w^2))^(-a)*(k)^b*Exp[Sqrt[16*Ca/9*Log[(4*qbar - MJpsi^2 + w^2)/(4*qbar)]*Log[Log[k/lambda]/Log[qo/lambda]]]], k] /. qbar -> qbarr /. qo -> qoo, {k, qoo, (w^2 - MJpsi^2)/4}] + 
         0.118/qbar/qo*Log[(qbar + qo)/qbar]*NN*((4*qbar)/((4*qbar - MJpsi^2) + w^2))^(-a)*(qo)^b*(2^(2*(D[Log[xg],Log[1/(((4*qbar))/((4*qbar - MJpsi^2) + w^2))]]) + 3)/Sqrt[Pi]* Gamma[D[Log[xg], Log[1/(((4*qbar))/((4*qbar - MJpsi^2) + w^2))]] +5/2]/Gamma[D[Log[xg], Log[1/(((4*qbar))/((4*qbar - MJpsi^2) + w^2))]] + 4]))^2) /. qbar -> qbarr /. qo -> qoo] 

chisq5[a_?NumericQ, b_?NumericQ, NN_?NumericQ] := 
  Sum[((MasData5[[k, 1, 2]] - F5[MasData5[[k, 1, 1]], a, b, NN])/
      MasData5[[k, 2, 1]])^2, {k, 1, Length[MasData5]}];

I've used a Module, in order to localise the valiues of qbar and qo, because you need to evaluate the derivatives involving them before plugging numbers in. I would definitely check that those derivatives are actually doing what you expect though, because differentiation with respect to an expression is something that I've had problems with in the past.

Having done that, you get a numerical value when you evaluate e.g. chisq5[1,1,1]. You'll need to use NMinimize though.

NMinimize[chisq5[a, b, NN], {a, b, NN}]

You may want to try something like the following definition of chisq5, to cache the best value that you've got so far so you can interrupt the NMinimize. You may also want to check out FindMinimum instead of NMinimize

min = {∞, {1, 1, 1}}
chisq5[a_?NumericQ, b_?NumericQ, 
   NN_?NumericQ] := (val = 
    Sum[((MasData5[[k, 1, 2]] - F5[MasData5[[k, 1, 1]], a, b, NN])/
        MasData5[[k, 2, 1]])^2, {k, 1, Length[MasData5]}]; 
   If[val < min[[1]], min = {val, {a, b, NN}}]; val);

The best value that running it for a few minutes for me is a=-0.141347,b=-0.362054,NN=1.16421, which has a value of 55.87, which seems a pretty good start.

Improve this answer
$\endgroup$
17
  • $\begingroup$ Oh thanks very much! I tried the code you provided and obtained values of 54.7837, a-> -0.128017, b->-0.456918 and NN->-1.24031. However, NMinimize failed to converge after 100 iterations. I just wonder is there a way to extend the number of iterations NMinimize can do so that maybe I can obtain the same numbers that you got? $\endgroup$
    – CAF
    Commented Nov 16, 2017 at 13:13
  • $\begingroup$ I switched to FindMinimum actually (once I had a reasonable value from the initial NMinimize I wanted to zoom in on it). Yes, you can change the maximum number of iterations (and various other things, such as the method used), look in the help files for how to do that. $\endgroup$
    – SPPearce
    Commented Nov 16, 2017 at 13:38
  • $\begingroup$ Thanks! Found it, the last thing I had a quick query on was in my determination of the covariance matrix. The code is Ucovinv = IdentityMatrix[3]; vecpar = {a, b, NN}; For[i = 1, i <= 3, i++, For[k = 1, k <= 3, k++, Ucovinv[[i]][[ k]] = (1/ 2 D[chisq5[a, b, NN], vecpar[[i]], vecpar[[k]]] /. %19[[2]])]]; Ucov = Inverse[Ucovinv]; Ucov // MatrixForm, but this returns a matrix with peculiar output such as chisq5 raised to a power (0,1,1). Do you know what that means? $\endgroup$
    – CAF
    Commented Nov 16, 2017 at 16:27
  • $\begingroup$ It is a symbolic derivative, but here you have a numerical function. That is really not Mathematica-style code though. Read mathematica.stackexchange.com/questions/134609/… and mathematica.stackexchange.com/questions/7924/…. $\endgroup$
    – SPPearce
    Commented Nov 17, 2017 at 6:28
  • $\begingroup$ I see, thanks +1, I also replaced alphas[k_] = 4*Pi/9/Log[k/lambda] with alphas[k_] = 4*Pi/9/Log[k/lambda] - 16*Pi*16/9/9/9*Log[Log[k/lambda]]/(Log[k/lambda])^2 (an improved functional form for alphas) - however this gives errors about non-numeric input and I don't see why they come about by simply redefining my functional form for alphas. $\endgroup$
    – CAF
    Commented Nov 18, 2017 at 15:25

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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