Search code examples
pythonnumpymatplotlibregressionpolynomials

how to use Numpy.polyfit to plot trend


Thanks for user Eduard Ilyasov help me few days ago

Now i got some result, but i hardly understood these

I was trying to calculate the trend of temperature from 1979 to 2016.

    #calculate trend
    ####numpy.ployfit
    nmon = nyr * 12
    tdum = MA.arange(0,nmon)
    ntimes, ny, nx = tempF.shape 
#ntimes is time, ny is latitude, nx is longitude 
print tempF.shape

trend = MA.zeros((ny,nx), dtype='2f')
#trend = MA.zeros((ny,nx),dtype=float)

print trend.shape

for y in range (0,ny):
    for x in range (0,nx):
        trend[y,x]= numpy.polyfit(tdum, tempF[:,y,x],1)

print trend.shape
print trend

these are Results:

(
(456, 241, 480)
(241, 480, 2)
(241, 480, 2)
[[[ 0.00854342 -1.94362879]
  [ 0.00854342 -1.94362879]
  [ 0.00854342 -1.94362879]
  ..., 
  [ 0.00854342 -1.94362879]
  [ 0.00854342 -1.94362879]
  [ 0.00854342 -1.94362879]]

 [[ 0.00824162 -1.87496781]
  [ 0.00824792 -1.87640166]
  [ 0.00825524 -1.87806702]
  ..., 
  [ 0.00822667 -1.87156749]
  [ 0.00823172 -1.87271607]
  [ 0.0082366  -1.87382615]]

 [[ 0.00767854 -1.7468679 ]
  [ 0.00769076 -1.74964726]
  [ 0.00770384 -1.75262356]
  ..., 
  [ 0.00764879 -1.74010038]
  [ 0.00765911 -1.74244869]
  [ 0.00766829 -1.74453557]]

 ..., 
 [[-0.0025295   0.57546186]
  [-0.00252633  0.57474071]
  [-0.00252274  0.57392275]
  ..., 
  [-0.00253488  0.57668549]
  [-0.00253269  0.57618785]
  [-0.00253125  0.57585901]]

 [[-0.00315533  0.71783835]
  [-0.00315261  0.71721852]
  [-0.00314936  0.71648043]
  ..., 
  [-0.00315671  0.71815109]
  [-0.00315621  0.71803892]
  [-0.00315584  0.71795386]]

 [[-0.00309109  0.7032221 ]
  [-0.00309109  0.7032221 ]
  [-0.00309109  0.7032221 ]
  ..., 
  [-0.00309109  0.7032221 ]
  [-0.00309109  0.7032221 ]
  [-0.00309109  0.7032221 ]]]

What i understood that the second value in every brackets should be the coefficients which the value of trend, but i don't understand the shape of trend. What is the meaning of first number in every[], and what is the value of trend that i should use to plot the trend map?


Solution

  • If you read the documentation for numpy.polyfit() further you will see the definition of this function

    The solution minimizes the squared error

    E = \sum_{j=0}^k |p(x_j) - y_j|^2

    in the equations:

    x[0]**n * p[0] + ... + x[0] * p[n-1] + p[n] = y[0]
    x[1]**n * p[0] + ... + x[1] * p[n-1] + p[n] = y[1]
    ...
    x[k]**n * p[0] + ... + x[k] * p[n-1] + p[n] = y[k]
    

    For your case where the trend is linear that means that trend[y,x,0] is the value of trend (also called the slope) and trend[y,x,1] is the intercept.

    For an illustration consider the following example:

    import numpy as np
    from matplotlib import pyplot as plt
    
    N = 10
    
    # create repeatable data
    np.random.seed(2023)
    
    # random x test points
    xs = np.random.random(N)
    
    # random y test points
    ys = np.random.random(N)
    
    # fit the model to a given degree, 1 in this case
    trend = np.polyfit(xs, ys, 1)
    
    # plot the scatter points
    plt.plot(xs, ys, 'o')
    
    # calculate the trendline
    trendpoly = np.poly1d(trend) 
    
    # plot the trend line
    plt.plot(xs, trendpoly(xs))
    

    enter image description here