Created
January 31, 2020 02:15
-
-
Save komasaru/283f9d0413b925cc7f902d249c32cea7 to your computer and use it in GitHub Desktop.
Fortran 95 source code to compute multiple regression equations.(3 explanatory variables)
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
!**************************************************** | |
! 重回帰式計算(説明(独立)変数3個限定) | |
! | |
! date name version | |
! 2019.12.08 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_3e | |
contains | |
! 重回帰式計算 | |
! * 説明変数3個限定 | |
! | |
! :param(in) real(8) x(:, 3): 説明変数配列 | |
! :param(in) real(8) y(:): 目的変数配列 | |
! :param(out) real(8) c: 定数 | |
! :param(out) real(8) v(3): 係数 | |
subroutine calc_reg_multi_3e(x, y, c, v) | |
implicit none | |
real(DP), intent(in) :: x(:, :), y(:) | |
real(DP), intent(out) :: c, v(3) | |
integer(SP) :: s_x1, s_x2, s_x3, s_y | |
real(DP) :: sum_x1, sum_x1x1, sum_x1x2, sum_x1x3 | |
real(DP) :: sum_x2, sum_x2x1, sum_x2x2, sum_x2x3 | |
real(DP) :: sum_x3, sum_x3x1, sum_x3x2, sum_x3x3 | |
real(DP) :: sum_y, sum_x1y, sum_x2y, sum_x3y | |
real(DP) :: mtx(4, 5) | |
s_x1 = size(x(:, 1)) | |
s_x2 = size(x(:, 2)) | |
s_x3 = size(x(:, 3)) | |
s_y = size(y) | |
if (s_x1 == 0 .or. s_x2 == 0 .or. s_x3 == 0 .or. s_y == 0) then | |
print *, "[ERROR] array size == 0" | |
stop | |
end if | |
if (s_x1 /= s_y .or. s_x2 /= s_y .or. s_x3 /= s_y) then | |
print *, "[ERROR] size(X) != size(Y)" | |
stop | |
end if | |
sum_x1 = sum(x(:, 1)) | |
sum_x2 = sum(x(:, 2)) | |
sum_x3 = sum(x(:, 3)) | |
sum_x1x1 = sum(x(:, 1) * x(:, 1)) | |
sum_x1x2 = sum(x(:, 1) * x(:, 2)) | |
sum_x1x3 = sum(x(:, 1) * x(:, 3)) | |
sum_x2x1 = sum_x1x2 | |
sum_x2x2 = sum(x(:, 2) * x(:, 2)) | |
sum_x2x3 = sum(x(:, 2) * x(:, 3)) | |
sum_x3x1 = sum_x1x3 | |
sum_x3x2 = sum_x2x3 | |
sum_x3x3 = sum(x(:, 3) * x(:, 3)) | |
sum_y = sum(y) | |
sum_x1y = sum(x(:, 1) * y) | |
sum_x2y = sum(x(:, 2) * y) | |
sum_x3y = sum(x(:, 3) * y) | |
mtx(1, :) = (/real(s_x1, DP), sum_x1, sum_x2, sum_x3, sum_y/) | |
mtx(2, :) = (/ sum_x1, sum_x1x1, sum_x1x2, sum_x1x3, sum_x1y/) | |
mtx(3, :) = (/ sum_x2, sum_x2x1, sum_x2x2, sum_x2x3, sum_x2y/) | |
mtx(4, :) = (/ sum_x3, sum_x3x1, sum_x3x2, sum_x3x3, sum_x3y/) | |
call gauss_e(4, mtx) | |
c = mtx(1, 5) | |
v = mtx(2:4, 5) | |
end subroutine calc_reg_multi_3e | |
! 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 regression_multi_3e | |
use const | |
use comp | |
implicit none | |
character(9), parameter :: F_INP = "input.txt" | |
integer(SP), parameter :: UID = 10 | |
real(DP) :: c, v(3) | |
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, 3)) | |
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(1) = (", x(:, 1), ")" | |
print f, "説明変数 X(2) = (", x(:, 2), ")" | |
print f, "説明変数 X(3) = (", x(:, 3), ")" | |
print f, "目的変数 Y = (", y, ")" | |
print '(A)', "---" | |
! IN ファイル CLOSE | |
close (UID) | |
call calc_reg_multi_3e(x, y, c, v) | |
print '(A, F14.8)', "定数項 = ", c | |
print '(A, F14.8)', "係数-1 = ", v(1) | |
print '(A, F14.8)', "係数-2 = ", v(2) | |
print '(A, F14.8)', "係数-3 = ", v(3) | |
! 配列用メモリ解放 | |
deallocate(x) | |
deallocate(y) | |
end program regression_multi_3e |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment