# Least square fitting
For data analysis sake, I read a text file but only know it contain 3 columns with unknown lines. I first go through the data and then I allocate the memory for three array to store the column separately.
Next, I need to preform least square fitting of the form of $y=P(x)=a_0+a_1x+a_2x^2$ by minimizing the weighted mean square error given by following equation.
\begin{equation}
S=\dfrac{1}{M}\sum_{i=1}^{M}\left(\dfrac{y_i - P(x_i)}{e_i}\right)^2
\tag{1}
\end{equation}
So in order to find the best cooeficient that minimize $S$, I differentiate by $a_0$, $a_1$ and $a_2$.
\begin{equation}
\dfrac{\partial s}{\partial a_0}=\dfrac{\partial s}{\partial a_1}=\dfrac{\partial s}{\partial a_2}=0
\tag{2}
\end{equation}
Then I can get the following relation in which I can transform these in to matrix form, so it can be easier to prefrom systematically on computer.
\begin{equation*}
\sum_{i=1}^{M}{\dfrac{a_2x_i^2+a_1x_i+a_0}{e_i^2}}=\sum_{i=1}^{M}{\dfrac{y_i}{e_i^2}}
\tag{3}
\end{equation*}
\begin{equation*}
\sum_{i=1}^{M}{\dfrac{a_2x_i^3+a_1x_i^2+a_0x_i}{e_i^2}}=\sum_{i=1}^{M}{\dfrac{y_ix_i}{e_i^2}}
\tag{4}
\end{equation*}
\begin{equation*}
\sum_{i=1}^{M}{\dfrac{a_2x_i^4+a_1x_i^3+a_0x^2}{e_i^2}}=\sum_{i=1}^{M}{\dfrac{y_ix_i^2}{e_i^2}}
\tag{5}
\end{equation*}
Finally, for solving inverse matrix, I use the simple Gaussian eleimination.
## Code
```fortran
!========================================================
! Purpose: Simple data analysis
!
! Date Programer Description of change
! ==== ========= =====================
! 10/14/20 Morris Original Code
!========================================================
PROGRAM MAIN
IMPLICIT NONE
INTEGER :: nmax=-1, stat=0, i=1
REAL(8), DIMENSION(:), ALLOCATABLE :: x, y, e
REAL(8) :: coef(3), a, b, c
open(4, file='exptA.dat', status='old', action='read')
print *, "exptA:"
! determine the lines in the data file
do while(stat == 0)
nmax = nmax + 1
read(4, *, iostat=stat)
end do
write (6, *) "number of line in data file:", nmax
! allocate the memory
allocate(x(nmax))
allocate(y(nmax))
allocate(e(nmax))
rewind(4)
! read the data into x, y, e
do
read(4, *, end=10) a, b, c
x(i) = a
y(i) = b
e(i) = c
i = i + 1
end do
10 close(4)
! calculate parameter
CALL QFIT(x, y, e, nmax, coef)
! write para into a file
open(99, file='exptA.coef', status='unknown')
write(99, *) coef
END PROGRAM MAIN
SUBROUTINE QFIT(x, y, e, nmax, coef)
!========================================================
! Fit a quadratic polynomail
!
! Methods: Least square fitting
!========================================================
IMPLICIT NONE
INTEGER :: nmax
REAL(8) :: x(nmax), y(nmax), e(nmax), d(3, 3), b(3)
REAL(8) :: inv(3, 3), coef(3), error
INTEGER :: i, j, k
d = 0.d0
do i = 1, 3
do j = 1, 3
do k = 1, nmax
d(i, j) = d(i, j) + x(k) ** (i+j-2) / e(k)**2
end do
end do
end do
b = 0.d0
do i = 1, 3
do j = 1, nmax
b(i) = b(i) + y(j) * x(j) ** (i-1) / e(j)**2
end do
end do
CALL INVERSE(d, inv, 3)
coef = MATMUL(inv, b)
! calculate error
error = 0.d0
do i = 1, nmax
error = error + ((y(i) - coef(3) * x(i)**2 - &
coef(2) * x(i) - coef(1)) / e(i))**2
end do
error = error / nmax
write(*, *) "mean-square error = ", error
END SUBROUTINE QFIT
SUBROUTINE INVERSE(a,c,n)
!========================================================
! Purpose: Inverse Matrix
!
! Methods: Gaussian elimintion
!--------------------------------------------------------
! input
! a(n, n) - input nxn matrix
! c - store inverse matrix
! n - dimension of matrix
!========================================================
IMPLICIT NONE
INTEGER :: n
REAL(8) :: a(n,n), c(n,n)
REAL(8) :: L(n,n), U(n,n), b(n), d(n), x(n)
REAL(8) :: coeff
INTEGER :: i, j, k
L=0.0
U=0.0
b=0.0
! step 1: forward elimination
do k=1, n-1
do i=k+1,n
coeff=a(i,k)/a(k,k)
L(i,k) = coeff
do j=k+1,n
a(i,j) = a(i,j)-coeff*a(k,j)
end do
end do
end do
! Step 2: prepare L and U matrices
! L matrix is a matrix of the elimination coefficient
! + the diagonal elements are 1.0
do i=1,n
L(i,i) = 1.0
end do
! U matrix is the upper triangular part of A
do j=1,n
do i=1,j
U(i,j) = a(i,j)
end do
end do
! Step 3: compute columns of the inverse matrix C
do k=1,n
b(k)=1.0
d(1) = b(1)
! Step 3a: Solve Ld=b using the forward substitution
do i=2,n
d(i)=b(i)
do j=1,i-1
d(i) = d(i) - L(i,j)*d(j)
end do
end do
! Step 3b: Solve Ux=d using the back substitution
x(n)=d(n)/U(n,n)
do i = n-1,1,-1
x(i) = d(i)
do j=n,i+1,-1
x(i)=x(i)-U(i,j)*x(j)
end do
x(i) = x(i)/u(i,i)
end do
! Step 3c: fill the solutions x(n) into column k of C
do i=1,n
c(i,k) = x(i)
end do
b(k)=0.0
end do
END SUBROUTINE INVERSE
```
## Result
For exptA and exptB data file, I get the following information regarding about numbers of line in data file and the mean-square error of least square fitting.
I then plot the data file with error bar and use the coefficient to comfirm the result of least square fittin which shows a pretty decent $\sigma$.
### exptA fitting result
### exptB fitting result
## Data
### exptA
```
-0.39420E+01 -0.20366E+02 0.18148E+01
-0.38304E+01 -0.16336E+02 0.25552E+01
-0.37188E+01 -0.12371E+02 0.39173E+01
-0.36071E+01 -0.15043E+02 0.19161E+01
-0.34955E+01 -0.89962E+01 0.52464E+01
-0.33839E+01 -0.12672E+02 0.27429E+01
-0.32723E+01 -0.11751E+02 0.22792E+01
-0.31606E+01 -0.34933E+00 0.67875E+01
-0.30490E+01 -0.10032E+02 0.42028E+01
-0.29374E+01 -0.25678E+01 0.41445E+01
-0.28258E+01 -0.13226E+02 0.56073E+00
-0.27141E+01 -0.11237E+02 0.95364E+00
-0.26025E+01 -0.11934E+02 0.30547E+00
-0.24909E+01 -0.51196E+01 0.44967E+01
-0.23792E+01 -0.67579E+01 0.28415E+01
-0.22676E+01 -0.13111E+02 -0.15046E+01
-0.21560E+01 -0.62505E+00 0.43766E+01
-0.20444E+01 -0.61363E+00 0.34564E+01
-0.19328E+01 -0.67593E+01 0.74104E+00
-0.18211E+01 -0.32825E+01 0.15868E+01
-0.17095E+01 -0.36012E+01 0.12858E+01
-0.15979E+01 -0.38088E+01 0.24391E+01
-0.14863E+01 0.30390E+01 0.37374E+01
-0.13746E+01 0.10038E+02 0.58238E+01
-0.12630E+01 -0.12983E+01 0.24274E+01
-0.11514E+01 0.31505E+01 0.35237E+01
-0.10398E+01 -0.21647E+01 0.61067E+00
-0.92813E+00 0.19698E+02 0.81555E+01
-0.81650E+00 0.63760E+01 0.53351E+01
-0.70488E+00 -0.11642E+01 0.53889E+00
-0.59325E+00 -0.10838E+01 0.36327E+00
-0.48163E+00 0.11133E+02 0.45123E+01
-0.37000E+00 0.99588E+01 0.40000E+01
-0.25838E+00 0.39754E+01 0.15726E+01
-0.14675E+00 -0.10716E+01 -0.37843E+00
-0.35125E-01 0.11123E+02 0.66361E+01
0.76500E-01 0.86299E+01 0.29585E+01
0.18813E+00 0.44622E+01 0.24125E+01
0.29975E+00 0.18929E+02 0.59853E+01
0.41137E+00 0.14078E+02 0.50938E+01
0.52300E+00 0.88920E+01 0.39000E+01
0.63462E+00 0.57576E+01 0.36613E+01
0.74625E+00 -0.11310E+01 -0.20916E+01
0.85787E+00 0.10038E+02 0.30618E+01
0.96950E+00 0.87802E+01 0.37976E+01
0.10811E+01 0.42462E+01 0.13839E+01
0.11928E+01 0.73228E+01 0.28725E+01
0.13044E+01 0.87032E+01 0.47868E+01
0.14160E+01 0.10245E+02 0.26515E+01
0.15276E+01 0.63178E+01 0.20043E+01
0.16393E+01 0.13180E+02 0.51823E+01
0.17509E+01 0.11510E+02 0.25522E+01
0.18625E+01 0.94104E+01 0.49381E+01
0.19741E+01 0.36609E+01 -0.63191E-01
0.20857E+01 0.58740E+01 0.19929E+01
0.21974E+01 0.76188E+01 0.13328E+01
0.23090E+01 0.49788E+01 0.71108E+00
0.24206E+01 0.79534E+01 0.16513E+01
0.25322E+01 0.66128E+01 0.10483E+01
0.26439E+01 0.73168E+01 0.26288E+01
0.27555E+01 0.79948E+01 0.21873E+01
0.28671E+01 0.13596E+02 0.40718E+01
0.29787E+01 0.89023E+01 0.39006E+01
0.30904E+01 0.90932E+01 0.24376E+01
0.32020E+01 0.75051E+01 0.31098E+01
0.33136E+01 0.75501E+01 0.31945E+01
0.34252E+01 0.11603E+02 0.30616E+01
0.35369E+01 0.22161E+01 -0.28639E+00
0.36485E+01 0.31837E+01 0.18716E+00
0.37601E+01 0.95248E+01 0.34178E+01
0.38717E+01 0.53823E+01 0.12591E+01
0.39834E+01 0.46315E+01 0.14847E+01
0.40950E+01 0.64340E+01 0.30438E+01
0.42066E+01 0.32120E+01 0.87419E+00
0.43183E+01 0.29900E+01 0.15097E+01
0.44299E+01 0.84406E+01 0.68364E+01
0.45415E+01 -0.19381E+01 -0.21139E+01
0.46531E+01 0.37380E+01 0.28566E+01
0.47647E+01 0.60862E+01 0.37179E+01
0.48764E+01 0.56171E+01 0.29501E+01
0.49880E+01 -0.20904E+00 0.75058E+00
0.50996E+01 0.45858E+01 0.51418E+01
0.52112E+01 -0.99358E+00 0.29860E+00
0.53229E+01 -0.34724E+01 -0.51248E+00
0.54345E+01 -0.80005E+00 0.99196E+00
0.55461E+01 -0.36412E+00 0.10651E+01
0.56578E+01 -0.18928E+01 0.70685E+00
0.57694E+01 0.66534E+01 0.57708E+01
0.58810E+01 0.14227E+01 0.37067E+01
0.59926E+01 -0.12939E+01 0.22748E+01
0.61042E+01 0.22355E+00 0.31367E+01
0.62159E+01 -0.33610E+01 0.10123E+01
0.63275E+01 0.53624E+01 0.45411E+01
0.64391E+01 -0.25897E+01 0.17435E+01
0.65508E+01 -0.35644E+00 0.26448E+01
0.66624E+01 -0.11984E+02 -0.10543E+01
0.67740E+01 -0.87079E+01 0.29719E+00
0.68856E+01 -0.13085E+01 0.34274E+01
0.69972E+01 -0.70956E+01 0.28106E+01
0.71089E+01 -0.15399E+01 0.54317E+01
0.72205E+01 -0.10095E+02 0.19596E+01
0.73321E+01 -0.10985E+02 0.12425E+01
0.74437E+01 -0.85746E+01 0.41862E+01
0.75554E+01 -0.11676E+02 0.28047E+01
0.76670E+01 -0.17860E+02 -0.11153E+01
```
### exptB
```
0.94200E+01 0.12003E+02 0.25237E+01
0.98802E+01 0.22413E+02 0.32783E+01
0.10340E+02 0.18672E+02 0.33360E+01
0.10801E+02 0.30008E+02 0.47507E+01
0.11261E+02 0.20125E+02 0.19533E+01
0.11721E+02 0.19832E+02 0.74258E+00
0.12181E+02 0.21018E+02 0.11917E+01
0.12642E+02 0.27243E+02 0.32780E+01
0.13102E+02 0.18685E+02 0.11419E+01
0.13562E+02 0.28819E+02 0.18491E+01
0.14022E+02 0.23723E+02 0.37038E+00
0.14482E+02 0.23557E+02 0.10863E+01
0.14943E+02 0.34416E+02 0.33841E+01
0.15403E+02 0.36641E+02 0.34528E+01
0.15863E+02 0.23023E+02 0.17781E+00
0.16323E+02 0.27936E+02 0.49830E+00
0.16784E+02 0.32528E+02 0.39793E+01
0.17244E+02 0.33982E+02 0.47126E+01
0.17704E+02 0.36867E+02 0.16400E+01
0.18164E+02 0.25480E+02 0.27535E+01
0.18625E+02 0.30230E+02 0.37615E+00
0.19085E+02 0.30579E+02 0.23139E+01
0.19545E+02 0.31657E+02 0.27682E+01
0.20005E+02 0.42230E+02 0.38381E+01
0.20465E+02 0.42455E+02 0.49446E+01
0.20926E+02 0.35572E+02 0.45570E+01
0.21386E+02 0.38299E+02 0.15124E+01
0.21846E+02 0.31688E+02 0.52855E+00
0.22306E+02 0.43937E+02 0.24275E+01
0.22767E+02 0.46227E+02 0.21311E+01
0.23227E+02 0.48397E+02 0.19415E+01
0.23687E+02 0.42365E+02 0.18236E+00
0.24147E+02 0.47931E+02 0.12596E+01
0.24608E+02 0.48006E+02 0.46134E+01
0.25068E+02 0.49225E+02 0.13542E+01
0.25528E+02 0.49326E+02 0.17348E+01
0.25988E+02 0.54900E+02 0.31581E+01
0.26448E+02 0.42037E+02 0.84000E+00
0.26909E+02 0.56182E+02 0.35686E+01
0.27369E+02 0.52539E+02 0.28789E+01
0.27829E+02 0.45868E+02 0.11636E+01
0.28289E+02 0.50228E+02 0.21532E+01
0.28750E+02 0.50245E+02 0.32942E+01
0.29210E+02 0.50579E+02 0.13452E+00
0.29670E+02 0.53725E+02 0.52055E+01
```
Sign in with Wallet
Connect another wallet