Polynomial regression

Find an approximating polynom of known degree for a given data.

Example

For input data:

x = {0,  1,  2,  3,  4,  5,  6,   7,   8,   9,   10};
y = {1,  6,  17, 34, 57, 86, 121, 162, 209, 262, 321};

The approximating polynom is:

3 x² + 2 x + 1

Here, the polynom's coefficients are (3, 2, 1).

This task is intended as a subtask for Measure relative performance of sorting algorithms implementations.

Fortran

Library: LAPACK

module fitting
contains

  function polyfit(vx, vy, d)
    implicit none
    integer, intent(in)                   :: d
    real(8), dimension(d+1)               :: polyfit
    real(8), dimension(:), intent(in)     :: vx, vy
    
    real(8), dimension(:,:), allocatable  :: X
    real(8), dimension(:,:), allocatable  :: XT
    real(8), dimension(:,:), allocatable  :: XTX
    
    integer :: i, j
    
    integer     :: n, lda, lwork
    integer :: info
    integer, dimension(:), allocatable :: ipiv
    real(8), dimension(:), allocatable :: work
    
    n = d+1
    lda = n
    lwork = n
    
    allocate(ipiv(n))
    allocate(work(lwork))
    allocate(XT(n, size(vx)))
    allocate(X(size(vx), n))
    allocate(XTX(n, n))
    
    ! prepare the matrix
    do i = 0, d
       do j = 1, size(vx)
          X(j, i+1) = vx(j)**i
       end do
    end do
    
    XT  = transpose(X)
    XTX = matmul(XT, X)
    
    ! calls to LAPACK subs DGETRF and DGETRI
    call DGETRF(n, n, XTX, lda, ipiv, info)
    if ( info /= 0 ) then
       print *, "problem"
       return
    end if
    call DGETRI(n, XTX, lda, ipiv, work, lwork, info)
    if ( info /= 0 ) then
       print *, "problem"
       return
    end if
    
    polyfit = matmul( matmul(XTX, XT), vy)
    
    deallocate(ipiv)
    deallocate(work)
    deallocate(X)
    deallocate(XT)
    deallocate(XTX)
   
  end function
  
end module

Example

program PolynomalFitting
  use fitting
  implicit none
  
  ! let us test it
  integer                 :: i
  real(8), dimension(11)  :: x = (/ (i,i=0,10) /)
  real(8), dimension(11)  :: y = (/ 1,   6,  17,  34, &
                                    57,  86, 121, 162, &
                                    209, 262, 321 /)
  real(8), dimension(3)   :: a
  
  a = polyfit(x, y, 2)
  
  print *,a

end program

Output (lower powers first, so this seems the opposite of the Python output):

 0.999999999999813        2.00000000000002        2.99999999999998

Python

Library: numpy

>>> x = [0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10]
>>> y = [1,   6,  17,  34,  57,  86, 121, 162, 209, 262, 321]
>>> coeffs = numpy.polyfit(x,y,deg=2)
>>> coeffs
array([ 3.,  2.,  1.])

Substitute back received coefficients.

>>> yf = numpy.polyval(numpy.poly1d(coeffs), x)
>>> yf
array([   1.,    6.,   17.,   34.,   57.,   86.,  121.,  162.,  209., 262.,  321.])

Find max absolute error.

>>> '%.1g' % max(y-yf)
'1e-013'

Example

For input arrays `x' and `y':

>>> x = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9]
>>> y = [2.7, 2.8, 31.4, 38.1, 58.0, 76.2, 100.5, 130.0, 149.3, 180.0]

>>> p = numpy.poly1d(numpy.polyfit(x, y, deg=2), variable='N')
>>> print p
       2
1.085 N + 10.36 N - 0.6164

Thus we confirm once more that for already sorted sequences the considered quick sort implementation has quadratic dependence on sequence length (see Example section for Python language on Query Performance page).