【統計学入門（東京大学出版会）】第3章練習問題解答

前回に引き続き、東京大学出版会から出版されている統計学入門（基礎統計学Ⅰ）について第3章の練習問題の解答を書いていきます。

本章以外の解答
3.1
3.3
3.4
- (i)
- (ii)
- (iii)

本章以外の解答

本章以外の練習問題の解答は別の記事で公開しています。
必要に応じて参照してください。

3.1

自民得票率（tokuhyo）をX軸に、持ち家比率（motiie）をY軸にした散布図を描くPythonプログラムを次に示します。
matplotlibの図の背景はデフォルトで透明であるため、fig.patch.set_alpha(1) で不透明にしている点に注意してください。

import numpy as np
import matplotlib.pyplot as plt

tokuhyo = np.array([
    41.4, 76.3, 59.2, 51.8, 52.5, 53.2, 62.4, 55.0, 57.7,
    63.2, 37.5, 48.5, 32.4, 20.5, 47.9, 68.9, 68.5, 52.5, 63.3, 58.8,
    59.7, 48.4, 40.7, 51.0, 50.9, 34.3, 25.8, 32.1, 34.4, 55.1, 60.3,
    57.0, 45.6, 54.2, 55.1, 55.7, 70.3, 61.8, 47.6, 42.5, 71.3, 55.2,
    65.2, 42.9, 54.7, 62.0, 48.2
])
motiie = np.array([
    52.8, 71.2, 72.6, 63.7, 81.3, 81.8, 70.9, 74.0, 73.2,
    72.9, 66.7, 65.7, 43.7, 55.5, 79.6, 85.7, 75.3, 80.5, 73.0, 77.0,
    77.5, 69.2, 60.0, 78.2, 79.5, 61.8, 49.6, 59.6, 72.1, 71.0, 76.3,
    72.8, 71.8, 60.7, 67.0, 71.8, 71.2, 68.3, 68.5, 54.8, 76.0, 65.8,
    69.4, 66.9, 69.7, 71.2, 59.6
])

fig = plt.figure(figsize=(6, 6))
fig.patch.set_alpha(1)

plt.subplot(1, 1, 1)
plt.scatter(tokuhyo, motiie)
plt.xlabel("tokuhyo")
plt.ylabel("motiie")
plt.title("3.1")

上記のプログラムを実行すると、次の散布図が描画されます。

f:id:nuka137:20200822151604p:plain

変数 $x, y$ の間の相関係数 $r_{xy}$ は、次の式によって求められます。

$\displaystyle r_{xy} = \frac{\sum (x_{i} - \bar{x})(y_{i} - \bar{y})}{\sqrt{\sum (x_{i} - \bar{x})^{2}}\sqrt{\sum (y_{i} - \bar{y})^{2}}}$

変数 $x$ を自民得票率、変数 $y$ を持ち家比率として相関係数 $r_{xy}$ を求めるPythonプログラムを次に示します。

import numpy as np

tokuhyo = np.array([
    41.4, 76.3, 59.2, 51.8, 52.5, 53.2, 62.4, 55.0, 57.7,
    63.2, 37.5, 48.5, 32.4, 20.5, 47.9, 68.9, 68.5, 52.5, 63.3, 58.8,
    59.7, 48.4, 40.7, 51.0, 50.9, 34.3, 25.8, 32.1, 34.4, 55.1, 60.3,
    57.0, 45.6, 54.2, 55.1, 55.7, 70.3, 61.8, 47.6, 42.5, 71.3, 55.2,
    65.2, 42.9, 54.7, 62.0, 48.2
])
motiie = np.array([
    52.8, 71.2, 72.6, 63.7, 81.3, 81.8, 70.9, 74.0, 73.2,
    72.9, 66.7, 65.7, 43.7, 55.5, 79.6, 85.7, 75.3, 80.5, 73.0, 77.0,
    77.5, 69.2, 60.0, 78.2, 79.5, 61.8, 49.6, 59.6, 72.1, 71.0, 76.3,
    72.8, 71.8, 60.7, 67.0, 71.8, 71.2, 68.3, 68.5, 54.8, 76.0, 65.8,
    69.4, 66.9, 69.7, 71.2, 59.6
])

n = tokuhyo.size
tokuhyo_ave = np.mean(tokuhyo)
motiie_ave = np.mean(motiie)
rxy = np.sum((tokuhyo - tokuhyo_ave) * (motiie - motiie_ave)) / (np.sqrt(np.sum((tokuhyo - tokuhyo_ave)**2)) * np.sqrt(np.sum((motiie - motiie_ave)**2)))
print(f"Correlation Coefficient: {rxy}")

上記プログラムの実行結果を次に示します。

Correlation Coefficient: 0.6387815655787518

3.3

スピアマンの順位相関係数 $r_{s}$ は、次式で求められます。

