polynomial curve fitting {Fortran 90 / Python} [algorithm] by 바죠

Least Squares Fitting--Polynomial (최소 자승 다항식 맞추기)


http://mathworld.wolfram.com/LeastSquaresFittingPolynomial.html

Polynomial regression (다항식 회귀)
http://en.wikipedia.org/wiki/Polynomial_regression
Curve fitting (곡선 맞추기)
http://en.wikipedia.org/wiki/Curve_fitting

다항함수로 피팅하여 데이터 표현하기.
아래와 같은 {x,y} 형식의 데이터가 y= a0 + a1 *x + a2 *x**2 +......와 같은 다항식으로 적절히 표현할 수 있다고 생각할 때, 해당 에러와 a0, a1, a2,.... 계수들을 구하는 문제이다.
아래에서 처럼 포트란 90 또는 Python의 Numpy 모듈을 사용하면 쉽게 위의 문제를 해결할 수 있다.
구체적인 의미와 이론은 참고 URL을 통해서 얻을수 있다.

cat inp1.dat
0.0d0   51.49866347d0
0.1d0   23.92508775d0
0.2d0   10.35390700d0
0.3d0    2.58426990d0
0.4d0   -2.18351656d0
0.5d0   -5.41745387d0
0.6d0   -7.62452181d0
0.7d0   -9.25455804d0
0.8d0  -10.45592989d0
0.9d0  -11.39244138d0
1.0d0  -12.12433704d0

 polfit.x
 Enter the number of points, _ndata:
11
 Enter the file name, fname:
inp1.dat
 Enter the _ndeg :
10
February 22 2012  11:06:17.527 PM
  X matrix :           11 by           11
    11.00000000000000         5.500000000000000         3.850000000000000     
    3.025000000000000         2.533300000000000         2.208250000000000     
    1.978405000000000         1.808042500000000         1.677313330000000     
    1.574304985000000         1.491434192500000         5.500000000000000     
    3.850000000000000         3.025000000000000         2.533300000000000     
    2.208250000000000         1.978405000000000         1.808042500000000     
    1.677313330000000         1.574304985000000         1.491434192500000     
    1.423643196250000         3.850000000000000         3.025000000000000     
    2.533300000000000         2.208250000000000         1.978405000000000     
    1.808042500000000         1.677313330000000         1.574304985000000     
    1.491434192500000         1.423643196250000         1.367428536133000     
    3.025000000000000         2.533300000000000         2.208250000000000     
    1.978405000000000         1.808042500000000         1.677313330000000     
    1.574304985000000         1.491434192500000         1.423643196250000     
    1.367428536133000         1.320286076114500         2.533300000000000     
    2.208250000000000         1.978405000000000         1.808042500000000     
    1.677313330000000         1.574304985000000         1.491434192500000     
    1.423643196250000         1.367428536133000         1.320286076114500     
    1.280378029534450         2.208250000000000         1.978405000000000     
    1.808042500000000         1.677313330000000         1.574304985000000     
    1.491434192500000         1.423643196250000         1.367428536133000     
    1.320286076114500         1.280378029534450         1.246324856379625     
    1.978405000000000         1.808042500000000         1.677313330000000     
    1.574304985000000         1.491434192500000         1.423643196250000     
    1.367428536133000         1.320286076114500         1.280378029534450     
    1.246324856379625         1.217070613200973         1.808042500000000     
    1.677313330000000         1.574304985000000         1.491434192500000     
    1.423643196250000         1.367428536133000         1.320286076114500     
    1.280378029534450         1.246324856379625         1.217070613200973     
    1.191793189353773         1.677313330000000         1.574304985000000     
    1.491434192500000         1.423643196250000         1.367428536133000     
    1.320286076114500         1.280378029534450         1.246324856379625     
    1.217070613200973         1.191793189353773         1.169842891165485     
    1.574304985000000         1.491434192500000         1.423643196250000     
    1.367428536133000         1.320286076114500         1.280378029534450     
    1.246324856379625         1.217070613200973         1.191793189353773     
    1.169842891165485         1.150699451020125         1.491434192500000     
    1.423643196250000         1.367428536133000         1.320286076114500     
    1.280378029534450         1.246324856379625         1.217070613200973     
    1.191793189353773         1.169842891165485         1.150699451020125     
    1.133941318588326    


            0  ierror, normal return from findinv
 ai; i=0, 1, 2,...., k
    51.48422225565688        -436.9223538385205         2669.208069765820     
   -15779.54914660206         70748.66086736298        -214985.7450545331     
    430337.0343733301        -557305.2334284838         448565.3468817059     
   -203744.9585273829         39869.01621239795    

   5.0000000000000000E-003  sec
  xmgrace fort.11 fort.12

