Last active
March 29, 2019 02:08
-
-
Save komasaru/e0e28a15062786d91819afd97137b283 to your computer and use it in GitHub Desktop.
Fortran 95 source code to compute a correlation coefficient.
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
!**************************************************** | |
! 相関係数計算 | |
! | |
! 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 |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment