***LEAST SQUARE REGRESSION**********************************************
*                  PROGRAM FOR LEAST SQUARE REGRESSION                 *
************************************************************************
*
      program LEASTREGRESS
      implicit doubleprecision(a-h,o-z)
      parameter (nd = 40)
      doubleprecision XM(nd),YM(nd)
      integer option
*
      sx = 0
      sy = 0
      x2 = 0
      xy = 0
      st = 0
      sr = 0
*
      print*, 'The relationship between "X & Y" data:'
      write(6,101) '[1]  LINEAR','[2]  EXPONENTIAL','[3]  POWER'
101   format(''/a32/a37/a31)
      write(6,*)
      write(6,*)
      write(6,*)' Select < 1 / 2 / 3 >:'
      read(*,*) option
      print*, 'Enter the number of "X & Y" data pairs, n:'
      read(*,*) n
      print*, 'Enter "X" data:'
      read(*,*) (XM(i),i = 1,n)
      print*, 'Enter "Y" data:'
      read(*,*) (YM(i),i = 1,n)
*
      if (option == 2) then
         do i = 1, n
            YM(i) = dlog(YM(i))
         enddo
      elseif (option == 3) then
         do i = 1, n
            YM(i) = dlog10(YM(i))
            XM(i) = dlog10(XM(i))
         enddo
      endif
*
      do i = 1, n
         sx = sx + XM(i)
         sy = sy + YM(i)
         x2 = x2 + (XM(i)**2)
         xy = xy + XM(i) *YM(i)
      enddo
      xmean = sx/n
      ymean = sy/n
      a1 = (n*xy - sx*sy)/(n*x2 - sx**2)
      a0 = ymean - (a1 *xmean)
*
      do i = 1, n
         st = st + (YM(i) - ymean)**2
         sr = sr + (YM(i) - a0 - a1 *XM(i))**2
      enddo
*
      sd = dsqrt(st/(n - 1))
      if (n.ne.2) then
         serror = dsqrt(sr/(n - 2))
      endif
      perror = (st - sr)/st *100
      correl = dsqrt(dabs(perror))/10d0
*
      write(6,*)'------------------------------------------------------'
      if (option == 1) then
         write(6,102) 'The least square fit is:    Y   =',a0,'+',a1,'X'
      elseif (option == 2) then
         write(6,103) 'The least square fit is:    ln Y  = ',a1,'X +',a0
         c0 = dexp(a0)
         write(6,104) 'The exponential equation for the fit is:   Y  = '
     $,c0,'exp(',a1,'X )'
      elseif (option == 3) then
         write(6,105) 'The least square fit is:  log Y  = ',a1,' log X +
     $',a0
         c0 = 10**a0
         write(6,106) 'The power equation for the fit is:  Y  = ',c0,'X 
     $**',a1
      endif
*
102   format('',a33,f9.4,3x,a1,f8.4,a2)
103   format('',a36,f9.4,1x,a3,f8.4)
104   format('',a48,f9.4,a5,f7.3,1x,a3)
105   format('',a35,f9.4,a9,f8.4)
106   format('',a40,f8.4,a5,f7.4)
*
      write(6,*)
      write(6,*) 'The total standard deviation:',sd
      if (n.ne.2) then
         write(6,*)'The standard error of estimate:',serror
      endif
      write(6,*)'The correlation coefficient:',correl
      write(6,*)'Percentage of error reduction due to straight line fit:
     $',perror,' %'
      write(6,*)'------------------------------------------------------'
*
      stop
      end