xmgrace : http://incredible.egloos.com/4486070



----------------------------------------------------------------------




!234567890
       module polfit
       implicit none
       private
       save
!      Written by In-Ho Lee (KIAS/KRISS) 2010.
!      Reference : http://mathworld.wolfram.com/LeastSquaresFittingPolynomial.html
       integer ndeg,ndata
       integer npt
       real*8, allocatable :: xdata(:),ydata(:)
       real*8, allocatable :: ak(:),xmat(:,:),vec2(:),zmat(:,:),amat(:,:),bmat(:,:)

       public :: initial, final, polfit_do

       contains

       subroutine initial(_ndeg,_ndata,_npt,fname)
       implicit none
       character*60 fname
       integer _ndeg,_ndata,_npt
       integer k1,i

       ndeg=_ndeg
       ndata=_ndata
       npt=_npt
       if(_npt <= 0) npt=30*ndata

       allocate(xdata(ndata),ydata(ndata))

       open(31,file=trim(fname),form='formatted')
       do i=1,ndata
       read(31,*) xdata(i),ydata(i)
       enddo
       close(31)

       allocate(ak(0:ndeg))
       k1=ndeg+1
       allocate(xmat(ndata,k1)) ; allocate(zmat(k1,ndata))
       allocate(amat(k1,k1)) ; allocate(bmat(k1,k1)) ; allocate(vec2(k1))
       end subroutine initial

       subroutine final()
       implicit none
       deallocate(vec2) ; deallocate(bmat,amat,zmat,xmat) ; deallocate(ak)
       deallocate(xdata,ydata)
       end subroutine final

       subroutine polfit_do()
       implicit none
       integer i,j,k1,ierror
       real*8 tmp,tmq,tmr
!
       open(11,file='fort.11',form='formatted')
       do i=1,ndata
       write(11,'(2f18.8)') xdata(i),ydata(i)
       enddo
       close(11)
!
       k1=ndeg+1
       xmat=0.d0
       do i=1,ndata
       do j=1,k1
       if(j == 1)then
       xmat(i,j)=1.d0
                 else
       xmat(i,j)=xdata(i)**(j-1)
                 endif
       enddo
       enddo
       zmat=transpose(xmat)
       amat=matmul(zmat,xmat)
       print*, ' X matrix :', (ndeg+1), 'by', (ndeg+1)
       print*, amat
       print*

       call findinv(amat, bmat, k1, ierror)
       if(ierror == 0)then
       print*, ierror,' ierror, normal return from findinv'
                      else
       print*, ierror,' ierror, abnormal return from findinv'
                      endif
       vec2=matmul(matmul(bmat,zmat),ydata)

       print*, 'ai; i=0, 1, 2,...., k'
       print*, vec2
       print*
       do i=0,ndeg
       ak(i)=vec2(i+1)
       enddo
!
       open(12,file='fort.12',form='formatted')
       tmr=(xdata(ndata)-xdata(1))/float(npt-1)
       do i=1,npt
       tmp=xdata(1)+tmr*float(i-1)
       tmq=ak(0)
       do j=1,ndeg
       tmq=tmq+ak(j)*tmp**j
       enddo
       write(12,'(2f18.8)') tmp,tmq
       enddo
       close(12)
!
       end subroutine polfit_do
       end module polfit



       program test
       USE polfit, ONLY : initial,final,polfit_do
       implicit none
       integer _ndeg,_ndata,_npt
       character*60 fname
       real*8 tt1,tt2


       _ndeg=3
       _ndeg=5
       _ndeg=10
       _npt=0
       _ndata=16
       fname='cv.dat'

       write(6,*) 'Enter the number of points, _ndata: '
       read(5,*) _ndata
       write(6,*) 'Enter the file name, fname: '
       read(5,*) fname
       write(6,*) 'Enter the _ndeg : '
       read(5,*) _ndeg

       call timestamp()
       call real_time (tt1)

       call initial(_ndeg,_ndata,_npt,fname)
       call polfit_do()
       call final()

       call real_time (tt2)
       write(6,*) tt2-tt1,' sec'
       write(6,*) ' xmgrace fort.11 fort.12'
!      call timestamp()
       end program test
