Polynomial regression: Difference between revisions
(<code>) |
m (<lang>) |
||
| Line 14: | Line 14: | ||
=={{header|Ada}}== |
=={{header|Ada}}== |
||
< |
<lang ada> |
||
with Ada.Numerics.Real_Arrays; use Ada.Numerics.Real_Arrays; |
with Ada.Numerics.Real_Arrays; use Ada.Numerics.Real_Arrays; |
||
| Line 27: | Line 27: | ||
return Solve (A * Transpose (A), A * Y); |
return Solve (A * Transpose (A), A * Y); |
||
end Fit; |
end Fit; |
||
</ |
</lang> |
||
The function Fit implements least squares approximation of a function defined in the points as specified by the arrays ''x''<sub>''i''</sub> and ''y''<sub>''i''</sub>. The basis φ<sub>''j''</sub> is ''x''<sup>''j''</sup>, ''j''=0,1,..,''N''. The implementation is straightforward. First the plane matrix A is created. A<sub>ji</sub>=φ<sub>''j''</sub>(''x''<sub>''i''</sub>). Then the linear problem AA<sup>''T''</sup>''c''=A''y'' is solved. The result ''c''<sub>''j''</sub> are the coefficients. Constraint_Error is propagated when dimensions of X and Y differ or else when the problem is ill-defined. |
The function Fit implements least squares approximation of a function defined in the points as specified by the arrays ''x''<sub>''i''</sub> and ''y''<sub>''i''</sub>. The basis φ<sub>''j''</sub> is ''x''<sup>''j''</sup>, ''j''=0,1,..,''N''. The implementation is straightforward. First the plane matrix A is created. A<sub>ji</sub>=φ<sub>''j''</sub>(''x''<sub>''i''</sub>). Then the linear problem AA<sup>''T''</sup>''c''=A''y'' is solved. The result ''c''<sub>''j''</sub> are the coefficients. Constraint_Error is propagated when dimensions of X and Y differ or else when the problem is ill-defined. |
||
===Example=== |
===Example=== |
||
< |
<lang ada> |
||
with Fit; |
with Fit; |
||
with Ada.Float_Text_IO; use Ada.Float_Text_IO; |
with Ada.Float_Text_IO; use Ada.Float_Text_IO; |
||
| Line 46: | Line 46: | ||
Put (C (2), Aft => 3, Exp => 0); |
Put (C (2), Aft => 3, Exp => 0); |
||
end Fitting; |
end Fitting; |
||
</ |
</lang> |
||
Sample output: |
Sample output: |
||
<pre> |
<pre> |
||
| Line 54: | Line 54: | ||
=={{header|Fortran}}== |
=={{header|Fortran}}== |
||
{{libheader|LAPACK}} |
{{libheader|LAPACK}} |
||
< |
<lang fortran> |
||
module fitting |
module fitting |
||
contains |
contains |
||
| Line 119: | Line 119: | ||
end module |
end module |
||
</ |
</lang> |
||
===Example=== |
===Example=== |
||
< |
<lang fortran> |
||
program PolynomalFitting |
program PolynomalFitting |
||
use fitting |
use fitting |
||
| Line 142: | Line 142: | ||
end program |
end program |
||
</ |
</lang> |
||
Output (lower powers first, so this seems the opposite of the Python output): |
Output (lower powers first, so this seems the opposite of the Python output): |
||
| Line 155: | Line 155: | ||
{{libheader|numpy}} |
{{libheader|numpy}} |
||
< |
<lang python> |
||
>>> x = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10] |
>>> x = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10] |
||
>>> y = [1, 6, 17, 34, 57, 86, 121, 162, 209, 262, 321] |
>>> y = [1, 6, 17, 34, 57, 86, 121, 162, 209, 262, 321] |
||
| Line 161: | Line 161: | ||
>>> coeffs |
>>> coeffs |
||
array([ 3., 2., 1.]) |
array([ 3., 2., 1.]) |
||
</ |
</lang> |
||
Substitute back received coefficients. |
Substitute back received coefficients. |
||
< |
<lang python> |
||
>>> yf = numpy.polyval(numpy.poly1d(coeffs), x) |
>>> yf = numpy.polyval(numpy.poly1d(coeffs), x) |
||
>>> yf |
>>> yf |
||
array([ 1., 6., 17., 34., 57., 86., 121., 162., 209., 262., 321.]) |
array([ 1., 6., 17., 34., 57., 86., 121., 162., 209., 262., 321.]) |
||
</ |
</lang> |
||
Find max absolute error. |
Find max absolute error. |
||
< |
<lang python> |
||
>>> '%.1g' % max(y-yf) |
>>> '%.1g' % max(y-yf) |
||
'1e-013' |
'1e-013' |
||
</ |
</lang> |
||
===Example=== |
===Example=== |
||
For input arrays `x' and `y': |
For input arrays `x' and `y': |
||
< |
<lang python> |
||
>>> x = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9] |
>>> 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] |
>>> 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') |
>>> p = numpy.poly1d(numpy.polyfit(x, y, deg=2), variable='N') |
||
>>> print p |
>>> print p |
||
2 |
2 |
||
1.085 N + 10.36 N - 0.6164 |
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 [[Query Performance|'''Example''' section for Python language on ''Query Performance'' page]]). |
Thus we confirm once more that for already sorted sequences the considered quick sort implementation has quadratic dependence on sequence length (see [[Query Performance|'''Example''' section for Python language on ''Query Performance'' page]]). |
||
Revision as of 13:41, 31 January 2009
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).