daily pastebin goal
6%
SHARE
TWEET

Untitled

a guest Mar 22nd, 2019 70 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. !****************************************************
  2. ! 単回帰曲線(2次回帰)計算
  3. ! : y = a + b * x + c * x^2
  4. ! : 連立方程式を ガウスの消去法で解く方法
  5.  
  6. !   date          name            version
  7. !   2019.03.17    mk-mode.com     1.00 新規作成
  8. !
  9. ! Copyright(C) 2019 mk-mode.com All Rights Reserved.
  10. !****************************************************
  11. !
  12. module const
  13.   ! SP: 単精度(4), DP: 倍精度(8)
  14.   integer,     parameter :: SP = kind(1.0)
  15.   integer(SP), parameter :: DP = selected_real_kind(2 * precision(1.0_SP))
  16. end module const
  17.  
  18. module comp
  19.   use const
  20.   implicit none
  21.   private
  22.   public :: calc_reg_curve
  23.  
  24. contains
  25.   ! 単回帰曲線(2次回帰)計算
  26.   !
  27.   ! :param(in)  real(8) x(:): 説明変数配列
  28.   ! :param(in)  real(8) y(:): 目的変数配列
  29.   ! :param(out) real(8)    a: 係数 a
  30.   ! :param(out) real(8)    b: 係数 b
  31.   ! :param(out) real(8)    b: 係数 c
  32.   subroutine calc_reg_curve(x, y, a, b, c)
  33.     implicit none
  34.     real(DP), intent(in)  :: x(:), y(:)
  35.     real(DP), intent(out) :: a, b, c
  36.     integer(SP) :: size_x, size_y, i
  37.     real(DP)    :: sum_x, sum_x2, sum_x3, sum_x4
  38.     real(DP)    :: sum_y, sum_xy, sum_x2y
  39.     real(DP)    :: mtx(3, 4)
  40.  
  41.     size_x = size(x)
  42.     size_y = size(y)
  43.     if (size_x == 0 .or. size_y == 0) then
  44.       print *, "[ERROR] array size == 0"
  45.       stop
  46.     end if
  47.     if (size_x /= size_y) then
  48.       print *, "[ERROR] size(X) != size(Y)"
  49.       stop
  50.     end if
  51.  
  52.     sum_x   = sum(x)
  53.     sum_x2  = sum(x * x)
  54.     sum_x3  = sum(x * x * x)
  55.     sum_x4  = sum(x * x * x * x)
  56.     sum_y   = sum(y)
  57.     sum_xy  = sum(x * y)
  58.     sum_x2y = sum(x * x * y)
  59.     mtx(1, :) = (/real(size_x, DP),  sum_x, sum_x2,   sum_y/)
  60.     mtx(2, :) = (/           sum_x, sum_x2, sum_x3,  sum_xy/)
  61.     mtx(3, :) = (/          sum_x2, sum_x3, sum_x4, sum_x2y/)
  62.     call solve_ge(3, mtx)
  63.     a = mtx(1, 4)
  64.     b = mtx(2, 4)
  65.     c = mtx(3, 4)
  66.   end subroutine calc_reg_curve
  67.  
  68.   ! 連立方程式を解く(ガウスの消去法)
  69.   !
  70.   ! :param(in)    integer(4)     n: 元数
  71.   ! :param(inout) real(8) a(n,n+1): 係数配列
  72.   subroutine solve_ge(n, a)
  73.     implicit none
  74.     integer(SP), intent(in)    :: n
  75.     real(DP),    intent(inout) :: a(n, n + 1)
  76.     integer(SP) :: i, j
  77.     real(DP)    :: d
  78.  
  79.     ! 前進消去
  80.     do j = 1, n - 1
  81.       do i = j + 1, n
  82.         d = a(i, j) / a(j, j)
  83.         a(i, j+1:n+1) = a(i, j+1:n+1) - a(j, j+1:n+1) * d
  84.       end do
  85.     end do
  86.  
  87.     ! 後退代入
  88.     do i = n, 1, -1
  89.       d = a(i, n + 1)
  90.       do j = i + 1, n
  91.         d = d - a(i, j) * a(j, n + 1)
  92.       end do
  93.       a(i, n + 1) = d / a(i, i)
  94.     end do
  95.   end subroutine solve_ge
  96. end module comp
  97.  
  98. program regression_line
  99.   use const
  100.   use comp
  101.   implicit none
  102.   real(DP)      :: a, b, c
  103.   integer(SP)   :: n, i
  104.   character(20) :: f
  105.   real(DP), allocatable :: x(:), y(:)
  106.  
  107.   read *, n
  108.   allocate(x(n))
  109.   allocate(y(n))
  110.   do i = 1, n
  111.     read *, x(i), y(i)
  112.   end do
  113.   write (f, '("(A, ", I0, "F8.2, A)")') n
  114.   print f, "説明変数 X = (", x, ")"
  115.   print f, "目的変数 Y = (", y, ")"
  116.   print '(A)', "---"
  117.   call calc_reg_curve(x, y, a, b, c)
  118.   print '(A, F12.8)', "a = ", a
  119.   print '(A, F12.8)', "b = ", b
  120.   print '(A, F12.8)', "c = ", c
  121.   deallocate(x)
  122.   deallocate(y)
  123. end program regression_line
RAW Paste Data
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy. OK, I Understand
 
Top