!Reference : Algorithm has been well explained in:
!http://math.uww.edu/~mcfarlat/inverse.htm          
!http://www.tutor.ms.unimelb.edu.au/matrix/matrix_inverse.html
        SUBROUTINE FINDInv(matrix, inverse, n, errorflag)
        IMPLICIT NONE
        !Declarations
        INTEGER, INTENT(IN) :: n
        INTEGER, INTENT(OUT) :: errorflag  !Return error status. -1 for error, 0 for normal
        REAL*8, INTENT(IN), DIMENSION(n,n) :: matrix  !Input matrix
        REAL*8, INTENT(OUT), DIMENSION(n,n) :: inverse !Inverted matrix

        LOGICAL :: FLAG = .TRUE.
        INTEGER :: i, j, k, l
        REAL*8 :: m
        REAL*8, DIMENSION(n,2*n) :: augmatrix !augmented matrix

        !Augment input matrix with an identity matrix
        DO i = 1, n
                DO j = 1, 2*n
                        IF (j <= n ) THEN
                                augmatrix(i,j) = matrix(i,j)
                        ELSE IF ((i+n) == j) THEN
                                augmatrix(i,j) = 1
                        Else
                                augmatrix(i,j) = 0
                        ENDIF
                END DO
        END DO

        !Reduce augmented matrix to upper traingular form
        DO k =1, n-1
                IF (augmatrix(k,k) == 0) THEN
                        FLAG = .FALSE.
                        DO i = k+1, n
                                IF (augmatrix(i,k) /= 0) THEN
                                        DO j = 1,2*n
                                                augmatrix(k,j) = augmatrix(k,j)+augmatrix(i,j)
                                        END DO
                                        FLAG = .TRUE.
                                        EXIT
                                ENDIF
                                IF (FLAG .EQV. .FALSE.) THEN
                                        PRINT*, "Matrix is non - invertible"
                                        inverse = 0
                                        errorflag = -1
                                        return
                                ENDIF
                        END DO
                ENDIF
                DO j = k+1, n                  
                        m = augmatrix(j,k)/augmatrix(k,k)
                        DO i = k, 2*n
                                augmatrix(j,i) = augmatrix(j,i) - m*augmatrix(k,i)
                        END DO
                END DO
        END DO
       
        !Test for invertibility
        DO i = 1, n
                IF (augmatrix(i,i) == 0) THEN
                        PRINT*, "Matrix is non - invertible"
                        inverse = 0
                        errorflag = -1
                        return
                ENDIF
        END DO
       
        !Make diagonal elements as 1
        DO i = 1 , n
                m = augmatrix(i,i)
                DO j = i , (2 * n)                             
                           augmatrix(i,j) = (augmatrix(i,j) / m)
                END DO
        END DO
       
        !Reduced right side half of augmented matrix to identity matrix
        DO k = n-1, 1, -1
                DO i =1, k
                m = augmatrix(i,k+1)
                        DO j = k, (2*n)
                                augmatrix(i,j) = augmatrix(i,j) -augmatrix(k+1,j) * m
                        END DO
                END DO
        END DO                         
       
        !store answer
        DO i =1, n
                DO j = 1, n
                        inverse(i,j) = augmatrix(i,j+n)
                END DO
        END DO
        errorflag = 0
        END SUBROUTINE FINDinv
subroutine real_time ( seconds )

!*******************************************************************************
!
!! REAL_TIME returns the real time in seconds.
!
!  Modified:
!
!    07 November 2006
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Output, real ( kind = 8 ) SECONDS, a reading of the real time clock,
!    in seconds.
!
  implicit none

  integer clock_count
  integer clock_max
  integer clock_rate
  real ( kind = 8 ) seconds

  call system_clock ( clock_count, clock_rate, clock_max )

  seconds = real ( clock_count, kind = 8 ) &
    / real ( clock_rate, kind = 8 )

  return
