Polynomial regression: Difference between revisions
Content added Content deleted
m (Moved to Math cat) |
(quick and dirty f90 code) |
||
Line 11: | Line 11: | ||
This task is intended as a subtask for [[Measure relative performance of sorting algorithms implementations]]. |
This task is intended as a subtask for [[Measure relative performance of sorting algorithms implementations]]. |
||
=={{header|Fortran}}== |
|||
{{libheader|LAPACK}} |
|||
<pre>module fitting |
|||
contains |
|||
function polyfit(vx, vy, d) |
|||
implicit none |
|||
integer, intent(in) :: d |
|||
real(8), dimension(d+1) :: polyfit |
|||
real(8), dimension(:), intent(in) :: vx, vy |
|||
real(8), dimension(:,:), allocatable :: X |
|||
real(8), dimension(:,:), allocatable :: XT |
|||
real(8), dimension(:,:), allocatable :: XTX |
|||
integer :: i, j |
|||
integer :: n, lda, lwork |
|||
integer :: info |
|||
integer, dimension(:), allocatable :: ipiv |
|||
real(8), 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</pre> |
|||
===Example=== |
|||
<pre>program PolynomalFitting |
|||
use fitting |
|||
implicit none |
|||
! let us test it |
|||
integer :: i |
|||
real(8), dimension(11) :: x = (/ (i,i=0,10) /) |
|||
real(8), dimension(11) :: y = (/ 1, 6, 17, 34, & |
|||
57, 86, 121, 162, & |
|||
209, 262, 321 /) |
|||
real(8), dimension(3) :: a |
|||
a = polyfit(x, y, 2) |
|||
print *,a |
|||
end program</pre> |
|||
Output (lower powers first, so this seems the opposite of the Python output): |
|||
<pre> |
|||
0.999999999999813 2.00000000000002 2.99999999999998 |
|||
</pre> |
|||
=={{header|Python}}== |
=={{header|Python}}== |
||
Revision as of 23:54, 18 December 2008
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 real(8), dimension(d+1) :: polyfit real(8), dimension(:), intent(in) :: vx, vy real(8), dimension(:,:), allocatable :: X real(8), dimension(:,:), allocatable :: XT real(8), dimension(:,:), allocatable :: XTX integer :: i, j integer :: n, lda, lwork integer :: info integer, dimension(:), allocatable :: ipiv real(8), 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 :: i real(8), dimension(11) :: x = (/ (i,i=0,10) /) real(8), dimension(11) :: y = (/ 1, 6, 17, 34, & 57, 86, 121, 162, & 209, 262, 321 /) real(8), dimension(3) :: a a = polyfit(x, y, 2) print *,a end program
Output (lower powers first, so this seems the opposite of the Python output):
0.999999999999813 2.00000000000002 2.99999999999998
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).