Polynomial regression: Difference between revisions
Content added Content deleted
(→Example: output format and paramenter) |
m (bug fix and made portable) |
||
Line 14: | Line 14: | ||
=={{header|Fortran}}== |
=={{header|Fortran}}== |
||
{{libheader|LAPACK}} |
{{libheader|LAPACK}} |
||
<code fortran> |
|||
⚫ | |||
⚫ | |||
⚫ | |||
⚫ | |||
⚫ | |||
⚫ | |||
⚫ | |||
integer, parameter :: dp = selected_real_kind(15, 307) |
|||
⚫ | |||
⚫ | |||
⚫ | |||
real(dp), dimension(:), intent(in) :: vx, vy |
|||
⚫ | |||
⚫ | |||
real( |
real(dp), dimension(:,:), allocatable :: X |
||
real( |
real(dp), dimension(:,:), allocatable :: XT |
||
real( |
real(dp), dimension(:,:), allocatable :: XTX |
||
integer :: i, j |
integer :: i, j |
||
integer :: n, lda, lwork |
integer :: n, lda, lwork |
||
integer :: info |
integer :: info |
||
integer, dimension(:), allocatable :: ipiv |
integer, dimension(:), allocatable :: ipiv |
||
real( |
real(dp), dimension(:), allocatable :: work |
||
n = d+1 |
n = d+1 |
||
lda = n |
lda = n |
||
lwork = n |
lwork = n |
||
allocate(ipiv(n)) |
allocate(ipiv(n)) |
||
allocate(work(lwork)) |
allocate(work(lwork)) |
||
allocate(XT(n, size(vx))) |
allocate(XT(n, size(vx))) |
||
allocate(X(size(vx), n)) |
allocate(X(size(vx), n)) |
||
allocate(XTX(n, n)) |
allocate(XTX(n, n)) |
||
! prepare the matrix |
! prepare the matrix |
||
do i = 0, d |
do i = 0, d |
||
do j = 1, size(vx) |
do j = 1, size(vx) |
||
X(j, i+1) = vx(j)**i |
X(j, i+1) = vx(j)**i |
||
end do |
end do |
||
end do |
end do |
||
XT = transpose(X) |
XT = transpose(X) |
||
XTX = matmul(XT, X) |
XTX = matmul(XT, X) |
||
! calls to LAPACK subs DGETRF and DGETRI |
! calls to LAPACK subs DGETRF and DGETRI |
||
call DGETRF(n, n, XTX, lda, ipiv, info) |
call DGETRF(n, n, XTX, lda, ipiv, info) |
||
if ( info /= 0 ) then |
if ( info /= 0 ) then |
||
print *, "problem" |
print *, "problem" |
||
return |
return |
||
end if |
end if |
||
call DGETRI(n, XTX, lda, ipiv, work, lwork, info) |
call DGETRI(n, XTX, lda, ipiv, work, lwork, info) |
||
if ( info /= 0 ) then |
if ( info /= 0 ) then |
||
print *, "problem" |
print *, "problem" |
||
return |
return |
||
end if |
end if |
||
polyfit = matmul( matmul(XTX, XT), vy) |
polyfit = matmul( matmul(XTX, XT), vy) |
||
deallocate(ipiv) |
deallocate(ipiv) |
||
deallocate(work) |
deallocate(work) |
||
deallocate(X) |
deallocate(X) |
||
deallocate(XT) |
deallocate(XT) |
||
deallocate(XTX) |
deallocate(XTX) |
||
end function |
end function |
||
end module |
end module |
||
</code> |
|||
===Example=== |
===Example=== |
||
<code fortran> |
|||
program PolynomalFitting |
|||
use fitting |
use fitting |
||
implicit none |
implicit none |
||
! let us test it |
! let us test it |
||
integer, parameter :: degree = 2 |
integer, parameter :: degree = 2 |
||
integer |
integer, parameter :: dp = selected_real_kind(15, 307) |
||
integer :: i |
|||
⚫ | |||
real( |
real(dp), dimension(11) :: x = (/ (i,i=0,10) /) |
||
⚫ | |||
57, 86, 121, 162, & |
57, 86, 121, 162, & |
||
209, 262, 321 /) |
209, 262, 321 /) |
||
real( |
real(dp), dimension(degree+1) :: a |
||
a = polyfit(x, y, degree) |
a = polyfit(x, y, degree) |
||
write (*, '(F9.4)') |
write (*, '(F9.4)') a |
||
end program |
end program |
||
</code> |
|||
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): |
Revision as of 12:34, 27 January 2009
Polynomial regression
You are encouraged to solve this task according to the task description, using any language you may know.
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.
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).