How do you fit parameter values in nonlinear model using R or Python

Question

I've got a function and a set of points. The function is:

s(t) = m*g*t/k-(m*m/(k*k))*(exp(-k*t/m)-1)

where m and k are the parameters in question. I have the data for (t, s(t))

I need to find the best m and k to fit this equation to my data, and I have little programming experience.

It seemed easier in R:

  df<-read.table("freefall2.txt", header = FALSE, sep = '\t', 
               strip.white = TRUE, na.strings = "empty")
time<-c(df[1])
position<-c(df[3])
fallData<-data.frame(t=time, position=c(df[3]))
plot(fallData)
m<-60
k<-1

secs=c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20)
ft=c(16,62,138,  242,  366,  504,  652,  808,  971,
     1138, 1309, 1483, 1657, 1831, 2005, 2179, 2353, 2527, 2701, 2875)
jitter(secs)
jitter(ft)

fit<-nls(ft~(m*g*secs/k)+(m^2/k^2)*(exp(-k*secs/m)-1),start = list(k=k,m=m))

That was the best try, but it throws:

Error in nlsModel(formula, mf, start, wts) : 
  singular gradient matrix at initial parameter estimates

The Python code:

def func(t,m,k):
    g=32.174
    return m*g*t/k-(m*m/(k*k))*(np.exp(-k*t/m)-1)
def plotone():
    plt.plot(D[0], D[2], 'b-', label='data')
    popt, pcov = curve_fit(func, D[0], D[2])
    plt.plot(D[0], func(D[0], *popt), 'r-', label='fit')
    popt, pcov = curve_fit(func, D[0], D[2], bounds=(0, [120, 120]))
    plt.plot(D[0], func(D[0], *popt), 'g--', label='fit-with-bounds')
    plt.xlabel('x')
    plt.ylabel('y')
    plt.legend()
    plt.show()

It generated a linear model and did not return the needed parameter values.

Rebooted attempt in Python with extensive recipe sampling:

x = D[0]
y = D[2]
ycount = 0
fitfunc = lambda t, m, k: m*32.174*t/k-(m*m/(k*k))*(np.exp(-k*t/m)-1) # Target function
errfunc = lambda t, m, k: fitfunc(t, m, k) - y[ycount]; ycount+=1# Distance to the target function
p0 = [-15., 0.8, 0., -1.] # Initial guess for the parameters
p1, success = optimize.leastsq(errfunc, D[2], args=(x, y))
print(type(p1))
print(type(success))
print(p1, success)

This one looks more promising, but I understand maybe 1/2 of it. The parts I do not understand are what I am looking at when I plot it, how to feed it my model function, what success is, and a slew of other issues.

If someone says this will work if you do x, y, and z, I will do it but otherwise I'd like to steer clear of this one since it's more involved and I have little time.

Please help me fit these parameters to my model function. I am so confused.

score 0 · Answer 1 · answered Oct 24 '17 at 02:58

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
xdata2 = np.linspace(1,20,20)
xdata2

ydata2=([16,62,138,  242,  366,  504,  652,  808,  971,1138, 1309, 1483, 1657, 1831, 2005, 2179, 2353, 2527, 2701, 2875])

ydata2

def func2(t,m,k):
    g=32.174
    return m*g*t/k-(m*m/(k*k))*(np.exp(-k*t/m)-1)

plt.plot(xdata2, ydata2, 'b-', label='xdata2')
#plt.show()

#xdata2
#ydata2

popt2, pcov2 = curve_fit(func2, xdata2, ydata2)

plt.plot(xdata2, func2(xdata2, *popt2), 'r-', label='fit')
#plt.show()

#popt3, pcov3 = curve_fit(func2, xdata2, ydata2, bounds=(0, [3., 2.]))

#plt.plot(xdata2, func2(xdata2, *popt3), 'g--', label='fit-with-bounds')


plt.show()

pcov2

popt2

#popt3

https://docs.scipy.org/doc/scipy-0.19.1/reference/generated/scipy.optimize.curve_fit.html — Mookayama, Oct 24 '17 at 02:58

score 0 · Answer 2 · answered Oct 24 '17 at 13:03

Your equation for s depends on m and k only through their ratio. That is

s(t) = r*g*t-(r*r)*(exp(-t/r)-1)

where

r = m/k

This will give many solvers a problem, as the data doesn't, for example, distinguish between values m, k and 10*m, 10*k.

You'd be better of solving for r alone.

M Newville · Answer 3 · 2017-10-25T20:28:14.780

As @dmuir points out, m and k are completely correlated.

You might find the python lmfit module (http://lmfit.github.io/lmfit-py/) useful for this. This is slightly higher level than curve_fit with better abstraction of Parameter objects. Your example would look like this:

import numpy as np
import matplotlib.pyplot as plt
from lmfit import Model

secs = np.linspace(1, 20, 20)
ft = np.array([16, 62, 138, 242, 366, 504, 652, 808, 971, 1138, 1309, 1483,
          1657, 1831, 2005, 2179, 2353, 2527, 2701, 2875])

def func(t, m, k, g=32.174):
    mk = m / k
    return mk*g*t + (mk*mk)*(np.exp(-t/mk)-1)

model = Model(func)

# make parameters with initial values for parameters
params = model.make_params(m=-50, g=32.174, k=-1)

# you can now fix parameters or keep them varying in the fit:
params['g'].vary = False
params['m'].vary = False

# or set bounds on parameters (you might want to avoid divide by 0)
params['k'].max=-1.e-8

# now do the fit:
result = model.fit(ft, params, t=secs, nan_policy='omit')

# print results
print(result.fit_report())

# plot results:
plt.plot(secs, ft, '.', label='data' )
plt.plot(secs, result.best_fit, '--', label='fit')
plt.xlabel('secs')
plt.ylabel('ft')
plt.legend()
plt.show()

Your fit doesn't look very good to me, but this might help you be better able to explore the data and model function.

I think that it should be np.exp(-t/mk)-1 – dmuir Oct 25 '17 at 17:19 — dmuir, Oct 25 '17 at 17:19

How do you fit parameter values in nonlinear model using R or Python

3 Answers3