issue with numpy arrays while curve fitting with scipy.optimize

Question

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

You should show the error. Others identifed the issue without that, but the full error would have the problem more obvious to the rest of us. — hpaulj
– hpaulj, Commented Aug 16, 2019 at 16:01

Christian K. · Accepted Answer · 2019-08-16 17:07:08Z

2

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).

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()

edited Aug 16, 2019 at 17:07

answered Aug 16, 2019 at 13:56

Christian K.

2,84322 silver badges43 bronze badges

Sign up to request clarification or add additional context in comments.

3 Comments

Raymond gsh Over a year ago

2* was a mathematical operation, I was not trying to extend the list, ... .the issue is python does not have cos function using math.cos also broke the code so I just tried the np.cos which now I see is doing something else

Christian K. Over a year ago

The result of the multiplication of the list by a factor of 2 does not depend on whether you use math.cos or np.cos. Only when multiplying a np.array by a factor you will get what you expect, as you say: a mathermatical operation.

Raymond gsh Over a year ago

can you please share the code you used? to get this plot?

Yaroslav Kornachevskyi · Accepted Answer · 2019-08-16 13:40:30Z

0

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

answered Aug 16, 2019 at 13:40

Yaroslav Kornachevskyi

1,24814 silver badges15 bronze badges

Collectives™ on Stack Overflow

issue with numpy arrays while curve fitting with scipy.optimize

2 Answers 2

3 Comments

Comments

Your Answer

Hot Network Questions

Collectives™ on Stack Overflow

2 Answers 2

3 Comments

Comments

Your Answer

Sign up or log in

Post as a guest

Related