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_{1}x+a_2{x}^2$ by minimizing the weighted mean square error given by following equation.

$$\begin{array}{}\text{(1)}& S={\displaystyle \frac{1}{M}}\sum _{i=1}^M{\left({\displaystyle \frac{y_i-P(x_i)}{e_i}}\right)}^2\end{array}$$
So in order to find the best cooeficient that minimize

$S$, I differentiate by

$a_0$,

$a_1$ and

$a_2$.

$$\begin{array}{}\text{(2)}& {\displaystyle \frac{\partial s}{\partial a_0}}={\displaystyle \frac{\partial s}{\partial a_1}}={\displaystyle \frac{\partial s}{\partial a_2}}=0\end{array}$$
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{array}{}\text{(3)}& \sum _{i=1}^M{\displaystyle \frac{a_2{x}_i^2+a_1{x}_i+a_0}{e_i^2}}=\sum _{i=1}^M{\displaystyle \frac{y_i}{e_i^2}}\end{array}$$

$$\begin{array}{}\text{(4)}& \sum _{i=1}^M{\displaystyle \frac{a_2{x}_i^3+a_1{x}_i^2+a_0{x}_i}{e_i^2}}=\sum _{i=1}^M{\displaystyle \frac{y_i{x}_i}{e_i^2}}\end{array}$$

$$\begin{array}{}\text{(5)}& \sum _{i=1}^M{\displaystyle \frac{a_2{x}_i^4+a_1{x}_i^3+a_0{x}^2}{e_i^2}}=\sum _{i=1}^M{\displaystyle \frac{y_i{x}_i^2}{e_i^2}}\end{array}$$
Finally, for solving inverse matrix, I use the simple Gaussian eleimination.

Code

!========================================================
!   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