scikit-hep / iminuit

Jupyter-friendly Python interface for C++ MINUIT2
https://scikit-hep.org/iminuit
Other
285 stars 76 forks source link

Minuit gives Nan? #321

Closed Alexandra-Wernersson closed 5 years ago

Alexandra-Wernersson commented 5 years ago

Trying to fit some data to a function that I call Ffitr but when running the code below I only get "nan" to all values. The problem lies within the chiaf if you remove if from the return sum the code works. I suspect that it as something to do with the log part of Ffitr but then I think that python should display "invalid value encountered in log". What else could be the reason for minuit settting all parameters to nan?

P.s. sorry that my code is a bit hard to read.

import numpy as np
from iminuit import Minuit

hc=197.3269788
#Parameters to fit to:
Zsrgi=np.array([1.47, 1.5, 1.54, 1.58])

Za=np.array([ 0.9468,  0.9632,  0.9707,  0.9756])

amssr=np.array([ 0.35006,  0.34592,  0.33967,  0.3175 ,  0.31341,  0.34965,
        0.33393,  0.31033,  0.30807,  0.32813,  0.29895,  0.26546,
        0.29559,  0.29298,  0.26026,  0.29264,  0.29094,  0.29074,
        0.25928,  0.3888 ,  0.23213,  0.23284,  0.22768,  0.20825,
        0.22591,  0.20437,  0.22353,  0.2036 ,  0.18987,  0.20088,
        0.18743,  0.18801])

af=np.array([ 0.0575465 ,  0.05547301,  0.04983955,  0.05096624,  0.04737787,
        0.04716958,  0.04579672,  0.0454464 ,  0.0437895 ,  0.04202845,
        0.04717754,  0.04625286,  0.04245786,  0.04105158,  0.04052182,
        0.04077226,  0.03857616,  0.03776707,  0.03453072,  0.0414683 ,
        0.0395172 ,  0.03501315,  0.03553733,  0.03483842,  0.03207193,
        0.03173218,  0.03206222,  0.03104359,  0.03046799,  0.02715095,
        0.0273168 ,  0.02622413])

amudr=np.array([ 0.01737619,  0.01377279,  0.00808912,  0.00970136,  0.00580544,
        0.00496463,  0.00397211,  0.00383197,  0.0019517 ,  0.00102585,
        0.01227267,  0.01208733,  0.00673   ,  0.00537933,  0.00521267,
        0.00450733,  0.00328733,  0.002288  ,  0.000658  ,  0.00950714,
        0.00923442,  0.00537273,  0.00533831,  0.00534351,  0.0022987 ,
        0.00224545,  0.00126688,  0.00588165,  0.00582215,  0.00267405,
        0.00266899,  0.0013943 ])

a=np.array([ 0.0904,  0.0755,  0.0647,  0.0552])

#coefficients
afij=np.array([  1.09440000e+00,  -1.00000000e+00,  -4.73000000e-02,
        -1.90580000e+00,  -1.25000000e+00,  -2.44535000e+02,
        -1.54989000e+01,  -9.39460000e+00,  -3.45830000e+00])

def Ffitr(X,s,pp,k,B=2.61,Fc=92.8,mu=770): #fit function
    temp1 = (2*B*X)/(4*pi*Fc)**2
    temp2 = temp1*(afij[0]+afij[1]*np.log((2*B*X)/mu**2))
    temp3 = temp1**2*(afij[2]+afij[3]*np.log((2*B*X)/mu**2)+\
                   afij[4]*(np.log((2*B*X)/mu**2))**2)
    temp4 = temp1**3*(afij[5]+afij[6]*np.log((2*B*X)/mu**2)+\
                      afij[7]*(np.log((2*B*X)/mu**2))**2+\
                      afij[8]*(np.log((2*B*X)/mu**2))**3)
    return Fc*pp*(1+k*s)*(1+temp2+temp3+temp4)

def chiwork(a0,a1,a2,a3,b0,b1,b2,b3,b4,b5,b6,b7,b8,b9,b10,b11,b12,b13,b14,b15,b16,b17,b18,b19,b20,b21,b22,b23,b24,b25,b26,b27,b28,b29,b30,b31,\
            z0,z1,z2,z3,za0,za1,za2,za3,\
            am0,am1,am2,am3,am4,am5,am6,am7,am8,am9,am10,am11,am12,am13,am14,am15,am16,am17,am18,am19,am20,am21,am22,am23,am24,am25,am26,am27,am28,am29,am30,am31,\
            k,y1): #chisquare

    ana=np.array([a0,a1,a2,a3])
    zs=np.array([z0,z1,z2,z3])
    za=np.array([za0,za1,za2,za3])
    am=np.array([am0,am1,am2,am3,am4,am5,am6,am7,am8,am9,am10,am11,am12,am13,am14,am15,am16,am17,am18,am19,am20,am21,am22,am23,am24,am25,am26,am27,am28,am29,am30,am31])
    bss=np.array([b0,b1,b2,b3,b4,b5,b6,b7,b8,b9,b10,b11,b12,b13,b14,b15,b16,b17,b18,b19,b20,b21,b22,b23,b24,b25,b26,b27,b28,b29,b30,b31])
    betalat=[0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,1,2,2,2,2,2,2,2,2,3,3,3,3,3]
    chia=sum((ana-a)**2)+sum((bss*a[betalat]/hc-amssr)**2)
    chiz=sum((zs-1/Zsrgi)**2)+sum((za-Za)**2)
    chiamd=sum((amudr-(1+y1*ana[betalat]**2)*ana[betalat]*am*zs[betalat]/hc)**2)
    chiaf=sum((af-ana[betalat]/hc*Ffitr(am,bss**2-mssph**2,ana[betalat]/za[betalat],k))**2)
    #chiaB=sum(((ampir**2/(2*amudr))-ampm)**2)
    return chia+chiz+chiamd+chiaf
#minuit part
mc=Minuit(chiwork)

fminc,paramc=mc.migrad()
HDembinski commented 5 years ago

Hi, this is not a fully working example, the variable hc is undefined. If you want me to test your code, please provide a working example. But you can easily debug this yourself: just print the arguments and the return value of the chiwork on every call. If it ever returns NaN, then Minuit will stop with all parameters being NaN. There is also a keyword in the Minuit.__init__ to raise an exception if it encounters NaN, see the docs for details.

Your hunch is correct that logs should give trouble if you pass a non-positive value as the argument. It could also be that you divide by zero. Both do not raise an exception when numpy arrays are involved, but in both cases you should get a warning.

PS: Your code would become much more readable if you use the new iminuit.from_array_func interface, which allows you to pass parameters to the objective function as numpy arrays. Then you don't have to spell out all named parameters in your chiwork function. See the docs for details.

Alexandra-Wernersson commented 5 years ago

Have added the value of hc now, whenever I add a print(am) in the function all the printed values are just nan's.

Not sure how to use the iminuit.from_array_func if you could spell out and example that would be great.

Thanks for your comment!

HDembinski commented 5 years ago

Have added the value of hc now, whenever I add a print(am) in the function all the printed values are just nan's.

It is still not a working example. You should try the code you are posting.

NameError: global name 'amssr' is not defined
  File "foo.py", line 56, in chiwork
    chia=sum((ana-a)**2)+sum((bss*a[betalat]/hc-amssr)**2)

Not sure how to use the iminuit.from_array_func if you could spell out and example that would be great.

See the docs, e.g. the basic tutorial. https://iminuit.readthedocs.io/en/latest/tutorials.html

HDembinski commented 5 years ago

I assume that this has resolved itself.