RSS category feeds

RSS site feeds

More

Trend polynomials with SciPy and NumPy PDF Print E-mail

This calculates trend polynomials degree 1 (linear trend) to 6. Example is implemented based on data pair series residing in an excel sheet - could be changed of course.

Using SciPy and NumPy for (hopefully) best performance. Results (ie. polynomial coefficients) are then store back in a new sheet, together with an estimate on just how well each polynome approximates the data. I refrained from using regression factors - a sum of absolute deviations seemed to be adequate enough.

 
#--------------------------------------------------------------------------------
# Calculate trend polynomes
#--------------------------------------------------------------------------------
from scipy import linalg, array     # load required SciPy and NumPy
import numpy as np                  # packages
#--------------------------------------------------------------------------------
# Create lists for x and y values
#--------------------------------------------------------------------------------
x_reihe = []
x_shift = []
y_reihe = []
x_avg=0.0    # contains average for x values
x_dev=0.0    # contains (not quite) deviation of x values
 
for i in range(anzfiles):
    x = float(xl.ActiveSheet.Cells(i+2,3).Value)            # get x values from somewhere
    x_avg += x / float(anzfiles)                            # keep x average current - anzfiles = number of x values
    x_reihe.append(x)                                       # extend the list
    y_reihe.append(float(xl.ActiveSheet.Cells(i+2,4).Value))   # get y values from somewhere
 
# x_dev will contain an average deviation of the x' from their average
for i in range(anzfiles):
    x_dev += abs(x_reihe[i] - x_avg)/float(anzfiles)
 
# Now adjust x-vlues for numerical stability
for i in range(anzfiles):
    x_shift.append((x_reihe[i] - x_avg)/x_dev)
 
 
#--------------------------------------------------------------------------------
# This routine calculates polynome value in a numerical stable way.
# x = x value, c = tuple of polynome coefficients
#--------------------------------------------------------------------------------
def _Poly(x, c):
    n = np.shape(c)[0]          # # of coefficients (= polynome degree + 1)
    y = c[0]
    for i in range(1, n):
        y = y * x + c[i]
    return y
 
#--------------------------------------------------------------------------------
# This routine neutralizes the effect of shifting x values,
# after polynome was calculated.
# Because coding effort increases dramatically with polynomial degree,
# I have confined myself to degree 6.
# the formulas you see here have been checked and / or calculated
# with mathematica. m resp. s are the average x resp. their deviation.
# c is the tuple of polynomial coefficients. c is copied here into "a" in
# reversed sequence, such that:
# poly(x) = a[0] + a[1]*x**1 + ... + a[n]*x**n
#--------------------------------------------------------------------------------
def _XUnshift(c, m, s):
    n = np.shape(c)[0]
    a = []
    for x in c:
        a.append(x)
 
    a.reverse()
 
    if n == 2:
        a[0] = a[0] - a[1] * m / s
        a[1] = a[1] / s
    
    elif n == 3:
        a[0] = a[0] + a[2] * m**2 / s**2 - a[1] * m / s
        a[1] = -2 * a[2] * m / s**2 + a[1] / s
        a[2] = a[2] / s**2
 
    elif n == 4:
        a[0] = a[0] - a[3] * m**3 / s**3 + a[2] * m**2 / s**2 - a[1] * m / s
        a[1] = 3 * a[3] * m**2 / s**3 - 2 * a[2] * m / s**2 + a[1] / s
        a[2] = -3 * a[3] * m / s**3 + a[2] / s**2
        a[3] = a[3] / s**3
    
    elif n == 5:
        a[0] = a[0] + a[4] * m**4 / s**4 - a[3] * m**3 / s**3 + a[2] * m**2 / s**2 - a[1] * m / s
        a[1] = -4 * a[4] * m**3 / s**4 + 3 * a[3] * m**2 / s**3 - 2 * a[2] * m / s**2 + a[1] / s
        a[2] = 6 * a[4] * m**2 / s**4 - 3 * a[3] * m / s**3 + a[2] / s**2
        a[3] = -4 * a[4] * m / s**4 + a[3] / s**3
        a[4] = a[4] / s**4
 
    elif n == 6:
        a[0] = a[0] - a[5] * m**5 / s**5 + a[4] * m**4 / s**4 - a[3] * m**3 / s**3 + a[2] * m**2 / s**2 - a[1] * m / s
        a[1] = 5 * a[5] * m**4 / s**5 - 4 * a[4] * m**3 / s**4 + 3 * a[3] * m**2 / s**3 - 2 * a[2] * m / s**2 + a[1] / s
        a[2] = -10 * a[5] * m**3 / s**5 + 6 * a[4] * m**2 / s**4 - 3 * a[3] * m / s**3 + a[2] / s**2
        a[3] = 10 * a[5] * m**2 / s**5 - 4 * a[4] * m / s**4 + a[3] / s**3
        a[4] = -5 * a[5] * m / s**5 + a[4] / s**4
        a[5] = a[5] / s**5
    
    elif n == 7:
        a[0] = a[0] + a[6] * m**6 / s**6 - a[5] * m**5 / s**5 + a[4] * m**4 / s**4 - a[3] * m**3 / s**3 + a[2] * m**2 / s**2 - a[1] * m / s
        a[1] = -6 * a[6] * m**5 / s**6 + 5 * a[5] * m**4 / s**5 - 4 * a[4] * m**3 / s**4 + 3 * a[3] * m**2 / s**3 - 2 * a[2] * m / s**2 + a[1] / s
        a[2] = 15 * a[6] * m**4 / s**6 - 10 * a[5] * m**3 / s**5 + 6 * a[4] * m**2 / s**4 - 3 * a[3] * m / s**3 + a[2] / s**2
        a[3] = -20 * a[6] * m**3 / s**6 + 10 * a[5] * m**2 / s**5 - 4 * a[4] * m / s**4 + a[3] / s**3
        a[4] = 15 * a[6] * m**2 / s**6 - 5 * a[5] * m / s**5 + a[4] / s**4
        a[5] = -6 * a[6] * m / s**6 + a[5] / s**5
        a[6] = a[6] / s**6
 
    a.reverse()
    for i in range(len(a)):
        c[i] = a[i]
 
