Search code examples
pythonscipyleast-squares

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?


Solution

  • 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