CODEWITHSUNDEEP
pythonnumpymathematical-optimizationcurve-fitting
diegus
diegus

Reputation: 1188

curve fitting with a known function numpy

I have a x and y one-dimension numpy array and I would like to reproduce y with a known function to obtain "beta". Here is the code I am using:

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit

y = array([ 0.04022493,  0.04287536,  0.03983657,  0.0393201 ,  0.03810298,
    0.0363814 ,  0.0331144 ,  0.03074823,  0.02795767,  0.02413816,
    0.02180802,  0.01861309,  0.01632699,  0.01368056,  0.01124232,
    0.01005323,  0.00867196,  0.00940864,  0.00961282,  0.00892419,
    0.01048963,  0.01199101,  0.01533408,  0.01855704,  0.02163586,
    0.02630014,  0.02971127,  0.03511223,  0.03941218,  0.04280329,
    0.04689105,  0.04960554,  0.05232003,  0.05487037,  0.05843364,
    0.05120701])

x= array([ 0.,  0.08975979,  0.17951958,  0.26927937,  0.35903916,
    0.44879895,  0.53855874,  0.62831853,  0.71807832,  0.80783811,
    0.8975979 ,  0.98735769,  1.07711748,  1.16687727,  1.25663706,
    1.34639685,  1.43615664,  1.52591643,  1.61567622,  1.70543601,
    1.7951958 ,  1.88495559,  1.97471538,  2.06447517,  2.15423496,
    2.24399475,  2.33375454,  2.42351433,  2.51327412,  2.60303391,
    2.6927937 ,  2.78255349,  2.87231328,  2.96207307,  3.05183286,
    3.14159265])

def func(x,beta):
   return 1.0/(4.0*np.pi)*(1+beta*(3.0/2*np.cos(x)**2-1.0/2))

guesses = [20]
popt,pcov = curve_fit(func,x,y,p0=guesses)

y_fit =  1/(4.0*np.pi)*(1+popt[0]*(3.0/2*np.cos(x)**2-1.0/2))

plt.figure(1)
plt.plot(x,y,'ro',x,y_fit,'k-')
plt.show()

The code works but the fitting is completely off (see picture). Any idea why?

enter image description here It looks like the formula to use contains an additional parameter, i.e. p

def func(x,beta,p):
   return p/(4.0*np.pi)*(1+beta*(3.0/2*np.cos(x)**2-1.0/2))

guesses = [20,5]
popt,pcov = curve_fit(func,x,y,p0=guesses)

y_fit =  func(angle_plot,*popt)

plt.figure(2)
plt.plot(x,y,'ro',x,y_fit,'k-')
plt.show()

print popt # [ 1.23341604  0.27362069]

In the popt which one is beta and which one is p?

enter image description here

Upvotes: 2

Views: 2537

Answers (2)

Lee
Lee

Reputation: 31040

This is perhaps not what you want but, if you are just trying to get a good fit to the data, you could use np.polyfit:

fit = np.polyfit(x,y,4)
fit_fn = np.poly1d(fit)
plt.scatter(x,y,label='data',color='r')
plt.plot(x,fit_fn(x),color='b',label='fit')
plt.legend(loc='upper left')

enter image description here

Note that fit gives the coefficient values of, in this case, a 4th order polynomial:

>>> fit
array([-0.00877534,  0.05561778, -0.09494909,  0.02634183,  0.03936857])

Upvotes: 3

CT Zhu
CT Zhu

Reputation: 54330

This is going to be as good as you can get (assuming you get the equation right as @mdurant suggested), an additional intercept term is required to further improve the fit:

def func(x,beta, icpt):
   return 1.0/(4.0*np.pi)*(1+beta*(3.0/2*np.cos(x)**2-1.0/2))+icpt

guesses = [20, 0]
popt,pcov = curve_fit(func,x,y,p0=guesses)

y_fit =  func(x, *popt)

plt.figure(1)
plt.plot(x,y,'ro', x,y_fit,'k-')
print popt #[ 0.33748816 -0.05780343]

enter image description here

Upvotes: 0

Related Questions