Polynomial regression: Difference between revisions
m (Taking "Example" off the TOC) |
(<code>) |
||
| Line 14: | Line 14: | ||
=={{header|Ada}}== |
=={{header|Ada}}== |
||
<ada> |
<code 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; |
||
</ |
</code> |
||
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=== |
||
<ada> |
<code 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; |
||
</ |
</code> |
||
Sample output: |
Sample output: |
||
<pre> |
<pre> |
||
| Line 155: | Line 155: | ||
{{libheader|numpy}} |
{{libheader|numpy}} |
||
<code python> |
|||
| ⚫ | |||
>>> x = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10] |
|||
| ⚫ | |||
| ⚫ | |||
>>> coeffs = numpy.polyfit(x,y,deg=2) |
|||
>>> coeffs |
|||
array([ 3., 2., 1.]) |
|||
</code> |
|||
Substitute back received coefficients. |
Substitute back received coefficients. |
||
<code python> |
|||
>>> yf = numpy.polyval(numpy.poly1d(coeffs), x) |
|||
>>> yf |
|||
array([ 1., 6., 17., 34., 57., 86., 121., 162., 209., 262., 321.]) |
|||
</code> |
|||
Find max absolute error. |
Find max absolute error. |
||
<code python> |
|||
>>> '%.1g' % max(y-yf) |
|||
'1e-013' |
|||
</code> |
|||
===Example=== |
===Example=== |
||
For input arrays `x' and `y': |
For input arrays `x' and `y': |
||
<code 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] |
|||
</code> |
|||
<code python> |
|||
>>> p = numpy.poly1d(numpy.polyfit(x, y, deg=2), variable='N') |
|||
| ⚫ | |||
>>> print p |
|||
2 |
|||
2 |
|||
1.085 N + 10.36 N - 0.6164 |
|||
</code> |
|||
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 17:18, 27 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
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;
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
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;
Sample output:
1.000 2.000 3.000
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
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
>>> 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).