2

Here's how I generate my data and the tried fit:

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

def f(t,a,b):
    return a*np.cos(b*t)

v = 0 
x = 0.03

t = 0
dt = 0.001

time = []
pos = []

while t<3:
    a = (-5*x)/0.1
    v = v + a*dt
    x = x + v*dt
    time.append(t)
    pos.append(x)
    t = t+dt

pop, pcov = optimize.curve_fit(f,time,pos)
print(pop)

Even when I indicate initial values for the parameters (such as 0.03 for "a" and "7" for b), the resulting fit is still way off (see below, dashed line is the fit function).

enter image description here

  • Am I using the wrong library? or have I made an obvious blunder?

Thanks for any hints.

user929304
  • 465
  • 1
  • 5
  • 21
  • Could you add a diagram of a cosine and your result? – peer Feb 23 '21 at 18:01
  • @peer Sure thing, please see the update. – user929304 Feb 23 '21 at 18:07
  • I don't know how helpful this is, and maybe you've figured that out already but looks like it's specifically struggling with estimating the value of b. My guess therefore would be that it's either something cosine related or maybe because cos is negative for half of the domain. Not sure this is helpful but that was a thought I had. – djvaroli Feb 23 '21 at 18:10
  • Sorry, asking to draw a cosine was a bad idea. Could you please print the data i.e. `pos` over time and a cosine with your estimated parameters, so we can see the difference. – peer Feb 23 '21 at 18:28
  • @peer the black dots are the plotted data of `time` and `pos` and not a generic cosine. The above code includes the data too (the two aforementioned lists), if you want to experiment. – user929304 Feb 23 '21 at 18:41
  • 1
    It works fine for me if it set initial values. I ran with `pop, pcov = optimize.curve_fit(f,time,pos,p0=(.03,7))` and got pop=[0.02999163 7.07627536]. Did you enter the initial parameters as default values for the function? – Tyberius Feb 23 '21 at 18:55
  • agree with @Tyberius , works fine if proper start values are given – mikuszefski Feb 24 '21 at 09:34
  • 3 additional things. 1: The moment you set `v!=0` you get a phase and your function does not work. 2: You may use [odeint](https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.odeint.html) for your time integration and 3: There is a good way to make an initial linear fit using the techniques described in [this post](https://stackoverflow.com/a/39436209/803359) and the link therein. – mikuszefski Feb 24 '21 at 09:41

1 Answers1

1

As Tyberius noted, you need to provide better initial values. Why is that? optimize.curve_fit uses least_squares which finds a local minimum of the cost function. I believe in your case you are stuck in such a local minimum (that is not the global minimum). If you look at your diagram, your fit is approximately y=0. (It is a bit wavy because it is a cosine) If you were to increase a a bit the error would go up, so a stays close to zero. And if you were to increase b to fit the frequency of the data, the cost function would go up as well so that one stays low as well.

If you don't provide initial values, the parameters start at 1 each so it looks like this:

plt.plot(time, pos, 'black', label="data")
a,b = 1,1
init = [a*np.cos(b*t) for t in time]
plt.plot(time, init, 'b', label="a,b=1,1")
plt.legend()
plt.show()

enter image description here

a will go down and b will stay behind. I believe the scale is an additional problem. If you normalized your data to have an amplitude of 1 the humps might be more pronounced and easier to fit. If you start with a convenient value for a, b can find its way from an initial value as low as 5:

plt.plot(time, pos, 'black', label="data")
for i in [1, 4.8, 4.9, 5]:
    pop, pcov = optimize.curve_fit(f,time,pos, p0=(0.035,i))
    a,b = pop
    fit = [a*np.cos(b*t) for t in time]
    plt.plot(time, fit, label=f"$b_0 = {i}$")

plt.legend()
plt.show()

enter image description here

peer
  • 4,171
  • 8
  • 42
  • 73