Polynomial regression: Difference between revisions
m
bug fix and made portable
(→Example: output format and paramenter) |
m (bug fix and made portable) |
||
Line 14:
=={{header|Fortran}}==
{{libheader|LAPACK}}
<code fortran>
contains▼
function polyfit(vx, vy, d)▼
▲<pre>module fitting
implicit none▼
▲contains
integer, parameter :: dp = selected_real_kind(15, 307)
▲ function polyfit(vx, vy, d)
▲ implicit none
▲ real(8), dimension(d+1) :: polyfit
▲ real(8), dimension(:), intent(in) :: vx, vy
real(
real(
real(
integer :: i, j
integer :: n, lda, lwork
integer :: info
integer, dimension(:), allocatable :: ipiv
real(
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
</code>
===Example===
<code fortran>
use fitting
implicit none
! let us test it
integer, parameter :: degree = 2
integer, parameter :: dp = selected_real_kind(15,
integer :: i
real(8), dimension(11) :: x = (/ (i,i=0,10) /)▼
real(
57, 86, 121, 162, &
209, 262, 321 /)
real(
a = polyfit(x, y, degree)
write (*, '(F9.4)')
end program
</code>
Output (lower powers first, so this seems the opposite of the Python output):
| |||