Cebbie
Cebbie

Reputation: 1889

Python- doing least square fitting on time series data?

I have a time series dataset pr11 (shape is (151,)) which looks like the graph below when plotted. Note the very small numbers. I want to find the average slope of the data by doing a least square fit to a straight line.

Time series of precipitation from model

I've tried two different methods from another StackExchange page to get the answer. I tried using scipy.optimize.curve_fit as below...

len = np.arange(pr11.shape[0])

def f(x, A, B):
return A*x + B

A,B = curve_fit(f,pr11,len)[0]

However, this gives me a slope (A) of 1.0, which I know isn't right, so something must be off here. The "fitted" data just ends up looking the exact same as my original data. I also tried scipy.stats...

slope, intercept, r_value, p_value, std_err = stats.linregress(len,pr11)

My slope this time is a number on the order of e-08. The issue with that is when I use the equation for a line slope*x + intercept, that number multiplies my time series data to a very low value (order e-15). Therefore when I plot the fitted line, the line is horizontal and doesn't fit my data at all.

How can I get a a fitted line for this data?

Upvotes: 0

Views: 2084

Answers (1)

Julien Spronck
Julien Spronck

Reputation: 15433

A package that I like to use for fitting is lmfit. Once you install it, you can do:

from lmfit import minimize, Parameters, Parameter, report_fit
import numpy as np

# create data to be fitted
x = np.arange(150)/100.
data = 2e-6*x-5e-7 + np.random.normal(size=len(x), scale=5e-7)

# define objective function: returns the array to be minimized
def fcn2min(params, x, data):
    """ model decaying sine wave, subtract data"""
    slope = params['slope'].value
    offset = params['offset'].value

    model = slope * x + offset
    return model - data

# create a set of Parameters
params = Parameters()
params.add('slope',   value= 1.,  min=0)
params.add('offset', value= 0.)

# do fit, here with leastsq model
result = minimize(fcn2min, params, args=(x, data))

# calculate final result
final = data + result.residual

# write error report
report_fit(result.params)

# [[Variables]]
#     slope:    2.1354e-06 +/- 9.33e-08 (4.37%) (init= 1)
#     offset:  -6.0680e-07 +/- 8.02e-08 (13.22%) (init= 0)
# [[Correlations]] (unreported correlations are <  0.100)
#     C(slope, offset)             = -0.865 

# plot results
import matplotlib.pyplot as plt
plt.plot(x, data, 'k+')
plt.plot(x, final, 'r')
plt.show()

enter image description here

Upvotes: 2

Related Questions