Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
Fortran 95 source code to compute a correlation coefficient.
!****************************************************
! 相関係数計算
!
! 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
You can’t perform that action at this time.