C Linear regression utility
C Written in a moment of weakness by Paul F. Johnson

C Source released under the eduware licence. 
C This is similar to freeware, except if used in an academic 
C environment, I'd like to know. 
C The source may not be distributed on a cover disc or CD or any 
C other form of magnetic media. This excludes the internet. 
C If diseminated on the internet, please let me know. 
C 
C email paul@acornusers.org to "register".

      PROGRAM LINREG
      COMMON/SETUP/X,Y,N
      COMMON/SIGMAS/SIGMAX,SIGMAX2,SIGMAY,SIGMAX22,AVEX,AVEY
      COMMON/BIGGIES/XLINE,YLINE,GRADIENT,STDDEV,R2
      DOUBLE PRECISION SIGMAX,SIGMAX2,SIGMAY,SIGMAX22,AVEX,AVEY
      DOUBLE PRECISION XLINE,YLINE,GRADIENT,STDDEV,R2
      DOUBLE PRECISION X,Y
      DIMENSION X(100),Y(100)
      CHARACTER*2 YN 
      INTEGER L,N      
      WRITE(*,*)
      WRITE(*,'('' LINEAR REGRESSION PROGRAM FOR THOSE WHO HATE '',
     1/ ,'' EXCEL OR ANYTHING PC!'')')
      WRITE(*,*)
      WRITE(*,'('' BY PAUL F. JOHNSON, JANUARY 1999'')')
      WRITE(*,*)
      WRITE(*,'('' VERSION 1.05. WRITTEN ON AN ACORN RISC PC700'')')
      WRITE(*,*)
      WRITE(*,'('' THIS PROGRAM WILL CALCULATE THE LINE OF BEST'',
     1/ ,'' FIT FOR UPTO 100 DATA POINTS.'')')
      WRITE(*,*)
      WRITE(*,'('' LETS BEGIN!'')')
      WRITE(*,*)
      DO 1 L=1,100
      WRITE(*,'('' ENTER X,Y (OR -999,0 TO START) : '',$)')
      READ *,X(L),Y(L)
      IF (X(L).EQ.-999.) GOTO 2
    1 CONTINUE
    2 N=L-1
      WRITE(*,*)
      WRITE(*,'('' OKAY, I WILL DO THE MATHS NOW'')')
      CALL CALCSIGMAS()
      CALL CALCOTHERS()
      WRITE(*,*)
      WRITE(*,'('' ANSWERS : '')')
      WRITE(*,*)
      PRINT *,'INTERCEPT ON X ',XLINE
      PRINT *,'INTERCEPT ON Y ',YLINE
      PRINT *,'GRADIENT       ',GRADIENT
      PRINT *,'STANDARD DEVN  ',STDDEV
      PRINT *,'R SQUARED      ',R2
      WRITE(*,*)
      WRITE(*,'('' DO YOU WANT A PRINT OUT? (Y/N) : '',$)')
      READ (*,'(A)') YN
      IF (YN.EQ.'Y'.OR.YN.EQ.'y') THEN
      OPEN (55,FORM='PRINTER',ERR=999)
      PRINT *,'X DATA','Y DATA'
      DO 3 L=1,N
      PRINT *,X(L),Y(L)
    3 CONTINUE
      WRITE (55,*)
      PRINT *,'INTERCEPT ON X ',XLINE
      PRINT *,'INTERCEPT ON Y ',YLINE
      PRINT *,'GRADIENT       ',GRADIENT
      PRINT *,'STANDARD DEVN  ',STDDEV
      PRINT *,'R SQUARED      ',R2
      CLOSE(55)
  999 WRITE(*,'('' ERROR WRITING TO PRINTER'')')
      ENDIF            
      WRITE(*,*)
      WRITE(*,'('' AND THATS WHY BANANAS HAVE SKINS ON!'')')
      END

C SUBROUTINE CALCULATE SIGMAS
C CALCULATES ALL THE BITS REQUIRED FOR THE MAIN PROGRAM
C IT'S A BIT EASY THIS ONE!

      SUBROUTINE CALCSIGMAS()
      COMMON/SETUP/X,Y,N
      COMMON/SIGMAS/SIGMAX,SIGMAX2,SIGMAY,SIGMAX22,SIGMAXY,AVEX,AVEY
      INTEGER LL,N
      DOUBLE PRECISION X(100),Y(100)
      DOUBLE PRECISION SIGMAX,SIGMAX2,SIGMAY,SIGMAX22,AVEX,AVEY,SIGMAXY
      DOUBLE PRECISION SX,SY,SXY,SX2
