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