Polynomial regression: Difference between revisions

Content added Content deleted
m (Fixed lang tags.)
Line 14: Line 14:


=={{header|Ada}}==
=={{header|Ada}}==
<lang ada>with Ada.Numerics.Real_Arrays; use Ada.Numerics.Real_Arrays;
<lang ada>
with Ada.Numerics.Real_Arrays; use Ada.Numerics.Real_Arrays;


function Fit (X, Y : Real_Vector; N : Positive) return Real_Vector is
function Fit (X, Y : Real_Vector; N : Positive) return Real_Vector is
Line 26: Line 25:
end loop;
end loop;
return Solve (A * Transpose (A), A * Y);
return Solve (A * Transpose (A), A * Y);
end Fit;
end Fit;</lang>
</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 &phi;<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>=&phi;<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 &phi;<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>=&phi;<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>
<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 45: Line 42:
Put (C (1), Aft => 3, Exp => 0);
Put (C (1), Aft => 3, Exp => 0);
Put (C (2), Aft => 3, Exp => 0);
Put (C (2), Aft => 3, Exp => 0);
end Fitting;
end Fitting;</lang>
</lang>
Sample output:
Sample output:
<pre>
<pre>
Line 59: Line 55:


<!-- {{does not work with|ELLA ALGOL 68|Any (with appropriate job cards AND formatted transput statements removed) - tested with release 1.8.8d.fc9.i386 - ELLA has no FORMATted transput}} -->
<!-- {{does not work with|ELLA ALGOL 68|Any (with appropriate job cards AND formatted transput statements removed) - tested with release 1.8.8d.fc9.i386 - ELLA has no FORMATted transput}} -->
<lang algol>MODE FIELD = REAL;
<lang algol68>MODE FIELD = REAL;


MODE
MODE
Line 262: Line 258:
=={{header|Fortran}}==
=={{header|Fortran}}==
{{libheader|LAPACK}}
{{libheader|LAPACK}}
<lang fortran>
<lang fortran>module fitting
contains
module fitting
contains


function polyfit(vx, vy, d)
function polyfit(vx, vy, d)
implicit none
implicit none
integer, intent(in) :: d
integer, intent(in) :: d
integer, parameter :: dp = selected_real_kind(15, 307)
integer, parameter :: dp = selected_real_kind(15, 307)
real(dp), dimension(d+1) :: polyfit
real(dp), dimension(d+1) :: polyfit
real(dp), dimension(:), intent(in) :: vx, vy
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)
real(dp), dimension(:,:), allocatable :: X
real(dp), dimension(:,:), allocatable :: XT
real(dp), dimension(:,:), allocatable :: XTX
deallocate(ipiv)
integer :: i, j
deallocate(work)
deallocate(X)
integer :: n, lda, lwork
deallocate(XT)
integer :: info
deallocate(XTX)
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
end function
</lang>
end module</lang>


===Example===
===Example===
<lang fortran>
<lang fortran>program PolynomalFitting
use fitting
program PolynomalFitting
implicit none
use fitting
implicit none
! let us test it
integer, parameter :: degree = 2
! let us test it
integer, parameter :: degree = 2
integer, parameter :: dp = selected_real_kind(15, 307)
integer, parameter :: dp = selected_real_kind(15, 307)
integer :: i
real(dp), dimension(11) :: x = (/ (i,i=0,10) /)
integer :: i
real(dp), dimension(11) :: x = (/ (i,i=0,10) /)
real(dp), dimension(11) :: y = (/ 1, 6, 17, 34, &
real(dp), dimension(11) :: y = (/ 1, 6, 17, 34, &
57, 86, 121, 162, &
57, 86, 121, 162, &
209, 262, 321 /)
real(dp), dimension(degree+1) :: a
209, 262, 321 /)
real(dp), dimension(degree+1) :: a
a = polyfit(x, y, degree)
a = polyfit(x, y, degree)
write (*, '(F9.4)') a
write (*, '(F9.4)') a


end program
end program</lang>
</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 389: Line 381:
=={{header|J}}==
=={{header|J}}==


X=:i.#Y=:1 6 17 34 57 86 121 162 209 262 321
<lang j> X=:i.#Y=:1 6 17 34 57 86 121 162 209 262 321
Y (%. (^/x:@i.@#)) X
Y (%. (^/x:@i.@#)) X
1 2 3 0 0 0 0 0 0 0 0
1 2 3 0 0 0 0 0 0 0 0</lang>


=={{header|Octave}}==
=={{header|Octave}}==
Line 402: Line 394:


{{libheader|numpy}}
{{libheader|numpy}}
<lang python>
<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]
>>> coeffs = numpy.polyfit(x,y,deg=2)
>>> coeffs = numpy.polyfit(x,y,deg=2)
>>> coeffs
>>> coeffs
array([ 3., 2., 1.])
array([ 3., 2., 1.])</lang>
</lang>
Substitute back received coefficients.
Substitute back received coefficients.
<lang python>
<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>
</lang>
Find max absolute error.
Find max absolute error.
<lang python>
<lang python>>>> '%.1g' % max(y-yf)
'1e-013'</lang>
>>> '%.1g' % max(y-yf)
'1e-013'
</lang>


===Example===
===Example===
For input arrays `x' and `y':
For input arrays `x' and `y':
<lang python>
<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]</lang>
>>> 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')
<lang python>
>>> 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>
</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]]).


Line 523: Line 505:


=={{header|TI-89 BASIC}}==
=={{header|TI-89 BASIC}}==
<lang ti-89>DelVar x
<lang ti89b>DelVar x
seq(x,x,0,10) → xs
seq(x,x,0,10) → xs
{1,6,17,34,57,86,121,162,209,262,321} → ys
{1,6,17,34,57,86,121,162,209,262,321} → ys
Retrieved from "https://rosettacode.org/wiki/Polynomial_regression"