• API
• FAQ
• Tools
• Archive
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. !
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.
108.   allocate(x(n))
109.   allocate(y(n))
110.   do i = 1, n
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.

Top