Skip to content

Instantly share code, notes, and snippets.

@komasaru
Created August 15, 2019 02:12
Show Gist options
  • Save komasaru/1a21421eaf0502c096940819295ee1e0 to your computer and use it in GitHub Desktop.
Save komasaru/1a21421eaf0502c096940819295ee1e0 to your computer and use it in GitHub Desktop.
Fortran 95 source code to compute multiple regression equations.(2d)
!****************************************************
! 重回帰式計算(説明(独立)変数2個、2次多項式モデル)
! * y = b0 + b1x1 + b2x2 + b3x1x2 + b4x1^2 + b5x2^2
! * y = b0 + b1x1 + b2x2 + b3x3 + b4x4 + b5x5
! (x3 = x1x2, x4 = x1^2, x5 = x2^2)
! ということ。
! date name version
! 2019.06.27 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_2d
contains
! 重回帰式計算
! * 説明変数2個、2次多項式モデル
! * y = b0 + b1x1 + b2x2 + b3x1x2 + b4x1^2 + b5x2^2
! * y = b0 + b1x1 + b2x2 + b3x3 + b4x4 + b5x5
! (x3 = x1x2, x4 = x1^2, x5 = x2^2)
!
! :param(in) real(8) x(:, 2): 説明変数配列
! :param(in) real(8) y(:): 目的変数配列
! :param(out) real(8) c: 定数
! :param(out) real(8) v(5): 係数
subroutine calc_reg_multi_2d(x, y, c, v)
implicit none
real(DP), intent(in) :: x(:, :), y(:)
real(DP), intent(out) :: c, v(5)
integer(SP) :: s_x1, s_x2, s_y
real(DP) :: mtx(6, 7)
real(DP), allocatable :: x1(:), x2(:), x3(:), x4(:), x5(:)
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
allocate(x1(s_x1))
allocate(x2(s_x1))
allocate(x3(s_x1))
allocate(x4(s_x1))
allocate(x5(s_x1))
x1 = x(:, 1)
x2 = x(:, 2)
x3 = x1 * x2
x4 = x1 * x1
x5 = x2 * x2
! 左辺・対角成分
mtx(1, 1) = s_x1
mtx(2, 2) = sum(x1 * x1)
mtx(3, 3) = sum(x2 * x2)
mtx(4, 4) = sum(x3 * x3)
mtx(5, 5) = sum(x4 * x4)
mtx(6, 6) = sum(x5 * x5)
! 左辺・右上成分
mtx(1, 2) = sum(x1)
mtx(1, 3) = sum(x2)
mtx(1, 4) = sum(x3)
mtx(1, 5) = sum(x4)
mtx(1, 6) = sum(x5)
mtx(2, 3) = sum(x1 * x2)
mtx(2, 4) = sum(x1 * x3)
mtx(2, 5) = sum(x1 * x4)
mtx(2, 6) = sum(x1 * x5)
mtx(3, 4) = sum(x2 * x3)
mtx(3, 5) = sum(x2 * x4)
mtx(3, 6) = sum(x2 * x5)
mtx(4, 5) = sum(x3 * x4)
mtx(4, 6) = sum(x3 * x5)
mtx(5, 6) = sum(x4 * x5)
! 左辺・左下成分
mtx(2, 1) = mtx(1, 2)
mtx(3, 1) = mtx(1, 3)
mtx(3, 2) = mtx(2, 3)
mtx(4, 1) = mtx(1, 4)
mtx(4, 2) = mtx(2, 4)
mtx(4, 3) = mtx(3, 4)
mtx(5, 1) = mtx(1, 5)
mtx(5, 2) = mtx(2, 5)
mtx(5, 3) = mtx(3, 5)
mtx(5, 4) = mtx(4, 5)
mtx(6, 1) = mtx(1, 6)
mtx(6, 2) = mtx(2, 6)
mtx(6, 3) = mtx(3, 6)
mtx(6, 4) = mtx(4, 6)
mtx(6, 5) = mtx(5, 6)
! 右辺
mtx(1, 7) = sum( y)
mtx(2, 7) = sum(x1 * y)
mtx(3, 7) = sum(x2 * y)
mtx(4, 7) = sum(x3 * y)
mtx(5, 7) = sum(x4 * y)
mtx(6, 7) = sum(x5 * y)
deallocate(x1)
deallocate(x2)
deallocate(x3)
deallocate(x4)
deallocate(x5)
call gauss_e(6, mtx)
c = mtx(1, 7)
v = mtx(2:6, 7)
end subroutine calc_reg_multi_2d
! 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_2d
use const
use comp
implicit none
character(9), parameter :: F_INP = "input.txt"
integer(SP), parameter :: UID = 10
real(DP) :: c, v(5)
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_2d(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)
print '(A, F14.8)', "係数-4 = ", v(4)
print '(A, F14.8)', "係数-5 = ", v(5)
! 配列用メモリ解放
deallocate(x)
deallocate(y)
end program regression_multi_2d
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment