Chi square minimisation wrt variables in an integral?

[CAF123]

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.

 

[phyzguy]

Ignore what I wrote below. I see that you are trying to differentiate with respect to a defined function in your expression:
D[Log[xg], Log[1/(((4*qbar))/((4*qbar - MJpsi^2) + w^2))]]

I don't see how this is ever going to work.Mathematica should be able to do this. I got it to run by making the following changes:

(1) You were using the same variable k for the integration variable and the index of the MasData5 list. I suspect this isn't what you wanted.
(2) I explicitly defined q0 and qbar, rather than assigning them. This may not be necessary.
(3) I made xg, F5, and chisq5 explicit functions with delayed evaluation. (the :=) I think you need to do this in order for Minimize to work.
(4) You had alphas[k_] on the RHS of F5. This won't work. It should be alphas[k].

At any rate, it now runs, but returns all zeros for the minimum values. It may not be doing what you want, but at least it runs now.

In[1]:= Needs["ErrorBarPlots`"]
In[2]:= 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*)
In[3]:= 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*)
qbar=2.4025;
qo=2.0;

xg[w_, a_, b_, NN_]:=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*)
In[12]:=
F5[w_, a_, b_, NN_] :=
 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[w, a, b, NN]],
                  Log[1/(((4*qbar))/((4*qbar - MJpsi^2) + w^2))]]) +
               3)/Sqrt[Pi]*
           Gamma[D[Log[xg[w, a, b, NN]],
               Log[1/(((4*qbar))/((4*qbar - MJpsi^2) + w^2))]] + 5/2]/
            Gamma[D[Log[xg[w, a, b, NN]],
               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[w, a, b, NN]],
               Log[1/(((4*qbar))/((4*qbar - MJpsi^2) + w^2))]]) + 3)/
         Sqrt[Pi]*
        Gamma[D[Log[xg[w, a, b, NN]],
            Log[1/(((4*qbar))/((4*qbar - MJpsi^2) + w^2))]] + 5/2]/
         Gamma[D[Log[xg[w, a, b, NN]],
            Log[1/(((4*qbar))/((4*qbar - MJpsi^2) + w^2))]] +
         
           4]))^2;(*my model function-involves a numerical \
integration over variable k*)
In[17]:= chisq5[a_,b_,NN_]:=Sum[((MasData5[[j]][[1]][[2]]-F5[MasData5[[j]][[1]][[1]],a,b,NN])/MasData5[[j]][[2]][[1]])^2,{j,1,Length[MasData5]}];(*def of my chi square function*)
In[18]:= rr=Minimize[chisq5,{a,b,NN}] (*minimise chi square wrt a,b,NN*)

 

[CAF123]

Hi phyzguy,

Ignore what I wrote below. I see that you are trying to differentiate with respect to a defined function in your expression:
D[Log[xg], Log[1/(((4*qbar))/((4*qbar - MJpsi^2) + w^2))]]

I don't see how this is ever going to work.

Oh, why wouldn't it work? The derivative is just equal to this function ->

 In[1] := D[Log[xg], Log[1/(((4*qbar))/((4*qbar - MJpsi^2) + w^2))]]

Out[1]:= (2 Log[Log[11.1111 qbar]/Log[11.1111 qo]])/(Sqrt[3] Sqrt[ Log[(-9.61 + 4 qbar + w^2)/(4 qbar)] Log[Log[11.1111 qbar]/ Log[11.1111 qo]]])

Mathematica should be able to do this. I got it to run by making the following changes:

(1) You were using the same variable k for the integration variable and the index of the MasData5 list. I suspect this isn't what you wanted.
(2) I explicitly defined q0 and qbar, rather than assigning them. This may not be necessary.
(3) I made xg, F5, and chisq5 explicit functions with delayed evaluation. (the :=) I think you need to do this in order for Minimize to work.
(4) You had alphas[k_] on the RHS of F5. This won't work. It should be alphas[k].

At any rate, it now runs, but returns all zeros for the minimum values. It may not be doing what you want, but at least it runs now.

Thanks. I input these changes but at the very last step I replaced chisq5 with chisq5[a_,b_,NN_] and it still gave me same errors. I think the returning of zeroes was that it didn't know what chisq5 was.

 

[phyzguy]

So why don't you try replacing those derivatives with the actual functions if you know what they are, and then past in your new code and let me take a look. If you do that, can you define q0 and qbar ahead of time?

 

[CAF123]

Sure, here it is:

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;

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

deriv = D[Log[xg], Log[1/(((4*qbar))/((4*qbar - MJpsi^2) + w^2))]]

alphas[k_] = 4*Pi/9/Log[k/lambda]; 
lambda = 0.09;
Ca = 3;
qbar = 2.4025;
qo = 2;

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*(deriv) + 3)/Sqrt[Pi]*
           Gamma[deriv + 5/2]/ Gamma[deriv + 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*(deriv) + 3)/Sqrt[Pi]*
        Gamma[deriv + 5/2]/ Gamma[deriv + 4]))^2;

chisq5 := 
 Sum[((MasData5[[j]][[1]][[2]] - F5[MasData5[[j]][[1]][[1]]])/
     MasData5[[j]][[2]][[1]])^2, {j, 1, Length[MasData5]}];

rr = Minimize[chisq5, {a, b, NN}];

I called the derivative 'deriv' in the above and defined qbar and qo ahead. Thanks for your help.

 

[phyzguy]

Well, I'm not sure if I'm helping, but at least I got it to give a numeric answer. You need to make sure each of the functions gives a numeric (not symbolic) answer before NIntegrate and Minimize can deal with them. See if this helps. I was just putting in random values to test each function. You would know better than I what numbers make sense. Also, I just included two of the data points so it wouldn't take so long.

In[1]:= Needs["ErrorBarPlots`"]

In[2]:= 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*)
In[3]:= gamma=5.55*^-6;
MJpsi=3.1;
alphaem=1/137;

In[6]:= 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]]]];
In[7]:= deriv=D[Log[xg],Log[1/(((4*qbar))/((4*qbar-MJpsi^2)+w^2))]]

Out[7]= (2 Ca Log[Log[qbar/lambda]/Log[qo/lambda]])/(3 Sqrt[Ca Log[(-9.61+4 qbar+w^2)/(4 qbar)] Log[Log[qbar/lambda]/Log[qo/lambda]]])
In[8]:= deriv1[ww_]:=deriv/.{Ca->3.0,qo->2.0,qbar->2.4025,lambda->0.09,w->ww};
In[9]:= deriv1[7.0]
Out[9]= 0.216839
In[10]:= deriv2 = D[(2^(2*(deriv)+3)/Sqrt[Pi]*Gamma[deriv+5/2]/Gamma[deriv+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];
In[11]:= deriv3[ww_, aa_, bb_, NNN_, kk_] := deriv2/.{Ca->3.0,qo->2.0,qbar->2.4025, lambda->0.09, w->ww, a->aa, b->bb, NN->NNN, k->kk};
In[12]:= deriv3[1.0,2.0,3.0,4.0,5.0]
Out[12]= -1.05609+3.19444 I
In[13]:= gamma=5.55*^-6;
MJpsi=3.1;
alphaem=1/137;
lambda=0.09;
Ca=3;
qbar=2.4025;
qo=2;
alphas[k_]=4*Pi/9/Log[k/lambda];

In[21]:= alphas[7.0]
Out[21]= 0.320696
In[22]:= F5[w_, a_, b_, NN_]:=3.89379*^5*1/4.5/16*4*Pi^3*MJpsi^3*gamma/12/alphaem*(NIntegrate[alphas[k]/qbar/(qbar+k)*deriv3[w,a,b,NN,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*(deriv1[w])+3)/Sqrt[Pi]*Gamma[deriv1[w]+5/2]/Gamma[deriv1[w]+4]))^2;
In[23]:= F5[1.0,2.0,3.0,4.0]
During evaluation of In[23]:= NIntegrate::ncvb: NIntegrate failed to converge to prescribed accuracy after 9 recursive bisections in k near {k} = {0.0860216}. NIntegrate obtained 165912. -25668.6 I
and 169613.28983553138` for the integral and error estimates. >>
Out[23]= 3.40175*10^13-1.0784*10^13 I
In[24]:= chisq5[a_,b_,NN_]:=Sum[((MasData5[[j]][[1]][[2]]-F5[MasData5[[j]][[1]][[1]],a,b,NN])/MasData5[[j]][[2]][[1]])^2,{j,1,2}]
In[25]:= chisq5[.1,.2,.3]
Out[25]= 376033.
In[26]:= rr=Minimize[chisq5[a,b,NN],{a,b,NN}]
Out[26]:={1.50641*10^-10, {a -> -0.305193, b -> 0.112788, NN -> -0.630418}}!

 

[CAF123]

Well, I'm not sure if I'm helping, but at least I got it to give a numeric answer. You need to make sure each of the functions gives a numeric (not symbolic) answer before NIntegrate and Minimize can deal with them. See if this helps. I was just putting in random values to test each function. You would know better than I what numbers make sense. Also, I just included two of the data points so it wouldn't take so long.

Thanks! I'm just trying to understand the change you made, however, that allowed the function to give a numeric answer? I see you made replacements of the form a->aa,b->bb,NN->NNN and k->kk. What does this replacement do?

I see that it returns best fit parameters for a,b,NN, but did you get the same errors as what I got in the attached pdf file page 3? (saved as pdf because the forums do not allow .nb extensions )

 

[phyzguy]

Thanks! I'm just trying to understand the change you made, however, that allowed the function to give a numeric answer? I see you made replacements of the form a->aa,b->bb,NN->NNN and k->kk. What does this replacement do?

I see that it returns best fit parameters for a,b,NN, but did you get the same errors as what I got in the attached pdf file page 3? (saved as pdf because the forums do not allow .nb extensions )

It's a little tricky. On the one hand, you want the functions to give a numeric value so numeric functions like NIntegrate and Minimize can work with them. On the other hand, you want to do symbolic differentiation of the functions. You also need to understand the difference between a function defined by "=" and a function defined by ":=". A function defined by "=" will be evaluated immediately, while a function defined by ":=" will be defined, but not evaluated until called for. In this sense, a function defined by ":=" is like a method definition in a more conventional language like Python. So what I did was do the symbolic differentiation with deriv = ...,
and then define a numeric function with deriv1[ww_] := ... The "./" symbol means to replace every occurrence of one symbol with another, so "/. {Ca → 3.0, qo → 2.0, qbar → 2.4025, lambda → 0.09, w → ww}" replaces the symbol qbar with the numeric value 2.4025, and replaces the symbol w with the function argument ww, and so on.

Yes, I did get the same type of errors you did. Since it still returned an answer, I interpreted this to mean that the integrand was ill-defined over part of the region. Perhaps it has the log of a negative number or something like that. I'm really not sure. You might try exploring F5 function an see where these failures occur.

Hope this helps.

 

[phyzguy]

Two other thoughts that occurred to me. First, if you have some prior knowledge of reasonable values for a, b, and NN, introducing those constraints might eliminate the problem we were seeing and/or speed the minimization. Second, I see that Mathematica has a function NMinimize, which does the minimization numerically. Since you've already made the problem numeric by using NIntegrate, NMinimize might be a better choice.

 

[CAF123]

Thanks for your comments, I understood some of the replacements but what do we achieve by replacing, for example, w->ww? It's like we just renamed the argument of the function. Likewise for a->aa, etc.

Yup I tried the NMinimize and it gave a new set of values for chi square and the three parameters, allow I had to manually set the number of iterations for it to converge. However, it still has those errors I gave in that pdf.

I'm also interested in determining a covariance matrix from this analysis, providing uncertainty estimates on the parameters and their correlations. The code I wrote for such a procedure is here but it returns some peculiar errors. Do you happen to see what's going on here? If I ask too much, sorry :P

Basically I add this on after the minimisation of the chi square

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[chisq[a_, b_, NN_], vecpar[[i]], vecpar[[k]]] /.
      rr[[2]])]];
Ucov = Inverse[Ucovinv];
Ucov // MatrixForm

rr[[2]] just replaces the parameters a,b,NN with their best fit estimates where rr = Minimise[chisq[a,b,NN],{a,b,NN}] Thanks!

 

[phyzguy]

Thanks for your comments, I understood some of the replacements but what do we achieve by replacing, for example, w->ww? It's like we just renamed the argument of the function. Likewise for a->aa, etc.

We did just rename the argument of the function. We have to somehow tell it that we want the same functional form, but to now interpret it as a delayed evaluation. Maybe there is an easier way to do this, but if so I don't know what it is. Look at the simple answer below if we don't do this. It returns a symbolic answer instead of a numeric answer, because it doesn't know to associate the argument of the new function with that variable. We could type in the function again with a new name, but this is easier.

On the other question, I'll try to look at it, but may not have time.

In[1]:= deriv = D[Exp[-w^2/qo],w]
Out[1]= -((2 E^(-(w^2/qo)) w)/qo)
In[2]:= deriv1[w_]:= deriv/.{qo->2.4}
In[3]:= deriv1[2.0]
Out[3]= -0.833333 E^(-0.416667 w^2) w
In[4]:= deriv1[ww_]:=deriv/.{qo->2.4, w->ww}
In[5]:= deriv1[2.0]
Out[5]= -0.314793