Instantly share code, notes, and snippets.

# komasaru/regression_multi_v2.f95

Created August 14, 2019 05:25
Show Gist options
• Save komasaru/bfba8cdd1531f47a2621fed1d8dcce12 to your computer and use it in GitHub Desktop.
Fortran 95 source code to compute multiple regression equations.(v2)
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 name version ! 2019.06.01 mk-mode.com 1.00 新規作成 ! ! Copyright(C) 2018 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, i 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 regression_multi use const use comp implicit none character(9), parameter :: F_INP = "input.txt" integer(SP), parameter :: UID = 10 real(DP) :: c, v(2) 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, 2)) 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, "目的変数 Y = (", y, ")" print '(A)', "---" ! IN ファイル CLOSE close (UID) call calc_reg_multi(x, y, c, v) print '(A, F14.8)', "定数項 = ", c print '(A, F14.8)', "係数-1 = ", v(1) print '(A, F14.8)', "係数-2 = ", v(2) ! 配列用メモリ解放 deallocate(x) deallocate(y) end program regression_multi