Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- !****************************************************
- ! 単回帰曲線(2次回帰)計算
- ! : y = a + b * x + c * x^2
- ! : 連立方程式を ガウスの消去法で解く方法
- ! date name version
- ! 2019.03.17 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
- contains
- ! 単回帰曲線(2次回帰)計算
- !
- ! :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) b: 係数 c
- subroutine calc_reg_curve(x, y, a, b, c)
- implicit none
- real(DP), intent(in) :: x(:), y(:)
- real(DP), intent(out) :: a, b, c
- integer(SP) :: size_x, size_y, i
- real(DP) :: sum_x, sum_x2, sum_x3, sum_x4
- real(DP) :: sum_y, sum_xy, sum_x2y
- real(DP) :: mtx(3, 4)
- 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_y = sum(y)
- sum_xy = sum(x * y)
- sum_x2y = sum(x * x * y)
- mtx(1, :) = (/real(size_x, DP), sum_x, sum_x2, sum_y/)
- mtx(2, :) = (/ sum_x, sum_x2, sum_x3, sum_xy/)
- mtx(3, :) = (/ sum_x2, sum_x3, sum_x4, sum_x2y/)
- call solve_ge(3, mtx)
- a = mtx(1, 4)
- b = mtx(2, 4)
- c = mtx(3, 4)
- end subroutine calc_reg_curve
- ! 連立方程式を解く(ガウスの消去法)
- !
- ! :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_line
- use const
- use comp
- implicit none
- real(DP) :: a, b, c
- integer(SP) :: n, i
- character(20) :: f
- real(DP), allocatable :: x(:), y(:)
- read *, n
- allocate(x(n))
- allocate(y(n))
- do i = 1, n
- read *, x(i), y(i)
- end do
- write (f, '("(A, ", I0, "F8.2, A)")') n
- print f, "説明変数 X = (", x, ")"
- print f, "目的変数 Y = (", y, ")"
- print '(A)', "---"
- call calc_reg_curve(x, y, a, b, c)
- print '(A, F12.8)', "a = ", a
- print '(A, F12.8)', "b = ", b
- print '(A, F12.8)', "c = ", c
- deallocate(x)
- deallocate(y)
- end program regression_line
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement