Created
May 8, 2019 02:24
-
-
Save komasaru/f3af965de671e3b60a93e5da8825cd5f to your computer and use it in GitHub Desktop.
Fortran 95 source code to calculate a simple regression curve.(exp)
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
!**************************************************** | |
! 単回帰曲線(指数(ab指数)回帰)計算 | |
! : y = a * b**x | |
! : 連立方程式は ガウスの消去法で解く | |
! 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_exp | |
contains | |
! 単回帰曲線(指数(ab指数)回帰)計算 | |
! | |
! :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_exp(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_x, sum_x2, sum_ly, sum_xly | |
real(DP) :: mtx(2, 3) | |
real(DP), allocatable :: 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(ly(size_y)) | |
ly = log(y) | |
sum_x = sum(x) | |
sum_x2 = sum(x * x) | |
sum_ly = sum(ly) | |
sum_xly = sum(x * ly) | |
deallocate(ly) | |
mtx(1, :) = (/real(size_x, DP), sum_x, sum_ly/) | |
mtx(2, :) = (/ sum_x, sum_x2, sum_xly/) | |
call solve_ge(2, mtx) | |
a = exp(mtx(1, 3)) | |
b = exp(mtx(2, 3)) | |
end subroutine calc_reg_curve_exp | |
! 連立方程式を解く(ガウスの消去法) | |
! | |
! :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_exp | |
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_exp(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_exp |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment