Polynomial regression

Polynomial regression
You are encouraged to solve this task according to the task description, using any language you may know.

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 x2 + 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.

Ada

<ada> with Ada.Numerics.Real_Arrays; use Ada.Numerics.Real_Arrays;

function Fit (X, Y : Real_Vector; N : Positive) return Real_Vector is

  A : Real_Matrix (0..N, X'Range);  -- The plane

begin

  for I in A'Range (2) loop
     for J in A'Range (1) loop
        A (J, I) := X (I)**J;
     end loop;
  end loop;
  return Solve (A * Transpose (A), A * Y);

end Fit; </ada> The function Fit implements least squares approximation of a function defined in the points as specified by the arrays x_i and y_i. The basis φ_j is x^j, j=0,1,..,N. The implementation is straightforward. First the plane matrix A is created. A_ji=φ_j(x_i). Then the linear problem AA^Tc=Ay is solved. The result c_j are the coefficients. Constraint_Error is propagated when dimensions of X and Y differ or else when the problem is ill-defined.

Example

<ada> with Fit; with Ada.Float_Text_IO; use Ada.Float_Text_IO;

procedure Fitting is

  C : constant Real_Vector :=
         Fit
         (  (0.0, 1.0,  2.0,  3.0,  4.0,  5.0,   6.0,   7.0,   8.0,   9.0,  10.0),
            (1.0, 6.0, 17.0, 34.0, 57.0, 86.0, 121.0, 162.0, 209.0, 262.0, 321.0),
            2
         );

begin

  Put (C (0), Aft => 3, Exp => 0);
  Put (C (1), Aft => 3, Exp => 0);
  Put (C (2), Aft => 3, Exp => 0);

end Fitting; </ada> Sample output:

 1.000 2.000 3.000

Fortran

Library: LAPACK

module fitting
contains

  function polyfit(vx, vy, d)
    implicit none
    integer, intent(in)                   :: d
    integer, parameter                    :: dp = selected_real_kind(15, 307)
    real(dp), dimension(d+1)              :: polyfit
    real(dp), dimension(:), intent(in)    :: vx, vy
   
    real(dp), dimension(:,:), allocatable :: X
    real(dp), dimension(:,:), allocatable :: XT
    real(dp), dimension(:,:), allocatable :: XTX
   
    integer :: i, j
   
    integer     :: n, lda, lwork
    integer :: info
    integer, dimension(:), allocatable :: ipiv
    real(dp), 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, parameter      :: degree = 2
  integer, parameter      :: dp = selected_real_kind(15, 307)
  integer                 :: i
  real(dp), dimension(11) :: x = (/ (i,i=0,10) /)
  real(dp), dimension(11) :: y = (/ 1,   6,  17,  34, &
                                   57,  86, 121, 162, &
                                   209, 262, 321 /)
  real(dp), dimension(degree+1) :: a
 
  a = polyfit(x, y, degree)
 
  write (*, '(F9.4)') a

end program

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

   1.0000
   2.0000
   3.0000

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).