Fortran - 2 つの配列から相関係数計算！
Updated:
Fortran 95 で、数値からなる同サイズの配列2つを2つの確率変数とみなして相関係数を計算する方法についての記録です。
0. 前提条件
- LMDE 3 (Linux Mint Debian Edition 3; 64bit) での作業を想定。
- GCC 6.3.0 (GFortran 6.3.0) でのコンパイルを想定。
1. アルゴリズムについて
当ブログ過去記事を参照のこと。
2. ソースコードの作成
File: correlation_coefficient.f95
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
!****************************************************
! 相関係数計算
!
! date name version
! 2018.12.18 mk-mode.com 1.00 新規作成
!
! Copyright(C) 2018 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_r
contains
! 相関係数計算
!
! :param(in) real(8) x(:): X 値配列
! :param(in) real(8) y(:): Y 値配列
! :return real(8) r: 相関係数
function calc_r(x, y) result(r)
implicit none
real(DP), intent(in) :: x(:), y(:)
real(DP) :: r
integer(SP) :: size_x, size_y, i
real(DP) :: mean_x, mean_y, cov, var_x, var_y
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
r = 0.0_DP
mean_x = sum(x) / size_x
mean_y = sum(y) / size_y
cov = sum((x(1:size_x) - mean_x) * (y(1:size_x) - mean_y))
var_x = sum((x(1:size_x) - mean_x) * (x(1:size_x) - mean_x))
var_y = sum((y(1:size_x) - mean_y) * (y(1:size_x) - mean_y))
r = (cov / sqrt(var_x)) / sqrt(var_y)
end function calc_r
end module comp
program correlation_coefficient
use const
use comp
implicit none
real(DP) :: x(10), y(10), r
integer(SP) :: i
x = (/(i, i=1,10)/)
y = (/(i, i=1,10)/)
r = calc_r(x, y)
print '("X = (", 10F6.2, ")")', x
print '("Y = (", 10F6.2, ")")', y
print '("r = ", F11.8)', r
print *
y = (/2, 3, 3, 4, 6, 7, 8, 9, 10, 11/)
r = calc_r(x, y)
print '("X = (", 10F6.2, ")")', x
print '("Y = (", 10F6.2, ")")', y
print '("r = ", F11.8)', r
print *
y = (/15, 13, 12, 12, 10, 10, 8, 7, 4, 3/)
r = calc_r(x, y)
print '("X = (", 10F6.2, ")")', x
print '("Y = (", 10F6.2, ")")', y
print '("r = ", F11.8)', r
end program correlation_coefficient
3. ソースコードのコンパイル
$ gfortran -o correlation_coefficient correlation_coefficient.f95
4. 動作確認
$ ./correlation_coefficient
X = ( 1.00 2.00 3.00 4.00 5.00 6.00 7.00 8.00 9.00 10.00)
Y = ( 1.00 2.00 3.00 4.00 5.00 6.00 7.00 8.00 9.00 10.00)
r = 1.00000000
X = ( 1.00 2.00 3.00 4.00 5.00 6.00 7.00 8.00 9.00 10.00)
Y = ( 2.00 3.00 3.00 4.00 6.00 7.00 8.00 9.00 10.00 11.00)
r = 0.99233730
X = ( 1.00 2.00 3.00 4.00 5.00 6.00 7.00 8.00 9.00 10.00)
Y = ( 15.00 13.00 12.00 12.00 10.00 10.00 8.00 7.00 4.00 3.00)
r = -0.98039069
以上。
Comments