#--------------------------------------------------------------------------------
# Calculate the matrix and result vector for linear equation system.
#--------------------------------------------------------------------------------
am = [float(anzfiles), 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.]
axy = [0., 0., 0., 0., 0., 0., 0.]
ix = 0.
 
for i in range(anzfiles):
    axy[0] = axy[0] + y_reihe[i]
    
    ix = x_shift[i]   # x**1
    axy[1] = axy[1] + ix * y_reihe[i]
    am[1] += ix
    
    ix *= x_shift[i] # x**2
    axy[2] += ix * y_reihe[i]
    am[2] += ix
    
    ix *= x_shift[i] # x**3
    axy[3] += ix * y_reihe[i]
    am[3] += ix
    
    ix *= x_shift[i] # x**4
    axy[4] += ix * y_reihe[i]
    am[4] += ix
    
    ix *= x_shift[i] # x**5
    axy[5] += ix * y_reihe[i]
    am[5] += ix
    
    ix *= x_shift[i] # x**6
    axy[6] += ix * y_reihe[i]
    am[6] += ix
    
    ix *= x_shift[i] # x**7
    am[7] += ix
    
    ix *= x_shift[i] # x**8
    am[8] += ix
    
    ix *= x_shift[i] # x**9
    am[9] += ix
    
    ix *= x_shift[i] # x**10
    am[10] += ix
    
    ix *= x_shift[i] # x**11
    am[11] += ix
    
    ix *= x_shift[i] # x**12
    am[12] += ix
 
#--------------------------------------------------------------------------------
axy.reverse()       # reverse result vector entries
#--------------------------------------------------------------------------------
AMat = np.arange(49).reshape(7,7).astype(np.float64) # define matrix
 
#--------------------------------------------------------------------------------
# ... and fill it with values.
# Matrix can be reused for all polynomial degrees. Eg. in the linear case, the
# 2*2 lower right sub matrix is used.
# Similar pertains to the result vector axy.
#--------------------------------------------------------------------------------
for i in range(7):
    for j in range(7):
        AMat[i, j] = am[12 - i - j]
 
#--------------------------------------------------------------------------------
# Calculate trend polynomes
#--------------------------------------------------------------------------------
xl.ActiveWorkBook.WorkSheets.Add()
xl.ActiveSheet.Name = "Trends"
xl.ActiveSheet.Cells(1,1).Value = 'Trend'
xl.ActiveSheet.Cells(1,2).Value = 'Deviation'
xl.ActiveSheet.Cells(2,1).Value = 'Poly deg 6'
xl.ActiveSheet.Cells(3,1).Value = 'Poly deg 5'
xl.ActiveSheet.Cells(4,1).Value = 'Poly deg 4'
xl.ActiveSheet.Cells(5,1).Value = 'Poly deg 3'
xl.ActiveSheet.Cells(6,1).Value = 'Poly deg 2'
xl.ActiveSheet.Cells(7,1).Value = 'linear'
xl.ActiveSheet.Cells(8,1).Value = 'exponential'     # not yet implemented
 
for i in range(6):
    L = linalg.cholesky(AMat[i:, i:], lower=True) # Cholesky triangular matrix
    X = linalg.cho_solve((L, True), axy[i:])  # solve equation system with it
    
    _XUnshift(X, x_avg, x_dev)  # correcting the poly coefficients
 
    # calculate the quality of this trend
    abw = 0.
    for j in range(len(x_reihe)):
        abw += abs(_Poly(x_reihe[j], X) - y_reihe[j])/float(anzfiles)
    xl.ActiveSheet.Cells(i+2,2).Value = abw         # store it in sheet
    for j in range(len(X)):
        xl.ActiveSheet.Cells(i+2,j+3).Value = X[j]  # store poly coeffs in sheet
 
xl.Activesheet.Columns(2).NumberFormat = "#.##0"    # some final formatting
for i in range(3, 9):
    xl.Activesheet.Columns(i).NumberFormat = "#.##0,0000"
 
for i in range(1, 10):
    xl.Activesheet.Columns(i).EntireColumn.Autofit
 
Last Updated ( Thursday, 13 March 2014 )
 
< Prev   Next >

Feedback

Advertisement

Comments

  • I have run your code on my own machine, but I just got the last paragraph of the... More...
  • Hello Gabriel, It seem that the linked site moved to github.com/.../weboutlook (... More...
  • There is no script to download, could you provide the right link More...
  • Interesting but doesn't really help me. I do not see what the 'xxxxx' in the Ses... More...
  • You can see more python editor comparison: sparkledge.com/.../ (http://sparkledg... More...

Login Form






Lost Password?

My prefered Python IDE

My prefered Python editor is Pyscripter from MMExperts. It is not only an editor. Pyscripter is a full Python IDE including (remote) debugging, a class browser, and all other nice helpers which a full featured IDE needs.

Do you have a script for me ?

Do you have an interesting Python script which does some really cool thing on Windows ? Please post them to this site. It`s very simple - simply copy&paste it to this form. No login is requiered.

Hint: For syntax highlighting and correct Python intendation place your code between html tags <pre> and </pre>.

My prefered web framework

My prefered web framework for developing web applications is Django. Django calls itself The web framework for perfectionists with deadlines. It is a really fast, scalable and (thanks Python) the sexiest web framework of the world.