!--------------------------------------------------------------------- ! ! COPYRIGHT NOTICE, DISCLAIMER, and LICENSE: ! ! This notice applies *only* to this specific expression of any ! algorithm, and does not imply ownership or invention of the ! implemented algorithm. ! ! If you modify this file, you may insert additional notices ! immediately following this sentence. ! ! (c) Copyright 2009 - 2013 Peter Charles Crouch. ! All rights reserved, except as noted herein. ! ! This computer program source file is supplied "AS IS". Peter Charles ! Crouch (hereinafter referred to as "Author") disclaims all warranties, ! expressed or implied, including, without limitation, the warranties ! of merchantability and of fitness for any purpose. The Author ! assumes no liability for direct, indirect, incidental, special, ! exemplary, or consequential damages, which may result from the use ! of this software, even if advised of the possibility of such damage. ! ! The Author hereby grants anyone permission to use, copy, modify, and ! distribute this source code, or portions hereof, for any purpose, ! without fee, subject to the following restrictions: ! ! 1. The origin of this source code must not be misrepresented. ! ! 2. Altered versions must be plainly marked as such and must not ! be misrepresented as being the original source. ! ! 3. This Copyright notice may not be removed or altered from any ! source or altered source distribution. ! ! The Author specifically permits (without fee) and encourages the use ! of this source code for entertainment, education, or decoration. If ! you use this source code in a product, acknowledgment is not required ! but would be appreciated. ! ! Acknowledgement: ! This license is based on the wonderful simple license that ! accompanies libpng and that used by Scott Robert Ladd. ! !----------------------------------------------------------------------- ! ! For more information or comments on this software, please contact Peter ! at pccrouch@bcs.org.uk ! !----------------------------------------------------------------------- ! module data_store ! Declaration of arrays used to store data for linear regression implicit none ! Declare real data type with 15 digits of precision and ! range of 10-307 to 10+307 integer, parameter, public :: double = selected_real_kind(15, 307) real(kind=double), dimension(:), allocatable, public :: xvals, yvals end module data_store module file_read private public :: read_file contains subroutine read_file(file_name, nval) ! Open data file and read in number of observations and x and y data ! Should really use iostat=ios and check that ios is zero after each read from file use data_store implicit none character(len=*), intent(in) :: file_name integer, intent(out) :: nval integer, parameter :: in_unit = 10 integer :: astat, ios nval = 0 open(unit=in_unit, status="old", action="read", file=file_name, & position="rewind", iostat=ios) if (ios == 0) then read(unit=in_unit, fmt=*) nval ! Now allocate space for x and y data, program will end if allocation fails allocate(xvals(nval), yvals(nval), stat=astat) if (astat /= 0) then print "(/,a)", "Allocating space for data failed" nval = 0 else read(unit=in_unit, fmt=*) xvals, yvals close(unit=in_unit) end if end if end subroutine read_file end module file_read module reg_calc private public :: calc_reg contains subroutine calc_reg(nval) ! Calculate linear regression for yvals upon xyals ! i.e. y = A + Bx where A is the intercept on the Y axis and B is the slope of ! the best fit straight line use data_store implicit none integer, intent(in) :: nval integer :: i, dastat real(kind=double) :: sumxy, sumxsq, sumysq, ssdureg, ssabreg, intercept, & slope, xbar, ybar, percent, meansq, fvalue, yest character(len=11) :: flabel ! First calculate means for x and y xbar = sum(xvals) / nval ybar = sum(yvals) / nval ! Replace original data with its deviation from means xvals = xvals - xbar yvals = yvals - ybar ! Calculate the corrected sums of squares and products sumxy = 0.0_double sumxsq = 0.0_double sumysq = 0.0_double do i = 1, nval sumxsq = sumxsq + xvals(i) * xvals(i) sumysq = sumysq + yvals(i) * yvals(i) sumxy = sumxy + xvals(i) * yvals(i) end do ! Now calculate regression parameters slope = sumxy / sumxsq intercept = ybar - slope * xbar ssdureg = (sumxy * sumxy) / sumxsq ssabreg = sumysq - ssdureg percent = (100.0_double * ssdureg) / sumysq meansq = ssabreg / (nval - 2) ! Variance ratio (F value) always calculated with larger estimate in the numerator if (ssdureg > meansq) then fvalue = ssdureg / meansq flabel = " F value" else fvalue = meansq / ssdureg flabel = " F' value" end if print "(/,a,f13.6)", "Intercept ", intercept print "(a,f13.6)", "Slope ", slope print "(a,f8.1,a)", "Percentage fit", percent, "%" print "(/,a,f13.6)", "Mean X ", xbar print "(a,f13.6)", "Mean Y ", ybar print "(/,a)", "ANALYSIS OF VARIANCE FOR REGRESSION" print "(/,a)", "Source of Variation Sum of Squares DoF Mean Square" print "(a,f13.6,i6,f16.6)", "Due to regression ", ssdureg, 1, ssdureg print "(a,f13.6,i6,f16.6,a)", "About regression ", ssabreg, nval - 2, meansq, & " Variance" print "(a,f13.6,i6,f16.6,a)", "Total ", sumysq, nval - 1, fvalue, & flabel ! Add means back to input data before calculating residuals xvals = xvals + xbar yvals = yvals + ybar print "(/,a)", "TABLE OF RESIDUALS" print "(/,a)", "Case No. Y Value Y Estimate Residual" do i = 1, nval yest = intercept + slope * xvals(i) print "(i5,3f15.6)", i, yvals(i), yest, yvals(i) - yest end do deallocate(yvals, xvals, stat=dastat) if (dastat /= 0) then print "(/,a)", "Deallocating space for data failed" end if end subroutine calc_reg end module reg_calc program linear ! Example program to calculate simple linear regression ! Reads input from data file, writes output to screen use file_read use reg_calc implicit none integer :: nval character(len=64) :: file_name ! Get the name of the input file from the command line if (command_argument_count() >= 1) then call get_command_argument(1, file_name) ! Open input file and read data into allocated arrays call read_file(trim(file_name), nval) ! Calculate regression and display results if (nval > 0) then call calc_reg(nval) end if end if end program linear