CODEWITHSUNDEEP
pythonnumpycurve-fitting
Raymond gsh
Raymond gsh

Reputation: 383

issue with numpy arrays while curve fitting with scipy.optimize

I am trying to fit a Cosine curve but I get an error regarding the shape of my arrays. The error only appears when I change the np.cos(phi - delta) to np.cos(2*phi - delta).

from scipy.optimize import fmin
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
import numpy as np
import math
axes = plt.gca()
axes.set_ylim([-5,0])
#bond CHARMM
def func(phi, kphi, delta):
     return kphi * (1 + np.cos(2*phi - delta))



class Dihedral:  
    def __init__(self):
        self.masses = {'H': 1, 'D': 2, 'C': 12, 'O': 16} 

    def Dihedral_fit (self,x,y):

        self.popt, pcov = curve_fit(func, x, y, p0 =(0,2.3),method='trf')
        print "Cosine fit"
        print  self.popt

        plt.plot(xdata, ydata, 'b-', label='data')

        diff=sum(abs(func(x,self.popt[0],self.popt[1])-y))/len(x)
        print "AAE is"
        print diff
        plt.plot(xdata, func(xdata, *self.popt), 'r-',label='Cos: kphi=%5.3f, delta=%5.3f' % tuple(self.popt))



if __name__ == "__main__":
    xdata = [0,15,30,45,60,75,90,105]
    ydata = [-4.24,-3.82,-3.08,-2.07,-1.04,-0.30,0,-30]
    x = np.array(xdata)
    y = np.array(ydata) 
    Harm=Dihedral()
    Harm.Dihedral_fit(x,y)

    plt.legend()
    plt.show()

any help will be grateful

Upvotes: 0

Views: 87

Answers (2)

Christian K.
Christian K.

Reputation: 2823

The problem is in the last line of Dihedral_fit. You are passing a list as first argument to func and as Yaroslav pointed out a multiplication of a list with a constant yields a list which is the repetition of the original list. So exchange xdata with x, which is a numpy array and the code should run as expected. Nevertheless, the model is not a good description of your data. You should add a fitting parameter for the period of the oscillation.

EDIT:

When allowing the period to adjust I get nice results (omitting the last point, which seems to be an outlier).

enter image description here

But the success of the fit depends very much (as always) on a good initial guess.

Here is is modified code:

from scipy.optimize import fmin
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
import numpy as np
import math
axes = plt.gca()
axes.set_ylim([-5,0])
#bond CHARMM
def func(phi, kphi, om, delta):
     return kphi * (1 + np.cos(om*phi - delta))



class Dihedral:  
    def __init__(self):
        self.masses = {'H': 1, 'D': 2, 'C': 12, 'O': 16} 

    def Dihedral_fit (self,x,y):
        self.popt, pcov = curve_fit(func, x, y, p0=(-3,.03,0), method='trf')
        print "Cosine fit"
        print  self.popt

        plt.plot(x, y, 'bo', label='data')

        diff=sum(abs(func(x,*self.popt)-y))/len(x)
        print "AAE is"
        print diff
        plt.plot(x, func(x, *self.popt), 'r-',label='Cos: kphi=%5.3f, delta=%5.3f, om=%5.3f' % tuple(self.popt))


if __name__ == "__main__":
    xdata = [0,15,30,45,60,75,90]
    ydata = [-4.24,-3.82,-3.08,-2.07,-1.04,-0.30,0]
    x = np.array(xdata)
    y = np.array(ydata) 
    Harm=Dihedral()
    Harm.Dihedral_fit(x,y)

    plt.legend()
    plt.show()

Upvotes: 2

Yaroslav Kornachevskyi
Yaroslav Kornachevskyi

Reputation: 1218

You are multiplying list by 2:

xdata = [0,15,30,45,60,75,90,105]
res = 2*xdata

and so 

res = [0, 15, 30, 45, 60, 75, 90, 105, 0, 15, 30, 45, 60, 75, 90, 105]

Either rework your code or use pandas:

import pandas as pd

...

#xdata = [0,15,30,45,60,75,90,105]
xdata = pd.Series([0,15,30,45,60,75,90,105])

Upvotes: 0

Related Questions