end
!! TIMESTAMP prints the current YMDHMS date as a time stamp.
!
!  Example:
!
!    May 31 2001   9:45:54.872 AM
!
!  Modified:
!
!    31 May 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    None
!
  implicit none

  character ( len = 8 ) ampm
  integer d
  character ( len = 8 ) date
  integer h
  integer m
  integer mm
  character ( len = 9 ), parameter, dimension(12) :: month = (/ &
    'January  ', 'February ', 'March    ', 'April    ', &
    'May      ', 'June     ', 'July     ', 'August   ', &
    'September', 'October  ', 'November ', 'December ' /)
  integer n
  integer s
  character ( len = 10 )  time
  integer values(8)
  integer y
  character ( len = 5 ) zone

  call date_and_time ( date, time, zone, values )

  y = values(1)
  m = values(2)
  d = values(3)
  h = values(5)
  n = values(6)
  s = values(7)
  mm = values(8)

  if ( h < 12 ) then
    ampm = 'AM'
  else if ( h == 12 ) then
    if ( n == 0 .and. s == 0 ) then
      ampm = 'Noon'
    else
      ampm = 'PM'
    end if
  else
    h = h - 12
    if ( h < 12 ) then
      ampm = 'PM'
    else if ( h == 12 ) then
      if ( n == 0 .and. s == 0 ) then
        ampm = 'Midnight'
      else
        ampm = 'AM'
      end if
    end if
  end if

  write ( *, '(a,1x,i2,1x,i4,2x,i2,a1,i2.2,a1,i2.2,a1,i3.3,1x,a)' ) &
    trim ( month(m) ), d, y, h, ':', n, ':', s, '.', mm, trim ( ampm )

  return
end


----------------------------------------------------------------------



Python에서 사용할 수 있는 Numpy 모듈을 활용하면 아주 손쉽게 폴리노미얼 피팅을 할 수 있다.
아래와 같이 입력을 하면 원하는 다항식을 얻어낼 수 있다.

Numpy : http://incredible.egloos.com/3607490



 python
Python 2.7.1 (r271:86832, Apr 11 2011, 16:56:03)
[GCC 4.1.2 20080704 (Red Hat 4.1.2-46)] on linux2
Type "help", "copyright", "credits" or "license" for more information.
>>> import numpy
>>> x=[0,1,2,3,4,5,6,7,8,9]
>>> y=[1.1,3.9,11.2,21.5,34.8,51.,70.2, 92.3, 117.4, 145.5]
>>> p=numpy.poly1d(numpy.polyfit(x,y,deg=2))
>>> print p
       2
1.517 x + 2.483 x + 0.4927
>>> x=[0.,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0]
>>> y=[51.4986,23.925,10.3539,2.5842,-2.1835,-5.4174,-7.6245,-9.2545,-10.04559,-11.3924,-12.124]
>>> p=numpy.poly1d(numpy.polyfit(x,y,deg=2))
>>> p=numpy.poly1d(numpy.polyfit(x,y,deg=10))
>>> print p
           10             9             8            7             6
9.332e+04 x  - 4.558e+05 x + 9.565e+05 x - 1.13e+06 x + 8.262e+05 x
              5             4             3        2
 - 3.883e+05 x + 1.184e+05 x - 2.366e+04 x + 3372 x - 463.2 x + 51.5
>>>



----------------------------------------------------------------------


Online regression tool:
http://www.xuru.org/rt/PR.asp


----------------------------------------------------------------------


    0.000000000000000        1.5413331334360181E-003
   5.2631578947368418E-002   1.3815039801161721E-002
   0.1052631578947368        3.8060233744356631E-002
   0.1578947368421053        7.3679917822953855E-002
   0.2105263157894737        0.1197970171999845    
   0.2631578947368421        0.1752759758349081    
   0.3157894736842105        0.2387507176420255    
   0.3684210526315789        0.3086582838174551    
   0.4210526315789473        0.3832773180720472    
   0.4736842105263158        0.4607704521360775    
   0.5263157894736842        0.5392295478639224    
   0.5789473684210527        0.6167226819279527    
   0.6315789473684210        0.6913417161825448    
   0.6842105263157894        0.7612492823579744    
   0.7368421052631579        0.8247240241650917    
   0.7894736842105263        0.8802029828000154    
   0.8421052631578947        0.9263200821770461    
   0.8947368421052631        0.9619397662556434    
   0.9473684210526315        0.9861849601988383    
    1.000000000000000        0.9984586668665640
    

!234567890
       implicit none
       integer nn
       integer i
       real*8 aa,bb
       real*8 pi
       real*8 tmp,tmr,du
       aa=0.d0
       bb=1.d0
       pi=4.d0*atan(1.d0)
       nn=20
       du=(bb-aa)/float(nn-1)
       do i=1,nn
       tmr=aa+du*(i-1)
       tmp=0.5d0*(aa+bb)-0.5d0*(bb-aa)*cos(pi*(2.d0*i-1.d0)/(2.d0*nn))
       write(6,*) tmr,tmp
       enddo
       stop
       end




----------------------------------------------------------------------
----------------------------------------------------------------------


핑백

덧글

댓글 입력 영역

최근 포토로그