Skip to content

Instantly share code, notes, and snippets.

@komasaru
Created May 8, 2019 01:46
Show Gist options
  • Save komasaru/ef3300e0649b2d921cb6c0c2fa9bc474 to your computer and use it in GitHub Desktop.
Save komasaru/ef3300e0649b2d921cb6c0c2fa9bc474 to your computer and use it in GitHub Desktop.
Fortran 95 source code to calculate a simple regression curve.(3d)
!****************************************************
! 単回帰曲線(3次回帰)計算
! : y = a + b * x + c * x^2 + d * x^3
! : 連立方程式を ガウスの消去法で解く方法
! date name version
! 2019.04.15 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_curve_3d
contains
! 単回帰曲線(3次回帰)計算
!
! :param(in) real(8) x(:): 説明変数配列
! :param(in) real(8) y(:): 目的変数配列
! :param(out) real(8) a: 係数 a
! :param(out) real(8) b: 係数 b
! :param(out) real(8) c: 係数 c
! :param(out) real(8) d: 係数 d
subroutine calc_reg_curve_3d(x, y, a, b, c, d)
implicit none
real(DP), intent(in) :: x(:), y(:)
real(DP), intent(out) :: a, b, c, d
integer(SP) :: size_x, size_y, i
real(DP) :: sum_x, sum_x2, sum_x3, sum_x4, sum_x5, sum_x6
real(DP) :: sum_y, sum_xy, sum_x2y, sum_x3y
real(DP) :: mtx(4, 5)
size_x = size(x)
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
sum_x = sum(x)
sum_x2 = sum(x * x)
sum_x3 = sum(x * x * x)
sum_x4 = sum(x * x * x * x)
sum_x5 = sum(x * x * x * x * x)
sum_x6 = sum(x * x * x * x * x * x)
sum_y = sum(y)
sum_xy = sum(x * y)
sum_x2y = sum(x * x * y)
sum_x3y = sum(x * x * x * y)
mtx(1, :) = (/real(size_x, DP), sum_x, sum_x2, sum_x3, sum_y/)
mtx(2, :) = (/ sum_x, sum_x2, sum_x3, sum_x4, sum_xy/)
mtx(3, :) = (/ sum_x2, sum_x3, sum_x4, sum_x5, sum_x2y/)
mtx(4, :) = (/ sum_x3, sum_x4, sum_x5, sum_x6, sum_x3y/)
call solve_ge(4, mtx)
a = mtx(1, 5)
b = mtx(2, 5)
c = mtx(3, 5)
d = mtx(4, 5)
end subroutine calc_reg_curve_3d
! 連立方程式を解く(ガウスの消去法)
!
! :param(in) integer(4) n: 元数
! :param(inout) real(8) a(n,n+1): 係数配列
subroutine solve_ge(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 solve_ge
end module comp
program regression_curve_3d
use const
use comp
implicit none
character(9), parameter :: F_INP = "input.txt"
integer(SP), parameter :: UID = 10
real(DP) :: a, b, c, d
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))
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 = (", x, ")"
print f, "目的変数 Y = (", y, ")"
print '(A)', "---"
! IN ファイル CLOSE
close (UID)
! 回帰曲線計算
call calc_reg_curve_3d(x, y, a, b, c, d)
print '(A, F12.8)', "a = ", a
print '(A, F12.8)', "b = ", b
print '(A, F12.8)', "c = ", c
print '(A, F12.8)', "d = ", d
! 配列用メモリ解放
deallocate(x)
deallocate(y)
end program regression_curve_3d
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment