生物工学会誌 –『間違いから学ぶ実践統計解析』をPythonで実装する #1

[『間違いから学ぶ実践統計解析』PDF版](https://www.sbj.or.jp/sbj/sbj_tokei_kaiseki.html) のRで書かれているコードをPythonで実装していきます。

実践統計解析【第1回】
---

```Python
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

# データ準備
data_A = np.array([15, 13, 11, 13, 18, 19, 22, 21, 20, 16, 11, 16, 18, 10, 11, 12, 11, 14, 12, 11])
data_B = np.array([15, 17, 21, 23, 28, 19, 22, 24, 30, 26, 21, 19, 18,20, 21, 22, 21, 24, 22, 21])

# 平均値を計算する
print('Aさんのデータの平均値：', np.mean(data_A))
print('Bさんのデータの平均値：', np.mean(data_B))
```
```
Aさんの平均値： 14.7
Bさんの平均値： 21.7
```
<br>

```Python
# Aさんのデータのヒストグラムを作成する
bins = list(range(np.min(data_A), np.max(data_A) + 1, 2))
freq = pd.Series(data_A).value_counts(bins=bins).values
plt.bar(bins[:-1], freq, width=np.diff(bins), align='edge', edgecolor='black', color='w')
plt.yticks(range(0, np.max(freq) + 1, 2))
plt.title('Histogram of A')
plt.xlabel('A')
plt.ylabel('Frequency')
plt.show()
```
![histgramA.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/469902/024c6cc5-d5c2-603d-516c-045d94582e01.png)

- テキストに合わせてビン幅を「2」とした
- 各ビン範囲が (a, b] となるよう[pandas.Series.value_counts](https://pandas.pydata.org/docs/reference/api/pandas.Series.value_counts.html)を使用した
- 見た目を調整した
<br>

```Python
# Bさんのデータのヒストグラムを作成する
value_count = pd.Series(data_B).value_counts(bins=8, sort=False)
bins = value_count.index.values.left
freq = value_count.values
width = np.diff(bins)[0]
plt.bar(bins, freq, width=width, align='edge', edgecolor='black', color='w')
plt.xticks(range(np.min(data_B), np.max(data_B)+1, 5))
plt.yticks(range(0, np.max(freq) + 1, 2))
plt.title('Histogram of B')
plt.xlabel('B')
plt.ylabel('Frequency')
plt.show()
```
![histgramB.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/469902/2bc8b6a7-614f-adc0-17bf-146823fa5c33.png)

- テキストに合わせてビン数を「8」とした
- 各ビン範囲が (a, b] となるよう[pandas.Series.value_counts](https://pandas.pydata.org/docs/reference/api/pandas.Series.value_counts.html)を使用した
- 見た目を調整した
<br>

```Python
# Aさんのデータの要約統計量を確認する
pd.Series(data_A).describe()
```
```
count    20.000000
mean     14.700000
std       3.812721
min      10.000000
25%      11.000000
50%      13.500000
75%      18.000000
max      22.000000
```

```Python
# Bさんのデータの要約統計量を確認する
pd.Series(data_B).describe()
```
```
count    20.000000
mean     21.700000
std       3.555574
min      15.000000
25%      19.750000
50%      21.000000
75%      23.250000
max      30.000000
```
<br>

```Python
# 箱ひげ図を作成する
plt.boxplot([data_A, data_B])
plt.yticks(range(np.min([np.min(data_A), np.min(data_B)]), np.max([np.max(data_A), np.max(data_B)])+1, 5))
plt.show()
```
![boxplot.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/469902/74d40adb-1028-e947-127a-a659351d9fb5.png)
<br>

```Python
# 不偏分散を計算する

# Aさんのデータ
print(np.var(data_A, ddof=1))
# Bさんのデータ
print(np.var(data_B, ddof=1))
```
```
14.53684210526316
12.642105263157896
```
<br>

```Python
# 不偏標準偏差を計算する

# Aさんのデータ
print(np.std(data_A, ddof=1))
# Bさんのデータ
print(np.std(data_B, ddof=1))
```
```
3.8127210893616597
3.5555738303623925
```
<br>

```Python
# 標準誤差を計算する

# Aさんのデータ
print(np.sqrt(np.var(data_A, ddof=1) / data_A.size))
# Bさんのデータ
print(np.sqrt(np.var(data_B, ddof=1) / data_B.size))
```
```
0.8525503535059721
0.7950504783709615
```

実行環境
---
&emsp;Python 3.9.12
&emsp;numpy 1.21.5
&emsp;pandas 1.4.2
&emsp;matplotlib 3.5.1