$\displaystyle r_{s} = 1 - \frac{6}{n^{3} -n} \sum_{i = 1}^{n} (R_{i} - R_{i}^{'})^{2}$

一方、ケンドールの順位相関係数 $r_{k}$ は、次式で求められます。

$\displaystyle r_{k} = \frac{G - H}{n (n - 1) / 2}$

$G$ は、 $R_{i} \gt R_{j}$ かつ $R_{i}^{'} \gt R_{j}^{'}$ または、 $R_{i} \lt R_{j}$ かつ $R_{i}^{'} \lt R_{j}^{'}$ を満たす個数。
$H$ は、 $R_{i} \gt R_{j}$ かつ $R_{i}^{'} \lt R_{j}^{'}$ または、 $R_{i} \lt R_{j}$ かつ $R_{i}^{'} \gt R_{j}^{'}$ を満たす個数。

これらの式に従って $r_{s}$ と $r_{k}$ を求めるPythonプログラムを次に示します。

import numpy as np
import itertools


def spearman(d1, d2):
    n = d1.size
    return 1 - 6 * np.sum((d1 - d2)**2) / (n**3 - n)

def kendall(d1, d2):
    G = H = 0
    n = d1.size
    for i in range(n):
        Ri = d1[i]
        Rid = d2[i]
        for j in range(i+1, n):
            Rj = d1[j]
            Rjd = d2[j]
            if (Ri > Rj and Rid > Rjd) or (Ri < Rj and Rid < Rjd):
                G += 1
            elif (Ri > Rj and Rid < Rjd) or (Ri < Rj and Rid > Rjd):
                H += 1

    return (G - H) / (n*(n-1) / 2)

data = {
    "1": np.array([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]),
    "2": np.array([1, 5, 2, 3, 6, 7, 15, 8, 4, 11, 10, 14, 18, 13, 22, 24, 16, 19, 30, 9, 25, 17, 26, 23, 12, 20, 28, 21, 27, 29]),
    "3": np.array([8, 3, 1, 4, 2, 5, 11, 7, 15, 9, 6, 13, 10, 22, 12, 14, 18, 19, 17, 22, 16, 24, 21, 20, 28, 30, 25, 26, 27, 29]),
    "4": np.array([20, 1, 4, 2, 6, 3, 12, 17, 8, 5, 18, 13, 23, 26, 29, 15, 16, 9, 10, 11, 30, 7, 27, 19, 14, 21, 28, 24, 22, 25]),
}

for d in itertools.combinations(data.keys(), 2):
    entry1 = d[0]
    entry2 = d[1]
    rs = spearman(data[entry1], data[entry2])
    rk = kendall(data[entry1], data[entry2])
    print(f"rs ({entry1} - {entry2}): {rs}")
    print(f"rk ({entry1} - {entry2}): {rk}")
    print("")

上記プログラムの実行結果を次に示します。

rs (1 - 2): 0.821579532814238
rk (1 - 2): 0.664367816091954

rs (1 - 3): 0.9272525027808676
rk (1 - 3): 0.767816091954023

rs (1 - 4): 0.5933259176863181
rk (1 - 4): 0.43908045977011495

rs (2 - 3): 0.671857619577308
rk (2 - 3): 0.5011494252873563

rs (2 - 4): 0.6373748609566184
rk (2 - 4): 0.46206896551724136

rs (3 - 4): 0.5343715239154616
rk (3 - 4): 0.3586206896551724

3.4

それぞれPythonのプログラムで解答を書きます。

(i)

乱数を発生させるPythonプログラムを示します。
1から11の乱数を発生させるために、random.randint を利用しました。

import random

random.seed(1)
d = int(random.randint(1, 11))

(ii)

変数 $x, y$ の間の相関係数 $r_{xy}$ は、次の式によって求められます。

$\displaystyle r_{xy} = \frac{\sum (x_{i} - \bar{x})(y_{i} - \bar{y})}{\sqrt{\sum (x_{i} - \bar{x})^{2}}\sqrt{\sum (y_{i} - \bar{y})^{2}}}$

相関係数 $r_{xy}$ を求めるPythonプログラムを次に示します。

import random
import numpy as np

random.seed(1)

x = np.array([71, 68, 66, 67, 70, 71, 70, 73, 72, 65, 66])
y = np.array([69, 64, 65, 63, 65, 62, 65, 64, 66, 59, 62])
x_mean = np.mean(x)
y_mean = np.mean(y)

x_sampled = np.zeros_like(x)
y_sampled = np.zeros_like(y)
for i in range(11):
    d = int(random.randint(1, 11)) - 1
    x_sampled[i] = x[d]
    y_sampled[i] = y[d]

rxy = np.sum((x_sampled - x_mean) * (y_sampled - y_mean)) / \
    (np.sqrt(np.sum((x_sampled - x_mean)**2)) * np.sqrt(np.sum((y_sampled - y_mean)**2)))
print(rxy)

上記プログラムの実行結果を次に示します。

0.49543369430686224

(iii)

(ii)を200回行い、 $r_{xy}$ の値をビン数10でヒストグラムにするPythonプログラムを次に示します。

import random
import numpy as np
import matplotlib.pyplot as plt

random.seed(1)

x = np.array([71, 68, 66, 67, 70, 71, 70, 73, 72, 65, 66])
y = np.array([69, 64, 65, 63, 65, 62, 65, 64, 66, 59, 62])
x_mean = np.mean(x)
y_mean = np.mean(y)

rxy_list = []
for _ in range(200):
    x_sampled = np.zeros_like(x)
    y_sampled = np.zeros_like(y)
    for i in range(11):
        d = int(random.randint(1, 11)) - 1
        x_sampled[i] = x[d]
        y_sampled[i] = y[d]

    rxy = np.sum((x_sampled - x_mean) * (y_sampled - y_mean)) / \
        (np.sqrt(np.sum((x_sampled - x_mean)**2)) * np.sqrt(np.sum((y_sampled - y_mean)**2)))
    rxy_list.append(rxy)

fig = plt.figure(figsize=(6, 6))
fig.patch.set_alpha(1)

plt.subplot(1, 1, 1)
plt.hist(rxy_list, bins=10)
plt.xlabel("rxy")
plt.ylabel("count")
plt.title("3.4 (iii)")

上記のプログラムを実行すると、次のような図が表示されます。

f:id:nuka137:20200822184620p:plain