Fortran 95 source code to calculate a simple regression curve.(power)
!**************************************************** | |
! 単回帰曲線(べき乗回帰)計算 | |
! : y = a * x**b | |
! : 連立方程式は ガウスの消去法で解く | |
! date name version | |
! 2019.04.10 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_pow | |
contains | |
! 単回帰曲線(べき乗回帰)計算 | |
! | |
! :param(in) real(8) x(:): 説明変数配列 | |
! :param(in) real(8) y(:): 目的変数配列 | |
! :param(out) real(8) a: 係数 a | |
! :param(out) real(8) b: 係数 b | |
subroutine calc_reg_curve_pow(x, y, a, b) | |
implicit none | |
real(DP), intent(in) :: x(:), y(:) | |
real(DP), intent(out) :: a, b | |
integer(SP) :: size_x, size_y, i | |
real(DP) :: sum_lx, sum_lx2, sum_ly, sum_lxly | |
real(DP) :: mtx(2, 3) | |
real(DP), allocatable :: lx(:), ly(:) | |
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 | |
allocate(lx(size_x)) | |
allocate(ly(size_y)) | |
lx = log(x) | |
ly = log(y) | |
sum_lx = sum(lx) | |
sum_lx2 = sum(lx * lx) | |
sum_ly = sum(ly) | |
sum_lxly = sum(lx * ly) | |
deallocate(lx) | |
deallocate(ly) | |
mtx(1, :) = (/real(size_x, DP), sum_lx, sum_ly/) | |
mtx(2, :) = (/ sum_lx, sum_lx2, sum_lxly/) | |
call solve_ge(2, mtx) | |
a = exp(mtx(1, 3)) | |
b = mtx(2, 3) | |
end subroutine calc_reg_curve_pow | |
! 連立方程式を解く(ガウスの消去法) | |
! | |
! :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_pow | |
use const | |
use comp | |
implicit none | |
character(9), parameter :: F_INP = "input.txt" | |
integer(SP), parameter :: UID = 10 | |
real(DP) :: a, b | |
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_pow(x, y, a, b) | |
print '(A, F12.8)', "a = ", a | |
print '(A, F12.8)', "b = ", b | |
! 配列用メモリ解放 | |
deallocate(x) | |
deallocate(y) | |
stop | |
end program regression_curve_pow |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment