0

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}}*)
CAF
  • 510
  • 2
  • 13
  • 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. – SPPearce Nov 15 '17 at 12:56
  • @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? – CAF Nov 15 '17 at 16:16
  • Yes, you can't call NIntegrate without giving numerical values for all the arguments. – SPPearce Nov 16 '17 at 10:53

1 Answers1

1

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.

SPPearce
  • 5,653
  • 2
  • 18
  • 40
  • 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? – CAF Nov 16 '17 at 13:13
  • 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. – SPPearce Nov 16 '17 at 13:38
  • 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? – CAF Nov 16 '17 at 16:27
  • It is a symbolic derivative, but here you have a numerical function. That is really not Mathematica-style code though. Read https://mathematica.stackexchange.com/questions/134609/why-should-i-avoid-the-for-loop-in-mathematica and https://mathematica.stackexchange.com/questions/7924/alternatives-to-procedural-loops-and-iterating-over-lists-in-mathematica. – SPPearce Nov 17 '17 at 06:28
  • I see, thanks +1, I also replaced alphas[k_] = 4Pi/9/Log[k/lambda] with alphas[k_] = 4Pi/9/Log[k/lambda] - 16Pi16/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. – CAF Nov 18 '17 at 15:25
  • That looks fine for me. Test it by going alphas[1], if that isn't a number that is your issue. Did you manage to sort out your derivatives? – SPPearce Nov 19 '17 at 20:46
  • I see that the code runs with this new alphas but there are errors about integrand evaluated to non-numeric values and something about 'General: 7.0432....' is not a valid variable. I don't really understand this. Would you be able to check if you get the same errors on your machine? – CAF Nov 19 '17 at 21:50
  • That error usually means that you are trying to integrate or differentiate with respect to a variable that you have a value assigned to. Have you tried closing the kernel and trying again? It works for me (although it is worse for the specific a,b,NN that was doing well earlier) – SPPearce Nov 20 '17 at 11:16
  • Thanks for the check. I've tried it multiple times and I tried again just now and I got the same errors. I've edited my OP ^^ with the code I'm working with and the errors that I'm getting (commented at the end). I've also commented the (worsened) chi square it returns and a,b,NN. So it runs, but just gives me these errors in addition - should I be worried about this? – CAF Nov 20 '17 at 12:03
  • Yes, you should be worried about that, you can't just ignore those type of errors. It says that the integral was non-numeric. – SPPearce Nov 20 '17 at 13:22
  • The code in the first post works for me. Have you closed the kernel properly, and make sure you aren't defining a variable somewhere earlier. k is the likely culprit, as you posted some code earlier using k as an iterator in a For loop. So type k and see if it has a value (and if it does you can clear it with Clear[k] or k=. ) – SPPearce Nov 20 '17 at 13:25
  • Yup I tried that , I’ll maybe try running the code on a desktop at the university rather than my personal computer and see what happens. Just out of interest, was the 160 chi square that my code returned in line with what you got ? I should also say when I resort to using the previous formula for alphas everything works fine. – CAF Nov 20 '17 at 13:38
  • Yes, I get the same 160ish value. I'd recommend you track down where that error is coming from though. – SPPearce Nov 20 '17 at 13:46
  • Ok Thanks for your help! The last thing I would ask is how to make sure my chi square function, which I'm differentiating with respect to the parameters NN,a,b for returning the covariance matrix, is symbolic and not numeric? – CAF Nov 20 '17 at 21:09
  • Did you catch my last comment? :) – CAF Dec 20 '17 at 12:01
  • You should probably ask that as a separate question. – SPPearce Dec 21 '17 at 12:04
  • thanks, I have done so, it is here if you are able to help on this one, https://mathematica.stackexchange.com/questions/162655/how-to-perform-a-symbolic-derivative-with-respect-to-numerically-defined-variabl – CAF Dec 27 '17 at 17:10