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
<lang 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; </lang> The function Fit implements least squares approximation of a function defined in the points as specified by the arrays xi and yi. The basis φj is xj, j=0,1,..,N. The implementation is straightforward. First the plane matrix A is created. Aji=φj(xi). Then the linear problem AATc=Ay is solved. The result cj are the coefficients. Constraint_Error is propagated when dimensions of X and Y differ or else when the problem is ill-defined.
Example
<lang 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; </lang> Sample output:
1.000 2.000 3.000
Fortran
<lang fortran>
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
</lang>
Example
<lang fortran>
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
</lang>
Output (lower powers first, so this seems the opposite of the Python output):
1.0000 2.0000 3.0000
Python
<lang python> >>> 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.]) </lang> Substitute back received coefficients. <lang python> >>> yf = numpy.polyval(numpy.poly1d(coeffs), x) >>> yf array([ 1., 6., 17., 34., 57., 86., 121., 162., 209., 262., 321.]) </lang> Find max absolute error. <lang python> >>> '%.1g' % max(y-yf) '1e-013' </lang>
Example
For input arrays `x' and `y': <lang python> >>> 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] </lang>
<lang python> >>> p = numpy.poly1d(numpy.polyfit(x, y, deg=2), variable='N') >>> print p
2
1.085 N + 10.36 N - 0.6164 </lang> 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).