Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
Fortran 95 source code to compute a simple regression curve(2d).
!****************************************************
! 単回帰曲線(2次回帰)計算
! : y = a + b * x + c * x^2
!
! date name version
! 2019.03.17 mk-mode.com 1.00 新規作成
!
! Copyright(C) 2019 mk-mode.com All Rights Reserved.
!****************************************************
!
module const
! SP: 単精度(4), DP: 倍精度(8)
integer, parameter :: SP = kind(1.0)
integer(SP), parameter :: DP = selected_real_kind(2 * precision(1.0_SP))
end module const
module comp
use const
implicit none
private
public :: calc_reg_line
contains
! 単回帰曲線(2次回帰)計算
!
! :param(in) real(8) x(:): 説明変数配列
! :param(in) real(8) y(:): 目的変数配列
! :param(out) real(8) a: 係数 a
! :param(out) real(8) b: 係数 b
! :param(out) real(8) b: 係数 c
subroutine calc_reg_line(x, y, a, b, c)
implicit none
real(DP), intent(in) :: x(:), y(:)
real(DP), intent(out) :: a, b, c
integer(SP) :: size_x, size_y, i
real(DP) :: m_x, m_x2, m_x3, m_x4, m_xy, m_x2y, m_y
real(DP) :: s_xx, s_xy, s_xx2, s_x2x2, s_x2y
size_x = size(x)
size_y = size(y)
if (size_x == 0 .or. size_y == 0) then
print *, "[ERROR] array size == 0"
stop
end if
if (size_x /= size_y) then
print *, "[ERROR] size(X) != size(Y)"
stop
end if
m_x = sum(x) / size_x ! avg(X)
m_y = sum(y) / size_y ! avg(Y)
m_x2 = sum(x * x) / size_x ! avg(X^2)
m_x3 = sum(x * x * x) / size_x ! avg(X^3)
m_x4 = sum(x * x * x * x) / size_x ! avg(X^4)
m_xy = sum(x * y) / size_x ! avg(X * Y)
m_x2y = sum(x * x * y) / size_x ! avg(X-2 * Y)
s_xx = m_x2 - m_x * m_x ! Sxx
s_xy = m_xy - m_x * m_y ! Sxy
s_xx2 = m_x3 - m_x * m_x2 ! Sxx2
s_x2x2 = m_x4 - m_x2 * m_x2 ! Sx2x2
s_x2y = m_x2y - m_x2 * m_y ! Sx2y
b = s_xy * s_x2x2 - s_x2y * s_xx2
b = b / (s_xx * s_x2x2 - s_xx2 * s_xx2)
c = s_x2y * s_xx - s_xy * s_xx2
c = c / (s_xx * s_x2x2 - s_xx2 * s_xx2)
a = m_y - b * m_x - c * m_x2
end subroutine calc_reg_line
end module comp
program regression_curve
use const
use comp
implicit none
character(9), parameter :: F_INP = "input.txt"
integer(SP), parameter :: UID = 10
real(DP) :: a, b, c
integer(SP) :: n, i
character(20) :: f
real(DP), allocatable :: x(:), y(:)
! IN ファイル OPEN
open (UID, file = F_INP, status = "old")
! データ数読み込み
read (UID, *) n
! 配列用メモリ確保
allocate(x(n))
allocate(y(n))
! データ読み込み
do i = 1, n
read (UID, *) x(i), y(i)
end do
write (f, '("(A, ", I0, "F8.2, A)")') n
print f, "説明変数 X = (", x, ")"
print f, "目的変数 Y = (", y, ")"
print '(A)', "---"
! IN ファイル CLOSE
close (UID)
! 回帰曲線計算
call calc_reg_line(x, y, a, b, c)
print '(A, F12.8)', "a = ", a
print '(A, F12.8)', "b = ", b
print '(A, F12.8)', "c = ", c
! 配列用メモリ解放
deallocate(x)
deallocate(y)
end program regression_curve
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.