相関係数の検定
相関係数の検定（岡本安晴「データ分析のための統計学入門」2009、おうふう、pp．82−84）を行うプログラムt_test_corPrg.pyを以下のように作成した。
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝ t_test_corPrg.py ＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
from Ccalc_t_df_p
import Ccalc
pairs = [[70, 168],
[68, 165], [62, 159], [56, 156], [70, 170],
[64, 162], [59, 162], [70, 169], [68, 168], [69, 170],
[67, 168], [63, 166], [58, 163], [74, 170], [65, 168],
[60, 160], [65, 165], [72, 175], [61, 164], [77, 174],
[68, 164], [64, 161], [70, 166], [47, 147], [66, 165],
[68, 170], [73, 174], [74, 172], [63, 164], [83, 182],
[58, 159], [64, 166], [67, 167], [75, 174], [67, 163],
[68, 166], [68, 167], [66, 165], [63, 165], [55, 154],
[61, 163], [60, 160], [62, 164], [62, 160], [59, 160],
[51, 155], [59, 163], [68, 170], [56, 157], [79, 180]]
calc_p = Ccalc(pairs)
t, df, p =
calc_p.calc()
print("\nt =
%f df =
%d"%(t, df))
print("Prob(x
< %f) = %f Prob(x
> %f) = %f df
= %d"%(t, p, t, 1.0 - p, df))
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
データ対
をリスト
として表し、これらのリストを要素とするリストpairsとして入力データを表している。データを表すリストpairsを実引数としてクラス型Ccalcのオブジェクトcalc_pを作成して、メソッドcalc_p.calc()によりｔ値などを算出している。クラス型Ccalcは以下のように宣言されている。
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝ Ccalc_t_df_p.py ＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
import math
from calc_p_from_t
import p_value
def sqr(x):
return x * x
class Ccalc:
def __init__(self, pairs):
self.pairs = pairs
def calc(self):
n = len(self.pairs)
print("n = %d"%n)
sum_x = 0.0
sum_y = 0.0
for
i in range(0, n):
sum_x += self.pairs[i][0]
sum_y += self.pairs[i][1]
mean_x = sum_x / n
mean_y = sum_y / n
print("Mean_x
= %f"%mean_x)
print("Mean_y = %f"%mean_y)
ssum_x = 0.0
ssum_y = 0.0
sum_xy = 0.0
for
i in range(0, n):
ssum_x += sqr(self.pairs[i][0] - mean_x)
ssum_y += sqr(self.pairs[i][1] - mean_y)
sum_xy += (self.pairs[i][0] - mean_x) * (self.pairs[i][1] - mean_y)
var_x = ssum_x / n
var_y = ssum_y / n
print("Var_x = %f"%var_x)
print("Var_y = %f"%var_y)
sd_x
= math.sqrt(var_x)
sd_y
= math.sqrt(var_y)
print("sd_x = %f"%sd_x)
print("sd_y = %f"%sd_y)
cor
= sum_xy / math.sqrt(ssum_x * ssum_y)
print("r = %f"%cor)
t =
cor * math.sqrt(n - 2.0) / math.sqrt(1.0 - (cor ** 2.0))
df =
n - 2
p = p_value(t, df)
return t, df, p
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
上のプログラムを実行すると図１のような結果を得る。
図１
両側検定でのｐ値は
である。
上の例とは異なるデータを次のように設定してみる。
＝＝＝＝＝＝＝＝＝＝＝＝＝ t_test_corPrg_a.py ＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
from Ccalc_t_df_p
import Ccalc
pairs = [[5, 12],
[7, 15],
[8, 21] ]
calc_p = Ccalc(pairs)
t, df, p =
calc_p.calc()
print("\nt =
%f df =
%d"%(t, df))
print("Prob(x
< %f) = %f Prob(x
> %f) = %f df
= %d"%(t, p, t, 1.0 - p, df))
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
上のプログラムを実行すると図２の結果を得る。
図２
図２の場合、両側検定でのｐ値は
である。
プログラムは圧縮ファイルt_test_corPrg.zipにまとめた。ファイル名t_test_corPrg.zipをクリックしてダウンロードすることができる。