C SIGMAX22 IS THE SAME AS (SIGMA X)^2
      SX=0.D0
      SY=0.D0
      SXY=0.D0
      SX2=0.D0
      DO 1 LL=1,N
      SIGMAX=SX+X(LL)
      SIGMAY=SY+Y(LL)
      SIGMAXY=SXY+(X(LL)*Y(LL))
      SIGMAX2=SX2+(X(LL)**2)
      SX=SIGMAX
      SY=SIGMAY
      SXY=SIGMAXY
      SX2=SIGMAX2
    1 CONTINUE
      SIGMAX22=SIGMAX**2
      AVEX=SIGMAX/N
      AVEY=SIGMAY/N
      RETURN
      END

C SUBROUTINE CALCULATE THE OTHERS
C CALCULATES ALL THE USEFUL BITS SUCH AS THE GRADIENTS ETC
C THIS WILL BE THE FUN BIT!

      SUBROUTINE CALCOTHERS()
      COMMON/SETUP/X,Y,N
      COMMON/SIGMAS/SIGMAX,SIGMAX2,SIGMAY,SIGMAX22,SIGMAXY,AVEX,AVEY
      COMMON/BIGGIES/XLINE,YLINE,GRADIENT,STDDEV,R2
      INTEGER LL,N
      DOUBLE PRECISION XLINE,YLINE,GRADIENT,STDDEV,R2
      DOUBLE PRECISION X(100),Y(100)
      DOUBLE PRECISION SIGMAX,SIGMAX2,SIGMAY,SIGMAX22,AVEX,AVEY,SIGMAXY
      DOUBLE PRECISION M,C,S,R,YEYC2,YEYC
      DOUBLE PRECISION XX,SIGMAXXBYYB,RDBIT1,RDBIT2
      DIMENSION YEYC(100),XXBYYB(100),XXB(100),YYB(100)
      DOUBLE PRECISION XXBYYB,XXB,YYB,YCALC
C OKAY, LET'S DO THE GRADIENT
      M=(SIGMAX*SIGMAY)-(N*SIGMAXY)
      GRADIENT=M/(SIGMAX22-(N*SIGMAX2))
C INTERCEPT ON Y
      C=(SIGMAX*SIGMAXY)-(SIGMAY*SIGMAX2)
      YLINE=C/(SIGMAX22-(N*SIGMAX2))
C INTERCEPT ON X
      XLINE=-YLINE/GRADIENT
C THESE ARE THE BRUTES!
C STANDARD DEVIATION
      M=GRADIENT
      C=YLINE
      YEYC2=0
      XX=XLINE
      DO 1 LL=1,N
      YEYC(LL)=Y(LL)-(M*X(LL)+C)
      YEYC(LL)=YEYC(LL)**2
    1 CONTINUE
      DO 2 LL=1,N
      YEYC2=YEYC2+YEYC(LL)
    2 CONTINUE
      S=YEYC2/(N-2.)
      STDDEV=DSQRT(S)
C R^2
      DO 3 LL=1,N
      XXB(LL)=X(LL)-AVEX
      YYB(LL)=Y(LL)-AVEY
      XXBYYB(LL)=XXB(LL)*YYB(LL)
    3 CONTINUE
      SIGMAXXBYYB=0
      DO 4 LL=1,N
      SIGMAXXBYYB=SIGMAXXBYYB+XXBYYB(LL)
    4 CONTINUE
      RDBIT1=0.D0
      RDBIT2=0.D0
      DO 5 LL=1,N
      RDBIT1=RDBIT1+XXB(LL)**2
      RDBIT2=RDBIT2+YYB(LL)**2
    5 CONTINUE
      R2=SIGMAXXBYYB/DSQRT(RDBIT1*RDBIT2)
      IF(R2.LT.0.) R2=R2*(-1)
      RETURN
      END