Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
Fortran 95 source code to compute multiple regression equations.
!****************************************************
! 重回帰式計算(2元限定)
!
! date name version
! 2018.12.18 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
! 重回帰式計算
!
! :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) :: size_x, size_y, i
real(DP) :: mtx(3, 3), mtx_2(2, 3)
size_x = size(x) / 2
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
mtx(1, 1) = sum_p(y, y)
mtx(2, 1) = sum_p(y, x(:, 1))
mtx(3, 1) = sum_p(y, x(:, 2))
do i = 1, 2
mtx(1, i + 1) = sum_p(x(:, i), y)
mtx(2, i + 1) = sum_p(x(:, i), x(:, 1))
mtx(3, i + 1) = sum_p(x(:, i), x(:, 2))
end do
mtx_2 = transpose(reshape((/ &
& mtx(2, 2), mtx(3, 2), mtx(1, 2), &
& mtx(2, 3), mtx(3, 3), mtx(1, 3) &
& /), (/3, 2/)))
call gauss_e(2, mtx_2)
v = (/mtx_2(1, 3), mtx_2(2, 3)/)
c = calc_const(reshape((/y(:), x(:, 1), x(:, 2)/), (/size_x, 3/)), v)
end subroutine calc_reg_multi
! Sum-of-producsts computation
!
! :param(in) real(8) x(:): 実数配列
! :param(in) real(8) y(:): 実数配列
! :retrun real(8) sum_p: sum of products
real(DP) function sum_p(x, y)
implicit none
real(DP), intent(in) :: x(:), y(:)
integer(SP) :: size_x, size_y, i
real(DP) :: avg_x, avg_y
size_x = size(x)
size_y = size(y)
avg_x = sum(x) / size_x
avg_y = sum(y) / size_y
sum_p = sum((/((x(i) - avg_x) * (y(i) - avg_y), i=1,size_x)/))
end function sum_p
! 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
! Constant term computation
!
! :param(in) real(8) a(:, 3): 配列(目的変数, 説明変数1, 説明変数2)
! :param(in) real(8) v(2): 係数
! :return real(8) c: 定数項
function calc_const(a, v) result(c)
implicit none
real(DP), intent(in) :: a(:, :), v(:)
real(DP) :: c
integer(SP) :: size_a, i
real(DP) :: s(3)
size_a = size(a(:,1))
s = (/(sum(a(:, i)), i=1,3)/)
c = s(1) / size_a
do i = 2, size_a
c = c - s(i) * v(i - 1) / size_a
end do
end function calc_const
end module comp
program regression_multi
use const
use comp
implicit none
real(DP) :: x(10, 2), y(10), c, v(2)
integer(SP) :: i
x = reshape((/ &
& 17.5_DP, 17.0_DP, 18.5_DP, 16.0_DP, 19.0_DP, &
& 19.5_DP, 16.0_DP, 18.0_DP, 19.0_DP, 19.5_DP, &
& 30.0_DP, 25.0_DP, 20.0_DP, 30.0_DP, 45.0_DP, &
& 35.0_DP, 25.0_DP, 35.0_DP, 35.0_DP, 40.0_DP &
& /), (/10, 2/))
y = (/ &
& 45.0_DP, 38.0_DP, 41.0_DP, 34.0_DP, 59.0_DP, &
& 47.0_DP, 35.0_DP, 43.0_DP, 54.0_DP, 52.0_DP &
& /)
do i = 1, 2
print '(A, I0, A, 10F8.2, A)', "説明変数 X(", i, ") = (", x(:, i), ")"
end do
print '(A, 10F8.2, A)', "目的変数 Y = (", y, ")"
print '(A)', "---"
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)
end program regression_multi
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.