Created
February 2, 2020 02:07
-
-
Save komasaru/8eb0f66097a0bbcd2dc4c0bad7d034b3 to your computer and use it in GitHub Desktop.
Fortran 95 source code to calculate an adjusted coefficent of determination for multiple regression.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
!**************************************************** | |
! 重回帰分析(説明変数2個)決定係数計算 | |
! | |
! DATE AUTHOR VERSION | |
! 2019.12.18 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_multi | |
contains | |
! 重回帰式計算 | |
! * 説明変数2個限定 | |
! | |
! :param(in) real(8) x(:, 2): 説明変数配列 | |
! :param(in) real(8) y(:): 目的変数配列 | |
! :param(out) real(8) c: 定数 | |
! :param(out) real(8) v(2): 係数 | |
subroutine calc_reg_multi(x, y, c, v) | |
implicit none | |
real(DP), intent(in) :: x(:, :), y(:) | |
real(DP), intent(out) :: c, v(2) | |
integer(SP) :: s_x1, s_x2, s_y | |
real(DP) :: sum_x1, sum_x1x1, sum_x1x2 | |
real(DP) :: sum_x2, sum_x2x1, sum_x2x2 | |
real(DP) :: sum_y, sum_x1y, sum_x2y | |
real(DP) :: mtx(3, 4) | |
s_x1 = size(x(:, 1)) | |
s_x2 = size(x(:, 2)) | |
s_y = size(y) | |
if (s_x1 == 0 .or. s_x2 == 0 .or. s_y == 0) then | |
print *, "[ERROR] array size == 0" | |
stop | |
end if | |
if (s_x1 /= s_y .or. s_x2 /= s_y) then | |
print *, "[ERROR] size(X) != size(Y)" | |
stop | |
end if | |
sum_x1 = sum(x(:, 1)) | |
sum_x2 = sum(x(:, 2)) | |
sum_x1x1 = sum(x(:, 1) * x(:, 1)) | |
sum_x1x2 = sum(x(:, 1) * x(:, 2)) | |
sum_x2x1 = sum_x1x2 | |
sum_x2x2 = sum(x(:, 2) * x(:, 2)) | |
sum_y = sum(y) | |
sum_x1y = sum(x(:, 1) * y) | |
sum_x2y = sum(x(:, 2) * y) | |
mtx(1, :) = (/real(s_x1, DP), sum_x1, sum_x2, sum_y/) | |
mtx(2, :) = (/ sum_x1, sum_x1x1, sum_x1x2, sum_x1y/) | |
mtx(3, :) = (/ sum_x2, sum_x2x1, sum_x2x2, sum_x2y/) | |
call gauss_e(3, mtx) | |
c = mtx(1, 4) | |
v = mtx(2:3, 4) | |
end subroutine calc_reg_multi | |
! Gaussian elimination | |
! | |
! :param(in) integer(4) n: 元数 | |
! :param(inout) real(8) a(n,n+1): 係数配列 | |
subroutine gauss_e(n, a) | |
implicit none | |
integer(SP), intent(in) :: n | |
real(DP), intent(inout) :: a(n, n + 1) | |
integer(SP) :: i, j | |
real(DP) :: d | |
! 前進消去 | |
do j = 1, n - 1 | |
do i = j + 1, n | |
d = a(i, j) / a(j, j) | |
a(i, j+1:n+1) = a(i, j+1:n+1) - a(j, j+1:n+1) * d | |
end do | |
end do | |
! 後退代入 | |
do i = n, 1, -1 | |
d = a(i, n + 1) | |
do j = i + 1, n | |
d = d - a(i, j) * a(j, n + 1) | |
end do | |
a(i, n + 1) = d / a(i, i) | |
end do | |
end subroutine gauss_e | |
end module comp | |
program adjusted_coefficient_of_determination | |
use const | |
use comp | |
implicit none | |
character(9), parameter :: F_INP = "input.txt" | |
integer(SP), parameter :: UID = 10 | |
real(DP) :: c, v(2) | |
real(DP) :: y_b, s_r, s_e, r_2 | |
integer(SP) :: n, p, 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, 2)) | |
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(1) = (", x(:, 1), ")" | |
print f, "説明変数 X(2) = (", x(:, 2), ")" | |
print f, "目的変数 Y = (", y, ")" | |
print '(A)', "---" | |
! IN ファイル CLOSE | |
close (UID) | |
! 重回帰式 | |
call calc_reg_multi(x, y, c, v) | |
print '(A, F12.8)', " 定数項 = ", c | |
print '(A, F12.8)', " 係数-1 = ", v(1) | |
print '(A, F12.8)', " 係数-2 = ", v(2) | |
! 自由度調整済み決定係数 | |
y_e = c + v(1) * x(:, 1) + v(2) * x(:, 2) ! 推定値配列 | |
y_b = sum(y) / size(y) ! 標本値 Y (目的変数)の平均 | |
p = size(x(1, :)) | |
n = size(y) | |
s_r = sum((y - y_b)**2) ! 推定値の変動 | |
s_e = sum((y - y_e)**2) ! 残差の変動 | |
r_2 = 1.0_DP - (s_e / (n - p - 1)) / (s_r / (n - 1)) | |
print '(A)', "自由度調整済み決定係数" | |
print '(A, F20.16)', " R2f = ", r_2 | |
! 配列メモリ解放 | |
deallocate(x) | |
deallocate(y) | |
deallocate(y_e) | |
end program adjusted_coefficient_of_determination |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment