「入門はじめての統計解析」をPythonで写経 Vol.2 ～ １章「はじめての平均・分散・標準偏差」②散布図、相関係数、順位相関係数、オッズ、オッズ比
１章「はじめての平均・分散・標準偏差」
書籍の著者 石村貞夫 先生
この記事は、書籍「入門はじめての統計解析」１章「はじめての平均・分散・標準偏差」の Python写経活動 を取り扱います。
書籍の図・表・計算を淡々とPython化する写経シリーズです。
後編では２変数の関係を探り、散布図、共分散、相関係数、順位相関係数、オッズ、オッズ比などを実践します。
では書籍を開いて統計解析の旅に出発です🚀
はじめに
このブログシリーズは、書籍「入門はじめての統計解析」（東京図書、「テキスト」と呼びます）の Python 写経を通じて得た「統計解析の楽しさ」をご紹介します。
書籍の紹介と引用表記はリンク先の記事に掲載しています。
１章 はじめての平均・分散・標準偏差
この記事は１章の以下のSectionを取り扱います。
1.3 散布図と相関係数
1.4 いろいろな順位相関係数
1.5 クロス集計表とオッズ比
記事に用いるデータは、テキストに掲載されたデータそのものを引用しています。
データ件数の少ないものはコード上でデータを登録し、データ件数の多いものはCSVファイル化してデータを読み込みしています。
１章で用いるライブラリをインポートします。
### インポート
# 数値計算
import collections # python標準ライブラリ
import numpy as np
import pandas as pd
# 統計
import scipy.stats as stats
import pingouin as pg
# 描画
import matplotlib.pyplot as plt
import seaborn as sns
plt.rcParams['font.family'] = 'Meiryo'Section 1.3 散布図と相関係数
「２つの量的変数」を取り扱うSectionです。
■ 散布図 － ２変量の関係を見る p.36
乗用車の価格と最大出力のデータを登録します。
p.37 表1.3.3「相関表」のデータをお借りします。
### 例 乗用車の価格と最大出力の相関表と散布図 p.37 表1.3.3
# 相関表の登録 x:価格、y:最大出力
data8 = pd.DataFrame({'価格': [140, 101, 149, 82, 118, 98, 143, 69, 137],
'最大出力': [120, 113, 160, 76, 105, 93, 140, 73, 140]},
index=range(1, 10))
# 結果表示
print('data8.shape: ', data8.shape)
data8【実行結果】
乗用車の種類は９つです。
散布図を描画します。
p.37 図1.3.2「散布図」に相当します。
seaborn ライブラリ（可視化ライブラリ）を利用します。
# 散布図の描画 p.37 図1.3.2
ax = sns.scatterplot(data=data8, x='価格', y='最大出力', s=70)
ax.set(xlim=(60, 170), ylim=(60, 170));【実行結果】
右上がりの傾向が見られます。
価格が高くなると最大出力がが大きくなるような関係に見えます。
相関係数を確認します。
pandas データフレームの corr() を利用します。
# 相関係数
data8.corr().round(3)【実行結果】
価格と最大出力の相関係数は $${0.936}$$ であり、強い正の相関関係があると言えそうです。
次のデータに移ります。
妊産婦受診率と新生児死亡率のデータを登録します。
p.37 表1.3.4「相関表」のデータをお借りします。
### 例 妊産婦受診率と新生児死亡率の相関表と散布図 p.37 表1.3.3
# 相関表の登録 x:受診率、y:死亡率
data9 = pd.DataFrame({
'受診率': [1.54, 2.18, 9.59, 5.16, 7.39, 2.08, 4.64, 3.81, 2.38, 9.07, 3.74,
1.28],
'死亡率': [4.26, 5.35, 3.68, 4.72, 3.46, 3.91, 3.85, 5.02, 4.36, 4.15, 5.79,
5.63]},
index=range(1, 13))
# 結果表示
print('data9.shape: ', data9.shape)
data9【実行結果】
１２件のデータです。
受診率と死亡率の単位は％です。
散布図を描画します。
p.37 図1.3.3「散布図」に相当します。
seaborn ライブラリ（可視化ライブラリ）を利用します。
# 散布図の描画 p.37 図1.3.3
ax = sns.scatterplot(data=data9, x='受診率', y='死亡率', s=70)
ax.set(xlim=(0, 11), ylim=(2, 7));【実行結果】
右下がりの傾向が見られます。
受診率が大きくなると死亡率が小さくなるような関係に見えます。
相関係数を確認します。
# 相関係数
data9.corr().round(3)【実行結果】
受診率と死亡率の相関係数は $${-0.534}$$ であり、負の相関関係があると言えそうです。
次のデータに移ります。
ダムの高さと総貯水量のデータを登録します。
p.38 表1.3.5「相関表」のデータをお借りします。
### 例 10カ所のダムの高さと総貯水量の相関表と散布図 p.38 表1.3.5
# 相関表の登録 x:高さ[m]、y:総貯水量[100万㎥]
data10 = pd.DataFrame({
'高さ': [131, 131, 133, 140, 145, 149, 155, 156, 157, 186],
'総貯水量': [204, 370, 46, 223, 497, 189, 123, 327, 601, 199]},
index=range(1, 11))
# 結果表示
print('data10.shape: ', data10.shape)
data10【実行結果】
データの件数は１０です。
散布図を描画します。
p.37 図1.3.4「散布図」に相当します。
seaborn ライブラリ（可視化ライブラリ）を利用します。
# 散布図の描画 p.38 図1.3.4
ax = sns.scatterplot(data=data10, x='高さ', y='総貯水量', s=70)
ax.set(xlim=(120, 200), ylim=(0, 650));【実行結果】
右上がり、右下がりの傾向が見られません。
相関係数を確認します。
# 相関係数
data10.corr().round(3)【実行結果】
ダムの高さと総貯水量の相関係数は $${0.081}$$ であり、ほぼ０です。
相関関係が無い「無相関」と言えそうです。
■ 相関係数 p.40
p.40 の相関係数の定義の数式を用いて、相関係数算出関数を【Python標準縛り】で作成します。
最終行でダムデータの相関係数を計算します。
相関係数の定義の数式はぜひテキストでご確認下さい。
### 相関係数の定義 p.40
# 相関係数算出関数
def calc_corr(list_x, list_y):
## 設定と準備
N = len(list_x) # 標本サイズ len(x)=len(y)
x_bar = sum(list_x) / N # xの標本平均
y_bar = sum(list_y) / N # yの標本平均
## 相関係数の計算に必要な分子、分母の計算
# 分子: xとyの偏差積和
numerator = sum([(x_i - x_bar)*(y_i - y_bar)
for x_i, y_i in zip(list_x, list_y)])
# 分母1: xの偏差平方和、分母2: yの偏差平方和
denominator1 = sum([(x_i - x_bar)**2 for x_i in list_x])**(1/2)
denominator2 = sum([(y_i - y_bar)**2 for y_i in list_y])**(1/2)
## 戻り値：相関係数
return numerator / (denominator1 * denominator2)
# ダムのデータでテスト
calc_corr(data10['高さ'], data10['総貯水量'])【実行結果】
ダムデータの相関係数と合っています。
■ 共分散 － それは２変量の広がり p.41
p.41 の共分散の定義の数式を用いて、共分散算出関数を【Python標準縛り】で作成します。
最終行でダムデータの共分散を計算します。
共分散の定義の数式はぜひテキストでご確認下さい。
### 共分散の定義 p.41
# 共分散算出関数
def calc_covariance(list_x, list_y):
## 設定と準備
N = len(list_x) # 標本サイズ len(x)=len(y)
x_bar = sum(list_x) / N # xの標本平均
y_bar = sum(list_y) / N # yの標本平均
## 共分散の計算に必要な分子: xとyの偏差積和の計算
numerator = sum([(x_i - x_bar)*(y_i - y_bar)
for x_i, y_i in zip(list_x, list_y)])
## 戻り値：共分散
return numerator / (N - 1)
# ダムのデータでテスト
calc_covariance(data10['高さ'], data10['総貯水量'])【実行結果】
pandasの共分散計算で検算します。
pandas データフレームの cov() を利用します。
# pandasで検算
data10.cov().iloc[0, 1]【実行結果】
自作関数の結果と一致しました。
■ 相関係数と $${\cos \theta}$$ の関係 p.42
テキストによると「相関係数＝$${\cos \theta}$$」です。
p.42の以下の式で $${\cos \theta}$$ を求めて、相関係数と一致するか確認してみましょう。
$$
\cos \theta = \cfrac{\boldsymbol{x \cdot y}}{\boldsymbol{\|x\| \cdot \|y\|}}
$$
### ダムのデータで検算
# xの偏差ベクトルとyの偏差ベクトルを作成
vec_x = data10['高さ'].values - data10['高さ'].mean()
vec_y = data10['総貯水量'].values - data10['総貯水量'].mean()
# cosθの算出 テキストp.42の数式を利用
cos_theta = vec_x @ vec_y / (np.linalg.norm(vec_x) * np.linalg.norm(vec_y))
cos_theta【実行結果】
pandas で相関係数を算出して検算します。
# pandasで検算
data10.corr().iloc[0, 1]【実行結果】
ほぼ一致しました。
■ データの標準化 － 単位にとらわれない統計処理？！ p.44
p.44 表1.3.6 「センチメートルとメートル？」のデータを登録します。
### 例 センチメートルとメートル p.44 表1.3.6
# 表の登録
data11 = pd.DataFrame({
'身長(cm)': [158, 154, 162, 160, 153, 155, 163, 157, 155, 148]},
index=range(1, 11))
data11['身長(m)'] = data11['身長(cm)'] / 100
# 結果表示
print('data11.shape: ', data11.shape)
data11.T【実行結果】
p.44の平均値$${\bar{x}}$$を求めます。
pandas の mean() を利用します。
### 平均値 p.44
print(f"平均 cm: {data11['身長(cm)'].mean()}")
print(f"平均 m : {data11['身長(m)'].mean():.3f}")【実行結果】
p.44の分散$${s^2}$$（不偏分散）を求めます。
pandas の var() を利用します。
### 分散 p.44
print(f"分散 cm: {data11['身長(cm)'].var(ddof=1):.2f}")
print(f"分散 m : {data11['身長(m)'].var(ddof=1):.6f}")【実行結果】
p.45 のデータの標準化の定義の数式を用いて、データの標準化関数を【Python標準縛り】で作成します。
### データの標準化の定義
def data_standardization(list):
## 計算要素の算出
# 標本サイズ
N = len(list)
# 標本平均
x_bar = sum(list) / N
# 標本標準偏差
x_std = (sum([(x_i - x_bar)**2 for x_i in list]) / (N - 1))**(1/2)
## 戻り値: 標準化データ
return (list - x_bar) / x_stdp.45 の表1.3.6のデータの標準化を実行します。
こちらは単位がcm。
### 表1.3.6のデータの標準化 p.45
# cm
std_cm = data_standardization(data11['身長(cm)']).values
std_cm【実行結果】
こちらは単位がmです。
# m
std_m = data_standardization(data11['身長(m)']).values
std_m【実行結果】
データの単位 cm、m を標準化したデータは同じ値になりました。
データの標準化によって、平均０、分散１、標準偏差１になることを確認します。
### 標準化されたデータの平均値 p.45
std_cm.mean()【実行結果】
平均は０です。
### 標準化されたデータの分散 p.45
std_cm.var(ddof=1).round(6)【実行結果】
分散は１です。
### 標準化されたデータの標準偏差 p.45
std_cm.std(ddof=1).round(6)【実行結果】
標準偏差は１です。
■ 無相関的関係？！ p.47
p.47 の表1.3.7「２変量の$${x, y}$$の関係」のデータを登録します。
### データの作成 p.47 表1.3.7
# xの値をデータフレーム化
data12 = pd.DataFrame({'x': [2, 5, 0, 4, 1, 6, 3]}, index=range(1, 8))
# y, x², y², xyの各列を追加
data12['y'] = data12['x']**2 - 6*data12['x'] + 10
data12['x²'] = data12['x']**2
data12['y²'] = data12['y']**2
data12['xy'] = data12['x'] * data12['y']
# 結果の表示
print('data12.shape: ', data12.shape)
data12【実行結果】
$${x, y}$$の相関係数を pandas の corr で算出します。
### 相関係数rの算出 pandasのcorrメソッドを利用
data12[['x', 'y']].corr().round(6)【実行結果】
相関係数は０です。$${x,y}$$は無相関です。
散布図を描画します。
### x,yの散布図の描画
# x,yの散布図の描画
sns.scatterplot(data=data12, x='x', y='y', s=70)
# y=x^2-6x+10の曲線の描画
x_val = np.linspace(-1, 10, 101)
plt.plot(x_val, x_val**2 - 6*x_val + 10, color='tab:red', ls='--')
# 修飾
plt.xlim(-0.5, 6.5)
plt.ylim(0, 11);【実行結果】
２次式$${y=x^2-6x+10}$$の関係があります。
無相関だけど関係がある、ということです。
相関係数は直線的でないデータ（非線形なデータ）の相関をうまく表せないようです。
■ 例題：相関係数・共分散の求め方 p.49
テキストの手順に沿って相関係数と共分散を求めます。
手順１：表1.3.9「統計量」のデータの作成します。
### 手順1 データの作成 p.49 表1.3.9
data13 = pd.DataFrame({'x': [158, 154, 162, 160, 153, 155, 163, 157, 155,148],
'y': [62, 45, 47, 46, 40, 41, 55, 53, 45, 48]},
index=range(1, 11))
data13['x²'] = data13.x**2
data13['y²'] = data13.y**2
data13['xy'] = data13.x * data13.y
print('data13.shape: ', data13.shape)
data13【実行結果】
列の合計を算出します。
# 合計チェック
data13.sum(axis=0).rename('合計').to_frame().T【実行結果】
手順２：$${x,y}$$の相関係数と共分散を算出します。
【Python標準縛り】です。
### 手順2 x,yの相関係数rと共分散Cov(x,y)を計算
# 計算要素の算出
N = len(data13) # 標本サイズ
x_bar = data13.x.sum() / N # xの標本平均
y_bar = data13.y.sum() / N # yの標本平均
# x,yの標準偏差（不偏分散の標準偏差）の算出
std_x = (sum((data13.x - x_bar)**2) / (N - 1))**(1/2)
std_y = (sum((data13.y - y_bar)**2) / (N - 1))**(1/2)
# 共分散Cov(x,y), 相関係数rの算出
cov_xy = sum((data13.x - x_bar) * (data13.y - y_bar)) / (N - 1)
r = cov_xy / (std_x * std_y)
# 結果の表示
print(f'相関係数r = {r:.3f}, 共分散Cov(x,y) = {cov_xy:.3f}')
print(f'xの分散Var(x) = {std_x**2:.3f}, yの分散Var(y) = {std_y**2:.3f}')【実行結果】
pandas で検算します。
### pandasで計算
vars = data13[['x', 'y']]
print(f'相関係数r = {vars.corr().iloc[0, 1]:.3f}, '
f'共分散Cov(x,y) = {vars.cov().iloc[0, 1]:.3f}')
print(f'xの分散Var(x) = {vars.x.var(ddof=1):.3f}, '
f'yの分散Var(y) = {vars.y.var(ddof=1):.3f}')【実行結果】
■ 理解度チェック：相関係数・共分散の求め方
【問題１】p.50
乗用車の価格と最大出力の相関係数と共分散を求めます。
表1.3.10「乗用車の価格と最大出力の統計量」を作成します。
### データの作成 p.50 表1.3.10
# x,yをデータフレーム化
data14 = pd.DataFrame(
{'x': [140, 101, 149, 82, 118, 98, 143, 69, 137], # x:価格
'y': [120, 113, 160, 76, 105, 93, 140, 73, 140]}, # y:最大出力
index=range(1, 10)
)
# x²,y²,xy列を追加
data14['x²'] = data14.x**2
data14['y²'] = data14.y**2
data14['xy'] = data14.x * data14.y
# 結果確認
print('data14.shape: ', data14.shape)
data14【実行結果】
合計を算出します。
# 合計チェック
data14.sum(axis=0).rename('合計').to_frame().T【実行結果】
データの散布図を見ましょう。
# 散布図
sns.scatterplot(data=data14, x='x', y='y', s=70)
plt.xlabel('価格')
plt.ylabel('最大出力');【実行結果】
正の相関関係がありそうです。
分散、共分散、相関係数を求めます。
【Python標準縛り】です。
### 価格の分散 p.50
N = len(data14)
x_bar = sum(data14.x) / N
var_x = sum((data14.x - x_bar)**2) / (N - 1)
print(f'価格の分散 = {var_x:.3f}')【実行結果】
### 出力の分散 p.50
y_bar = sum(data14.y) / N
var_y = sum((data14.y - y_bar)**2) / (N - 1)
print(f'出力の分散 = {var_y:.3f}')【実行結果】
### 価格と出力の共分散 p.50
cov_xy = sum((data14.x - x_bar)*(data14.y - y_bar)) / (N - 1)
print(f'価格と出力の共分散 = {cov_xy:.3f}')【実行結果】
### 価格と出力の相関係数 p.50
r = cov_xy / (var_x*var_y)**(1/2)
print(f'価格と出力の相関係数 = {r:.3f}')【実行結果】
pandas で検算します。
### pandasで答え合わせ
print(f'価格の分散 = {data14.x.var(ddof=1):.3f}')
print(f'出力の分散 = {data14.y.var(ddof=1):.3f}')
print(f"価格と出力の共分散 = {data14[['x', 'y']].cov(ddof=1).iloc[0, 1]:.3f}")
print(f"価格と出力の相関係数 = {data14[['x', 'y']].corr().iloc[0, 1]:.3f}")【実行結果】
【問題２】p.51
発展途上国の妊産婦受診率と新生児死亡率の相関係数と共分散を求めます。
表1.3.11「受診率と死亡率の統計量」を作成します。
### データの作成 p.51 表1.3.11
# x,yをデータフレーム化
data15 = pd.DataFrame({
# x: 受診率
'x': [1.54, 2.18, 9.59, 5.16, 7.39, 2.08, 4.64, 3.81, 2.38, 9.07, 3.74,
1.28],
# y: 死亡率
'y': [4.26, 5.35, 3.68, 4.72, 3.46, 3.91, 3.85, 5.02, 4.36, 4.15, 5.79,
5.63]},
index=range(1, 13)
)
# x²,y²,xy列を追加
data15['x²'] = data15.x**2
data15['y²'] = data15.y**2
data15['xy'] = data15.x * data15.y
# 結果確認
print('data15.shape: ', data15.shape)
data15【実行結果】
合計を算出します。
# 合計チェック
data15.sum(axis=0).rename('合計').to_frame().T【実行結果】
散布図を見ましょう。
# 散布図
sns.scatterplot(data=data15, x='x', y='y', s=70)
plt.xlabel('受診率')
plt.ylabel('死亡率');【実行結果】
負の相関関係がありそうです。
分散、共分散、相関係数を求めます。
【Python標準縛り】です。
### 受診率の分散 p.51
N = len(data15)
x_bar = sum(data15.x) / N
var_x = sum((data15.x - x_bar)**2) / (N - 1)
print(f'受診率の分散 = {var_x:.3f}')【実行結果】
### 死亡率の分散 p.51
y_bar = sum(data15.y) / N
var_y = sum((data15.y - y_bar)**2) / (N - 1)
print(f'死亡率の分散 = {var_y:.3f}')【実行結果】
### 受診率と死亡率の共分散 p.51
cov_xy = sum((data15.x - x_bar)*(data15.y - y_bar)) / (N - 1)
print(f'受診率と死亡率の共分散 = {cov_xy:.3f}')【実行結果】
### 受診率と死亡率の相関係数 p.51
r = cov_xy / (var_x*var_y)**(1/2)
print(f'価格と出力の相関係数 = {r:.3f}')【実行結果】
pandas で検算します。
### pandasで答え合わせ
print(f'受診率の分散 = {data15.x.var(ddof=1):.3f}')
print(f'死亡率の分散 = {data15.y.var(ddof=1):.3f}')
print(f"受診率と死亡率の共分散 = {data15[['x', 'y']].cov(ddof=1).iloc[0, 1]:.3f}")
print(f"受診率と死亡率の相関係数 = {data15[['x', 'y']].corr().iloc[0, 1]:.3f}")【実行結果】
Section 1.4 いろいろな順位相関係数
テキストは「２つの変量 $${x, y}$$ の測定値が数値ではなく、学習の習熟度、病気の重症度、好みの程度などのような順位で与えられているとき、$${x, y}$$ の関係を」順位相関係数で「数量的に表現してみましょう」としています。
■ スピアマンの順位相関係数 － 順位による相関係数 p.52
テキストの計算式をお借りします。
$${N}$$ 組のデータ $${A,B}$$ の 順位 $${(a_1,b_1), (a_2,b_2), \ldots, (a_N, b_N)}$$（ただし同順位は無い）が与えられているとき、データ $${A,B}$$ のスピアマンの順位相関係数 $${r_s}$$ は次の式で求められます。
$$
r_s = 1 - \cfrac{6 \left((a_1-b_1)^2+(a_2-b_2)^2 + \cdots + (a_N-b_N)^2 \right)}{N(N^2-1)}
$$
統計WEBさんの記事を参考にして、同順位（タイ）がある場合を含む、スピアマンの順位相関係数の算出関数を作ってみました。
ありがとうございます！
もちろん【Python標準縛り】です。
### スピアマンの順位相関係数関数 ※list1, list2には順位を1から順に設定すること p.52
# 参考サイト 統計WEB: https://bellcurve.jp/statistics/course/26041.html
# スピアマンの順位相関係数の算出関数
def calc_spearman_rank_corr(list1, list2):
## 設定と準備
# 標本サイズ: ２つの標本サイズは同じ
N = len(list1)
# タイ（同順位）の件数カウント
tx = [v for v in collections.Counter(list1).values() if v > 1]
ty = [v for v in collections.Counter(list2).values() if v > 1]
# タイが無い場合の順位相関係数の算出
if (tx == []) & (ty == []):
# テキストの数式に従って、右辺の分母, rsを算出
numerator = 6 * sum((x - y)**2 for x, y in zip(list1, list2))
# 戻り値: スピアマンの順位相関係数rs
return 1 - numerator / (N * (N**2 - 1))
# タイがある場合の順位相関係数の算出: 統計WEBの数式を利用
else:
# データx、データyのいずれかにタイが無い場合、0のリストを設定
tx = tx if tx != [] else [0]
ty = ty if ty != [] else [0]
# 統計WEBの数式に従ってTx, Ty ,rsの右辺の分母, rs（戻り値）を算出
Tx = (N * (N**2 - 1) - sum([t * (t**2 - 1) for t in tx])) / 12
Ty = (N * (N**2 - 1) - sum([t * (t**2 - 1) for t in ty])) / 12
numerator = Tx + Ty - sum((x - y)**2 for x, y in zip(list1, list2))
# 戻り値: スピアマンの順位相関係数rs
return numerator / (2 * (Tx * Ty)**(1/2))p.53 の例、表1.4.2「手先の器用さと習熟度の順位」のデータを登録します。
### 例 手先の器用さと習熟度の順位 p.53 表1.4.2
# x,yをデータフレーム化
data16 = pd.DataFrame(
{'手先の器用さ': [5, 6, 9, 3, 2, 4, 8, 10, 7, 1],
'習熟度': [7, 5, 6, 3, 4, 1, 8, 9, 10, 2]},
index=range(1, 11)
)
# 結果の表示
print('data16.shape: ', data16.shape)
data16.T【実行結果】
p.53 表1.4.3「スピアマンの順位相関係数を求めるために…」のデータを作成します。
### スピアマンの順位相関係数を求めるために･･･ p.53 表1.4.3
# data16をコピー、列名を変更
data16_rank = data16.copy()
data16_rank.columns = ['順位ai', '順位bi']
# ai-bi, (ai-bi)²列の追加
data16_rank['ai-bi'] = data16_rank['順位ai'] - data16_rank['順位bi']
data16_rank['(ai-bi)²'] = data16_rank['ai-bi']**2
# 結果の表示
print('data16_rank.shape: ', data16_rank.shape)
data16_rank【実行結果】
さまざまな方法でスピアマンの順位相関係数を算出します。
① 自作関数利用
### スピアマンの順位相関係数の算出 自作関数
print(f"スピアマンの順位相関係数 = ",
f"{calc_spearman_rank_corr(data16['手先の器用さ'], data16['習熟度']):.4f}")【実行結果】
② p.52の数式利用
### スピアマンの順位相関係数rsの算出 p.52の数式利用
# 標本サイズの算出
N = len(data16_rank)
# スピアマンの順位相関係数の算出
rs = 1 - (6 * data16_rank['(ai-bi)²'].sum()) / (N * (N**2 - 1))
print(f'スピアマンの順位相関係数 = {rs:.4f}')【実行結果】
r にスピアマンの順位相関係数が表示されます。
③ pandas 利用
相関係数を計算する corr() の引数 method に 'spearman'（スピアマン）を与えることで、スピアマンの順位相関係数を算出できます。
### スピアマンの順位相関係数の算出 pandas
data16[['手先の器用さ', '習熟度']].corr(method='spearman')【実行結果】
④ scipy.stats 利用
spearmanr() でスピアマンの順位相関係数を算出できます。
### スピアマンの順位相関係数の算出 scipy.stats
# 帰無仮説「２つの標本に順位相関は無い」
rs, p_value = stats.spearmanr(data16['手先の器用さ'], data16['習熟度'])
print(f'スピアマンの順位相関係数 = {rs:.4f} (p値: {p_value:.4f})')【実行結果】
⑤ pingouin 利用
相関係数を計算する corr() の引数 method に 'spearman'（スピアマン）を与えることで、スピアマンの順位相関係数を算出できます。
### スピアマンの順位相関係数の算出 pingouin
# 帰無仮説「２つの標本に順位相関は無い」
pg.corr(data16['手先の器用さ'], data16['習熟度'], method='spearman')【実行結果】
■ ケンドールの順位相関係数 － ２組の順位の向きによる関係 p.54
ケンドールの順位相関係数は「順位の向き」を用いた相関です。
p.54 表1.4.4「２つのグループの好きなお酒の順位」のデータを作成します。
### 2つのグループの好きなお酒の順位 p.54 表1.4.4
# データフレームの作成
data17 = pd.DataFrame(
{'女子大生': [5, 2, 1, 4, 3], 'OL': [3, 1, 2, 5, 4]},
index=['日本酒', 'ビール', 'ワイン', 'ウィスキー', 'チューハイ'])
# 結果の表示
print('data17.shape: ', data17.shape)
data17.T【実行結果】
順位の向きを少々見てみます。
日本酒とビールに関しては・・・
・女子大生 日本酒 ＞ ビール
・ＯＬ 日本酒 ＞ ビール
順位の向きは同じです。
ビールとワインに関しては・・・
・女子大生 ビール ＞ ワイン
・ＯＬ ビール ＜ ワイン
順位の向きが逆です。
ケンドールの順位相関係数は、「同じ向きの組の数」と「逆の向きの組の数」を使って相関係数を算出します。
テキストのケンドールの順位相関係数 $\tau$ の計算式をお借りします。
$$
\tau = \cfrac{（同じ向きの組の数）-（逆の向きの組の数）}{（すべての組合せの数）}
$$
または、$${N}$$ 組のデータ $${A,B}$$ の 順位 $${(a_i,b_i), (a_j,b_j)}$$ について、同じ向きの組の数 $${P}$$、逆向きの組の数 $${Q}$$ とすると、ケンドールの順位相関係数 $${\tau}$$ は次の式で求められます。
$$
\tau = \cfrac{P-Q}{\cfrac{N(N-1)}{2}} = \cfrac{2(P-Q)}{N(N-1)}
$$
統計WEBさんの記事を参考にして、同順位（タイ）がある場合を含む、ケンドールの順位相関係数の算出関数を作ってみました。
ありがとうございます！
### ケンドールの順位相関係数関数 p.55
# list1, list2には順位が1から順に設定されている
# 参考サイト 統計WEB: https://bellcurve.jp/statistics/course/26041.html
def kendall_rank_corr(list1, list2):
## 設定
# 同じ向きの組P, 逆の向きの組Q, xのタイTx, yのタイTyのカウンターの初期化
P, Q, Tx, Ty = 0, 0, 0, 0
# 比較した組み合わせの総数
n = len(list1) # 標本サイズ
N = n * (n - 1) / 2 # 組み合わせの総数
## 同じ向きの組・逆の向きの組をカウント
# 組み合わせの特定
for i in range(n):
for j in range(i+1, n):
## x,yともに組み合わせ間の差を算出
x_comp = list1[i] - list1[j]
y_comp = list2[i] - list2[j]
## 向きの判定
# x,yの差が両方ともプラス or 両方ともマイナスの時、同じ向きPをカウント
if (x_comp > 0) & (y_comp > 0):
P += 1
elif (x_comp < 0) & (y_comp < 0):
P += 1
# x,yの差が一方がプラス、他方がマイナスの時、逆の向きQをカウント
elif (x_comp > 0) & (y_comp < 0):
Q += 1
elif (x_comp < 0) & (y_comp > 0):
Q += 1
# xがタイのとき、xのタイをカウント
elif x_comp == 0:
Tx += 1
# xがタイのとき、xのタイをカウント
elif y_comp == 0:
Ty += 1
## ケンドールの順位相関係数τの算出
# タイデータが無い場合
if (Tx == 0) & (Ty == 0):
tau = (P - Q) / N
# タイデータがある場合
else:
tau = (P - Q) / ((N - Tx)**(1/2) * (N - Ty)**(1/2))
## 戻り値: ケンドールの順位相関係数τ
return tau
### テスト
kendall_rank_corr(data17['女子大生'].values, data17['OL'].values)【実行結果】
テスト結果 0.4 はテキスト p.55 の計算結果と一致しています。
■ 理解度チェック：ケンドールの順位相関係数 p.57
p.57 表1.4.7 の 20代と40代の海外旅行先の好みデータを作成します。
### 20代女性と40代女性の海外旅行先の好み p.58 表1.4.7
# データフレームの作成
data18 = pd.DataFrame(
{'20代の順位': [2, 6, 4, 1, 3, 5], '40代の順位': [1, 5, 2, 3, 4, 6]},
index=['韓国', '中国', 'イギリス', 'フランス', 'イタリア', 'ロシア']
)
# 結果の表示
print('data18.shape: ', data18.shape)
data18.T【実行結果】
さまざまな方法でケンドールの順位相関係数を算出します。
① 自作関数利用
### ケンドールの順位相関係数の算出 自作関数利用
tau = kendall_rank_corr(data18['20代の順位'].values, data18['40代の順位'].values)
print(f'ケンドールの順位相関係数 = {tau:.6f}')【実行結果】
② pandas 利用
相関係数を計算する corr() の引数 method に 'kendall'（ケンドール）を与えることで、ケンドールの順位相関係数を算出できます。
### ケンドールの順位相関係数の算出 pandas
data18[['20代の順位', '40代の順位']].corr(method='kendall')【実行結果】
③ scipy.stats 利用
kendalltau() でケンドールの順位相関係数を算出できます。
### ケンドールの順位相関係数の算出 scipy.stats
# 帰無仮説「２つの標本に順位相関は無い」
tau, p_value = stats.kendalltau(data18['20代の順位'], data18['40代の順位'])
print(f'ケンドールの順位相関係数 = {tau:.6f} (p値: {p_value:.4f})')【実行結果】
④ pingouin 利用
相関係数を計算する corr() の引数 method に 'kendall'（ケンドール）を与えることで、ケンドールの順位相関係数を算出できます。
### ケンドールの順位相関係数の算出 pingouin
# 帰無仮説「２つの標本に順位相関は無い」
pg.corr(data18['20代の順位'], data18['40代の順位'], method='kendall')【実行結果】
r にケンドールの順位相関係数が表示されます。
同じデータでスピアマンの相関係数を算出してみます。
### 参考：スピアマンの順位相関係数の算出 pingouin
# 帰無仮説「２つの標本に順位相関は無い」
pg.corr(data18['20代の順位'], data18['40代の順位'], method='spearman')【実行結果】
r にスピアマンの順位相関係数が表示されます。
Section 1.5 クロス集計表とオッズ比
■ オッズ p.59
オッズの値が大きくなるほど、出来事 $${A}$$ の起こる確率 $${p}$$ も大きくなります。
$$
オッズ = \cfrac{p}{1-p}
$$
p.59 のオッズから確率 $${p}$$ を算出する関数を作ってみます。
【Python標準縛り】です。
### オッズから確率pを算出する関数 p.59
def calc_p_from_odds(odds):
return odds / (1 + odds)
# テスト
calc_p_from_odds(2)【実行結果】
■ オッズ比 p.60
p.60 のオッズ比の定義の数式を用いて、オッズ比算出関数を【Python標準縛り】で作成します。
### オッズ比算出関数 p.60
def calc_odds_ratio(p, q):
return (p / (1 - p)) / (q / (1 - q))
# テスト
calc_odds_ratio(p=0.5, q=0.3)【実行結果】
■ オッズ比と事象の独立との関係 p.61
テキストによると、オッズ比を用いることで、事象 $${A,B}$$ が独立かどうかを判定できるそうです。
$$
オッズ比が１ \Longleftrightarrow AとBは独立である
$$
■ 理解度チェック：オッズ比 p.63
p.63 表1.5.4「自動車事故の死傷者に関する調査結果」データを作成します。
### 自動車事故の死傷者に関する調査結果 p.63 表1.5.4
data19 = pd.DataFrame(
{'死亡した人': [167, 31], '死亡しなかった人': [8896, 10533]},
index=['シートベルト未着用', 'シートベルト着用']
)
print('data19.shape: ', data19.shape)
data19.T【実行結果】
オッズとオッズ比を算出します。
### オッズ、オッズ比の算出 p.63
# シートベルト未着用で死亡する事象A, シートベルト着用で死亡する事象Bとして確率p,qを算出
p, q = data19['死亡した人'] / data19.sum(axis=1)
print(f'シートベルト未着用者が死亡した確率 p = {p:.4f}\n'
f'シートベルト着用者が死亡した確率 q = {q:.4f}')
# オッズ
odds_a = p / (1 - p)
odds_b = q / (1 - q)
print(f'シートベルト未着用者のオッズ = {odds_a:.4f}')
print(f'シートベルト着用者のオッズ = {odds_b:.4f}')
# オッズ比
print(f'オッズ比 = {odds_a / odds_b:.4f}')【実行結果】
うさぎの「このオッズ比は次のようにしても求まります」を実践します。
### テキストのオッズ比の別解 p.63
(167 * 10533) / (31 * 8896)【実行結果】
今回の写経は以上です。
シリーズの記事
次の記事
前の記事
目次
ブログの紹介
note で７つのシリーズ記事を書いています。
ぜひ覗いていってくださいね！
１．のんびり統計
統計検定２級の問題集を手がかりにして、確率・統計をざっくり掘り下げるブログです。
雑談感覚で大丈夫です。ぜひ覗いていってくださいね。
統計検定２級公式問題集CBT対応版に対応しています。
Python、EXCELのサンプルコードの配布もあります。
２．実験！たのしいベイズモデリング１＆２をPyMC Ver.5で
書籍「たのしいベイズモデリング」・「たのしいベイズモデリング２」の心理学研究に用いられたベイズモデルを PyMC Ver.5で描いて分析します。
この書籍をはじめ、多くのベイズモデルは Ｒ言語＋Stanで書かれています。
PyMCの可能性を探り出し、手軽にベイズモデリングを実践できるように努めます。
身近なテーマ、イメージしやすいテーマですので、ぜひぜひPyMCで動かして、一緒に楽しみましょう！
３．実験！岩波データサイエンス1のベイズモデリングをPyMC Ver.5で
書籍「実験！岩波データサイエンスvol.1」の４人のベイジアンによるベイズモデルを PyMC Ver.5で描いて分析します。
この書籍はベイズプログラミングのイロハをざっくりと学ぶことができる良書です。
楽しくPyMCモデルを動かして、ベイズと仲良しになれた気がします。
みなさんもぜひぜひPyMCで動かして、一緒に遊んで学びましょう！
４．楽しい写経 ベイズ・Python等
ベイズ、Python、その他の「書籍の写経活動」の成果をブログにします。
主にPythonへの翻訳に取り組んでいます。
写経に取り組むお仲間さんのサンプルコードになれば幸いです🍀
５．RとStanではじめる心理学のための時系列分析入門 を PythonとPyMC Ver.5 で
書籍「RとStanではじめる心理学のための時系列分析入門」の時系列分析をPythonとPyMC Ver.5 で実践します。
この書籍には時系列分析のテーマが盛りだくさん！
時系列分析の懐の深さを実感いたしました。
大好きなPythonで楽しく時系列分析を学びます。
６．データサイエンスっぽいことを綴る
統計、データ分析、AI、機械学習、Pythonのコラムを不定期に綴っています。
統計・データサイエンス書籍にまつわる記事が多いです。
「統計」「Python」「数学とPython」「R」のシリーズが生まれています。
７．Python機械学習プログラミング実践記
書籍「Python機械学習プログラミング PyTorch & scikit-learn編」を学んだときのさまざまな思いを記事にしました。
この書籍は、scikit-learnとPyTorchの教科書です。
よかったらぜひ、お試しくださいませ。
最後までお読みいただきまして、ありがとうございました。
いいなと思ったら応援しよう！
応援ありがとうございます。これからもがんばって記事を作成します！