Piecewise Function using lmfit

[Sam]

Hello,

I am having trouble getting a correct fit using a piecewise function that is split up over three different intervals. I tried using if and elif statements but there was an issue with a part of the function being an array and the other part being a list... is there another way to do this without using if and elif statements???? Below is my code. Any help would be greatly appreciated! Thanks! - Sam

from numpy import sqrt, pi, exp, linspace, loadtxt, arctan, log

import math

from lmfit import Model, Minimizer

import scipy.special

import matplotlib.pyplot as plt

data = loadtxt('data.dat')

r = data[:,0]

v_tot = data[:,3]

def rotation1(r, r_d, M, r_c, rho_c):

    conv = 3.0857758e21

    Sigma_0 = (M/(2.0*pi*(r_d*r_d)))

    G = 0.0000000678

    R = r*conv

    x = (R/(2.0*r_d))

    Ve_sqrd_1 = (pi*G*Sigma_0)*(R**2/r_d)

    Ve_sqrd_2 = (scipy.special.i0(x)*scipy.special.k0(x)-scipy.special.i1(x)*scipy.special.k1(x))

    Ve_sqrd = Ve_sqrd_1*Ve_sqrd_2

    Ve = sqrt(Ve_sqrd)

    Ve_kms = Ve*0.00001

    xh = R/r_c

    c = 0.866

    e = sqrt(1-(c**2))     #c is the axis ratio

    for cookie in xh:

        if cookie < 0.3:

            xi = (sqrt((e**2-cookie**2)))

            psi = sqrt(cookie**2+e**2)

            Vh_1 = (pi*G*rho_c*(r_c**2)*c)

            Vh_2 = (-2.0/psi)*(arctan(xh)+arctan((e**2)/(cookie+(c*psi))))

            Vh_3 = (1.0/psi)*log((cookie**2+1.0)*(((xh+psi)**2)/((cookie+(c*psi))**2+e**4)))

            Vh_4 = (2.0/xi)*(arctan(xi/cookie)-arctan((c*xi)/(cookie+e**2)))

            Vh = Vh_1*(Vh_2+Vh_3+Vh_4)

            Vh_s = sqrt(Vh)

            Vh_kms = Vh_s*0.00001

        elif 0.3 < cookie < 1.2:

            dx = cookie - e

            Vh = sqrt(G*rho_c*(r_c**2))*(0.79273 + 1.18206*dx - 0.763731*(dx**2) + 0.175637*(dx**3) + 0.131187*(dx**4) - 0.125809*(dx**5) - 0.00628458*(dx**6) + 0.084855456*(dx**7) - 0.0606543127*(dx**8))

            Vh_kms = Vh*0.00001

        elif cookie > 1.2:

            xi_m = sqrt(cookie**2-e**2)

            xi_p = sqrt(cookie**2+e**2)

            Vh_1 = (pi*G*rho_c*(r_c**2)*c)

            Vh_2 = (-2.0/xi_p)*(arctan(cookie)+arctan((e**2)/(cookie+(c*xi_p))))

            Vh_3 = (1.0/xi_p)*log((cookie**2+1.0)*(((cookie+xi_p)**2)/((cookie+(c*xi_p))**2+e**4)))

            Vh_4 = (2.0/xi_m)*(log(cookie+1.0)*((cookie+xi_m)/(cookie+e**2+(c*xi_m))))

            Vh = Vh_1*(Vh_2+Vh_3+Vh_4)

            Vh_s = sqrt(Vh)

            Vh_kms = Vh_s*0.00001

    

    model = (sqrt(Ve_kms**2+Vh_kms**2))

    return model

rmod1 = Model(rotation1)

result1 = rmod1.fit(v_tot, r=r, r_d = 4.0e22, M = 9.0e43, r_c = 4.0e22, rho_c = 2.0e-24) 

print(result1.fit_report())

plt.plot(r, v_tot, 'bo', label = 'data)

plt.plot(r, result1.init_fit, 'k--', label = 'initial fit')

plt.plot(r, result1.best_fit, 'r-', label = 'best fit')

legend = plt.legend(loc='lower right')

plt.title('Data Total Fit')

plt.xlabel('Radius (kpc)')

plt.ylabel('Velocity (km/s)')

plt.show()

[Matt Newville]

Hi Sam,

On Tue, Feb 23, 2016 at 12:33 PM, Sam <samantha...@gmail.com> wrote:

Hello,

I am having trouble getting a correct fit using a piecewise function that is split up over three different intervals. I tried using if and elif statements but there was an issue with a part of the function being an array and the other part being a list... is there another way to do this without using if and elif statements???? Below is my code. Any help would be greatly appreciated! Thanks! - Sam

I'm not sure I know what the problem is off-hand -- including real error messages is always helpful.  Your script is kind of hard to read, but it looks like your intention is to iterate over a scaled value for r (cookie) and calculate a value for Vh_kms for each value, then do further processing with that.

Shouldn't Vh_kms be an array of the same length as r?    I think you'd want to evaluate each item in a list for Vh_kms, then turn that into a numpy array. 

But, as you suggest, a better approach is to remove the loop and if-elif (and since you are branching on the value of cookie this should be possible).   To do this, you can use numpy's where function, perhaps like
 

    vh_kms = numpy.zeros(len(cookie))

    vh_kms[numpy.where(cookie<0.3)]  =  scary_looking_calc1(...)
    vh_kms[numpy.where(cookie>1.2)]  =  scary_looking_calc2(...)

Hope that helps

--Matt

[Sam]

Hello Matt,

Sorry for the confusing script! Thank you for your response! The reason I didn't include an error message is that I don't get one... I know that the fit is not correct however because of the plot that is produced and because the chi^2 value should be 0 because I created mock data using the same function. Also, when I fit the three functions separately it works perfectly so I know it has to do with the if statements... I tried converting Vh_kms into an array but it wasn't in the same dimension as Ve_kms...

I have tried the numpy.where function:

def rotation1(r, r_d, M, r_c, rho_c):

    conv = 3.0857758e21

    R = conv*r   

    Ve_kms = number

    

    xh = R/r_c

    Vh_s_1 = lotsofmath1      

    Vh_s_2 = lotsofmath2

    Vh_s_3 = lotsofmath3

    Vh_kms = np.zeros(len(xh))

    Vh_kms[np.where(xh<0.3)] = Vh_s_1*0.00001

    Vh_kms[np.where((xh>0.3) and (xh<1.2))] = Vh_s_2*0.00001

    Vh_kms[np.where(xh>1.2)] = Vh_s_3*0.00001

    

    model = (sqrt(Ve_kms**2+Vh_kms**2))

    return model

but this doesn't give me a correct fit either I get 'nan' values in the fit report but again no error messages.

The issue is that I can't have all of the xh values go into all the functions because if xh is less than a certain value, i.e. if it is less than 0.3 in the function Vh_s_2, then you get something divided by zero and a nan number... I need to make it so that only the first couple xh go into the first function, then the next chunk of values go into the second function, and then the remaining functions into the remaining function... This is why an if statement would be ideal if I could create the same dimension array. 

Does this make sense?

Thanks!

-Sam

[Matt Newville]

Hi Sam,

Hm, not sure.  With a model function, you can just call that function with various input (that is, without fitting at all) and see what the results are.  It sounds like that's where the problems are.  If you're getting NaNs, a fit is not going to work very well.

--
You received this message because you are subscribed to the Google Groups "lmfit-py" group.
To unsubscribe from this group and stop receiving emails from it, send an email to lmfit-py+u...@googlegroups.com.
To post to this group, send email to lmfi...@googlegroups.com.

To view this discussion on the web visit https://groups.google.com/d/msgid/lmfit-py/8680b305-17fd-4bbc-967d-5bb4783f80ca%40googlegroups.com.

For more options, visit https://groups.google.com/d/optout.

--

--Matt Newville <newville at cars.uchicago.edu> 630-252-0431