Skip to content

Instantly share code, notes, and snippets.

@komasaru
Last active April 8, 2019 06:26
Show Gist options
  • Save komasaru/456945cf413eecbda4749360dd3c329d to your computer and use it in GitHub Desktop.
Save komasaru/456945cf413eecbda4749360dd3c329d to your computer and use it in GitHub Desktop.
Fortran source code to calculate a coefficent of determination for simple 2D regression.
!****************************************************
! 単回帰分析(2次曲線回帰)の決定係数計算!
!
! date name version
! 2019.03.31 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_curve
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_curve(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_curve
end module comp
program coefficient_of_determination_2d
use const
use comp
implicit none
character(9), parameter :: F_INP = "input.txt"
integer(SP), parameter :: UID = 10
real(DP) :: a, b, c
real(DP) :: y_b, s_r, s_y2, s_e
integer(SP) :: n, i
character(20) :: f
real(DP), allocatable :: x(:), y(:), y_e(:)
! IN ファイル OPEN
open (UID, file = F_INP, status = "old")
! データ数読み込み
read (UID, *) n
! 配列メモリ確保
allocate(x(n))
allocate(y(n))
allocate(y_e(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_curve(x, y, a, b, c)
print '(A, F12.8)', " a = ", a
print '(A, F12.8)', " b = ", b
print '(A, F12.8)', " c = ", c
! 決定係数
y_e = a + b * x + c * x * x ! 推定値
y_b = sum(y) / size(y) ! 標本値 Y (目的変数)の平均
s_r = sum((y_e - y_b)**2) ! 推定値の変動
s_y2 = sum((y - y_b)**2) ! 標本値の変動
s_e = sum((y - y_e)**2) ! 残差の変動
print '(A)', "決定係数"
! 解法-1. 決定係数 (= 推定値の変動 / 標本値の変動)
print '(A, F20.16)', " R2 (1) = ", s_r / s_y2
! 解法-2. 決定係数 (= 1 - 残差の変動 / 標本値の変動)
print '(A, F20.16)', " R2 (2) = ", 1.0 - s_e / s_y2
! 配列メモリ解放
deallocate(x)
deallocate(y)
deallocate(y_e)
end program coefficient_of_determination_2d
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment