高校数学を復習してデータ分析をはじめよう 〜データの分析編〜

はじめに

こんにちは。高校2年の樅山です。
2020/11 から始まった、ものづくりをする高校生のための新しいグループ、「Palettte」が主催する Palettte Advent Calendar 2020 の19日目の記事となります。

高校数学では、文系理系に関わらず多くの生徒が履修する科目、「数学I」ですが、その中に「データの分析」という単元があるのをご存知でしょうか。
この単元は、統計の基礎の基礎と言ってもよく、データを扱う上での基本的な事柄を抑えることができます。

しかし、数学という科目の中で扱うと、実際にグラフが表示されるわけでもなく、視覚的にはかなり退屈な（私見です）思い出があります。
そこで、データ分析に用いられるR言語を使って、実際にプログラムを書きながら学んでいけば楽しいかなと思い、この記事を書きました。

R言語

Rとは、統計の図表や関数などが提供された統合開発環境のことを指します。
そして、それに用いられるR言語は、C言語などの文法に似た、インタプリタ言語です。

R言語の代入

a <- 100
a = 100

R言語には、この2つの代入方法があります。
本記事では、どちらも利用するため、ご注意ください。

データの整理

データ

ある集団の特徴に関する値を数量的に表現したものを、データといいます。

4/1~4/16までの東京都の陽性患者数(人)

66  97  89  116
143 83  79  144
178 188 197 166
91  161 126 148 

開発でも往々にして扱うような、それこそデータベースに格納したくなるような集合が、データです。
これには、16日分の測定値が含まれています。この測定値の個数のことを、データの大きさ と呼びます。

東京都_新型コロナウイルス陽性患者発表詳細 でダウンロード可能なcsvファイルを利用して、現在までの東京都の陽性患者数のデータを取得してみます。

data.R

data = read.csv("https://stopcovid19.metro.tokyo.lg.jp/data/130001_tokyo_covid19_patients.csv")
data["公表_年月日"]

     公表_年月日
1     2020-01-24
2     2020-01-25
3     2020-01-30
4     2020-02-13
5     2020-02-14

... (省略)

3080  2020-04-19
3081  2020-04-19
3082  2020-04-19

このように、インターネット上のcsvデータを取得して、R上で扱えるようになりました。
また、data["公表_年月日"] の実行時点でデータの大きさもわかりますが、

data.R

data = read.csv("https://stopcovid19.metro.tokyo.lg.jp/data/130001_tokyo_covid19_patients.csv")
dim(data) # -> 3082   16

このcsvデータは 行列 という、いわば二次元配列のようなもので管理されているので、dim関数という行列の大きさを調べる関数を使えば、3082人分のデータがあるということがわかります。
つまり、データの大きさは 3082 です。

度数分布表

データの散らばりの様子を分布といいます。データの分布を見るための1つの方法として、度数分布表があります。

これが度数分布表です。Rで出力する方法がなかったので作りました。
度数分布表において、区切られた区間のことを 階級 、各区間の幅のことを 階級の幅 、各階級に入るデータの値の個数を 度数 、各階級の真ん中の値を 階級値 といいます。

ヒストグラム

データを柱状グラフで表現したものを、ヒストグラムといいます。
これをR言語でプロットしてみましょう。

histogram.R

person <- data["公表_年月日"]
plot(person)

実行すると、このようなウインドウが開きます。

柱はそれぞれ、各階級の度数に比例しています。これも度数分布表と同じようにデータの分布を見るための方法の1つですが、より視覚的に理解することができます。

データの代表値

データ全体の特徴を表す、適当なある1つの数値のことをデータの代表値といいます。
今回は、平均値、中央値、最頻値に注目して扱います。

前準備

これからデータの代表値を取得するために、csvファイルを整形して、何日に何人の陽性患者が公表されているかどうかを示すベクトルを作ります。
これ以降では、csvファイルを取得して得た data 変数と、何日に何人の陽性患者が公表されているか示す date_paticients を定義なしに利用します。

平均値

あるデータの大きさが $n$ 個、値が $x_1, x_2, x_3, ... x_n$ であるとします。
この時、値の総和を $n$ で割ったものを、データの平均値といいます。
つまり、平均値 $\bar{x}$ は

$$\bar{x} = \frac{1}{n}\sum_{k=1}^{n} x_k$$

と定義されます。

多くのプログラムでも、平均値を導出することはよくあると思います。
これをR上で表現します。

average.R

data = c(5.8, 5.5, 7.0, 5.2, 4.6, 6.1, 3.9, 4.2, 5.7, 7.1, 4.9, 6.7, 6.4, 5.8, 4.9)
mean(data) # -> 5.586667

Rでは、mean関数という組み込み関数を使うことで平均値を得ることができます。

中央値

あるデータを、値の大きさ順に並べ（=ソートし)、中央の位置にくる値を、中央値、またはメジアンと言います。データの大きさが偶数の場合は、中央に並ぶ2つの値の平均値を中央値とします。

median.R

data = c(5.8, 5.5, 7.0, 5.2, 4.6, 6.1, 3.9, 4.2, 5.7, 7.1, 4.9, 6.7, 6.4, 5.8, 4.9)
median(data) # -> 5.7

Rでは、median関数という組み込み関数を使うことで中央値を得ることができます。

最頻値

データの中で、最も個数の多い値のことを、最頻値、またはモードと言います。
先ほどの、度数分布表に整理した時は、度数が最も大きい階級の階級値のことを最頻値といいます。

Rでは、最頻値を求める組み込み関数は用意されていないので、table関数を使ってテーブルを出力して、その中から最頻値を絞り込みます。

table.R

data = c(5.8, 5.5, 7.0, 5.2, 4.6, 6.1, 3.9, 4.2, 5.7, 7.1, 4.9, 6.7, 6.4, 5.8, 4.9)
table(data)

出力結果

3.9 4.2 4.6 4.9 5.2 5.5 5.7 5.8 6.1 6.4 6.7   7 7.1 
  1   1   1   2   1   1   1   2   1   1   1   1   1

mord.R

table_data <- table(data)
table_data[table_data == max(table_data)]

# 出力結果
# 4.9 5.8 
# 2   2 

データの散らばりと四分位範囲

範囲

データの最大値と最小値の差のことを、範囲 といいます。
範囲は、データがどれくらい散らばっているか示す1つの値です。

では、Rでデータの範囲を調べてみます。

range.R

data = c(5.8, 5.5, 7.0, 5.2, 4.6, 6.1, 3.9, 4.2, 5.7, 7.1, 4.9, 6.7, 6.4, 5.8, 4.9)
max(data) - min(data) # -> 3.2

max関数と、min関数を使って、データの最大値と最小値を調べ、その差をとっています。

四分位数

範囲は、データの散らばりを示す1つの値でした。
しかし、これは最大値と最小値のみから計算しているので、例えば

Data1

1 10 10 10 10 10 10 10 10 10 10 10 101

Data2

1 10 20 30 40 50 60 70 80 90 100 101

この2つのデータはどちらも範囲が100です。しかし、この2つのデータは本当に同じだけ散らばっているのでしょうか。
最小値、最大値のみで考えるのではない、他の散らばりの表し方を考えてみましょう。
まず、データを昇順にソートして、4等分し、中央の50%のデータについて考えます。
この時、4等分する位置にある値を、四分位数といいます。
四分位数は、小さい方から、第1四分位数、第2四分位数、第3四分位数と呼んでいますが、長いので、 $Q_1, Q_2, Q_3$ で表しています。
また、これらを
$Q_1$ = 下位データの中央値
$Q_2$ = 中央値
$Q_3$ = 上位データの中央値
と定めます。

さて、Rで四分位数を計算してみます。

summary.R

data1 = c(1, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 101)
summary(data1)
data2 = c(1, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100, 101)
summary(data2)

Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
1.00   10.00   10.00   16.31   10.00  101.00 

Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
1.00   27.50   55.00   54.33   82.50  101.00

summary関数を使うと、最小値、最大値、平均値、中央値(第2四分位数)、そして第1四分位数、第3四分位数を求めることができます。
範囲は同じ値となるデータですが、かなり性質が異なることが見て取れると思います。

四分位範囲・四分位偏差

第3四分位数から第1四分位数を引いたものを 四分位範囲 といいます。

$$Q_3 - Q_1$$

四分位範囲はデータの中央50%の範囲とほとんど等しいので、飛び値 と言われるような、極端に小さい/大きい値の影響を受けにくいです。
四分位範囲を2で割った値を、 四分位偏差 といいます。

$$\frac{Q_3-Q_1}{2}$$

四分位範囲や四分位偏差は、大きいほど散らばっていて、小さいほど散らばっていないことを示しています。

では、Rで四分位範囲・四分位偏差を求めてみます。

quartile.R

data1 = c(1, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 101)
data2 = c(1, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100, 101)

getQuartileRange = function(data) {
  firstQuartile = summary(data)[['1st Qu.']]
  thirdQuartile = summary(data)[['3rd Qu.']]
  thirdQuartile - firstQuartile
}

getQuartileDeviation = function(data){
  getQuartileRange(data) / 2
}

getQuartileRange(data1) # -> 1.35
getQuartileRange(data2) # -> 55

getQuartileRange(data1) # -> 0.675
getQuartileRange(data2) # -> 27.5

summary関数から第1四分位数、第3四分位数を取得できるので、それを計算しています。

箱ひげ図

箱ひげ図ならば、なんとなく聞いたことや、習った覚えがあるのではないでしょうか？
データの分布を見るための図で、データの最小値、第1四分位数、中央値、第3四分位数、最大値を箱と線で表現しています。
つまり、先ほどの summary関数を視覚的に見ることができます。

では、Rでプロットしてみましょう。

BoxBeard.R

data2 = c(1, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100, 101)
boxplot(data2, main="data2")

重要な情報がひと目で理解できるので、summary関数の出力結果を見るよりも、わかりやすいと思います。
では、data1も一緒にプロットしてみましょう。

BoxBeard.R

data1 = c(1, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 101)
data2 = c(1, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100, 101)
boxplot(data1, data2, main="data")

data1の箱ひげ図がとんでもなく薄くプロットされてしまいました。
これを見れば一目瞭然ですが、最大値・最小値のみに依存して導出するのではなく、四分位数を用いて箱ひげ図にプロットすることで、データの散らばりをより正確に読み取ることができます。

分散と標準偏差

変量 $x$ についてのデータの値が、$x_1, x_2, x_3, ..., x_n$ で、平均値を $\bar{x}$ とします。
この時、 $x_1-\bar{x}, x_2-\bar{x}, x_3-\bar{x}, ..., x_n-\bar{x}$ をそれぞれ、 $x_1, x_2, x_3, ..., x_n$ の平均値からの 偏差 といいます。
ちなみに、偏差の平均値は、データの大きさを $n$ とする時、

$$\begin{aligned} \frac{1}{n}((x_1-\bar{x})+(x_2-\bar{x})+(x_3-\bar{x})+...+(x_n-\bar{x})) &= \frac{1}{n}((x_1+x_2+x_3+...+x_n)-n\bar{x})\\ &= \frac{1}{n}(x_1+x_2+x_3+...x_n)-\bar{x}\\ &= \bar{x}-\bar{x}\\ &= 0 \end{aligned}$$

と式変形できるので、常に0 になります。
つまり、偏差の平均値をとっても、データの散らばり度合いを表す事はできません。
そこで、まず偏差の2乗をとって、その平均値を考えることにします。つまり、

$$\frac{1}{n}((x_1-\bar{x})^2+(x_2-\bar{x})^2+(x_3-\bar{x})^2+...+(x_n-\bar{x})^2)$$

を考えます。
これを、データの 分散 と呼び、 $s^2$ で表現します。
分散が0である時、$x_1=x_2=x_3=...=x_n=\bar{x}$ が成り立つので、データの値は全て等しく、平均値にも等しいことがわかります。
分散が小さいという事は、データの平均値周辺の散らばり方が小さいことの目安になります。

しかし、これでは不都合な場合もあります。例えば、クラスメイトの身長のデータを用いて分散を算出したとします。
身長の単位はメートル($m$)ですが、分散の定義は $s^2$ でしたので、単位は平方メートル($m^2$)になってしまいます！

そこで、データの測定単位と同じ単位に直すために、$\sqrt(s^2)$ を散らばりの指標として使うことにしましょう。
この値を、おなじみ 標準偏差 といいます。
分散も標準偏差も、データの散らばりの指標の1つとして用いられます。

それでは、分散と標準偏差をRで算出してみましょう。使うのは、以下のデータです。

部署Aの社員の出社にかかる時間（分）

22 56 19 22 87 35 65 12 44 9 18 47

variance.R

data = = c(22,56,19,22,87,35,65,12,44,9,18,47)

getVariance <- function(data){
  sum = 0
  ave = mean(data)
  for(i in data) {
    sum <- sum + (i-ave)^2
  }
  sum / length(data)
}

getStandardDeviation <- function(data){
  sqrt(getVariance(data))
}

getVariance(data) # -> 528.0556
getStandartDeviation(data) # -> 22.97946

定義通りに求めました。
データ1つからは分散・標準偏差の値をみてもどういった性質なのかまだ理解し難いですが、複数のデータの分散・標準偏差の値を比較することで、散らばりについて考察することができるようになります。

さて、先ほどの分散の式をさらに整理してみましょう。

$$\begin{aligned} s^2 &= \frac{1}{n}((x_1-\bar{x})^2+(x_2-\bar{x})^2+(x_3-\bar{x})^2+...+(x_n-\bar{x})^2)\\ &= \frac{1}{n}((x_1^2+x_2^2+x_3^2+...+x_n^2)-2\bar{x}(x_1+x_2+x_3+...+x_n)+n(\bar{x})^2)\\ &= \frac{1}{n}(x_1^2+x_2^2+x_3^2+...+x_n^2) - 2\bar{x}\frac{1}{n}(z_1+x_2+x_3+...+x_n)+(\bar{x})^2\\ &= \bar{x^2}-2\bar{x}^2+(\bar{x})^2\\ &= x^2-(\bar{x})^2 \end{aligned}$$

丁寧に整理すると、なんと $x$のデータの分散 は、 $x^2$ のデータの平均値 から $x$のデータの平均値の2乗 であることがわかりました！
式を整理することで、より簡単に分散を導出できるようになりました。

データの相関

データを1つだけ扱うのではなく、2つ、3つ、たくさんのデータを比較したい時が往々にしてあります。
まずは、2つのデータの関係について調べる方法を学んでいきましょう！

散布図

まずは、このRスクリプトを実行してみましょう。

scatterPlot.R

height = c(181, 179, 167, 190, 166, 155, 162, 168, 164, 170, 181)
weight = c(99, 58, 68, 85, 57, 46, 60, 59, 60, 71, 77)

plot(height, weight)

すると、このような図表が表示されたのではないでしょうか。
このような、2つのデータについての関係を表現した図を 散布図 といいます。

相関関係

2つのデータについて、一方が増えると他方も増える傾向が認められる時、2つのデータの間に 正の相関関係 があるといいます。
逆に、一方が増えると他方が減る傾向が認められる時、負の相関関係 があるといいます。
どちらの傾向も認められない時、相関関係がない といいます。

例えば、先ほどの図では、外れ値が少しあるものの、正の相関関係があるといえると思います。

scatterPlot2.R

height2 = c(170,172,174,176,178,180,182,184,186,188,190)
weight2 = c(55, 56, 58, 59, 62, 65, 66, 66, 68, 74, 79)

plot(height2, weight2)

一方、新しく用意したこのデータで散布図を描画すると、よりわかりやすく 正の相関関係 が発見できると思います。
この傾向は、height, weight よりも、 height2, weight2 の方が顕著なので、より強い正の相関関係 があるということができます。

正の相関関係がある時、散布図の点は全体的に右上がりに分布します。
逆に、負の相関関係がある時、散布図の点は全体的に右下がりに分布します。

相関係数

データの散らばりを表現する値として 分散 や 標準偏差 を定義したように、相関関係を表現する値を定義してみましょう。

まずは、2つのデータが

$$(x_1, y_1), (x_2, y_2), (x_3, y_3), ..., (x_n, y_n)$$

のように与えられていたとします。
また、平均値をそれぞれ $\bar{x}, \bar{y}$ 、標準偏差を $s_x, s_y$ とします。
ここで、 $x$ の偏差と $y$ の偏差の積 $(x_i-\bar{x})(y_i-\bar{y})$ の平均値である、

$$\frac{1}{n}((x_1-\bar{x})(y_1-\bar{y})+(x_2-\bar{x})(y_2-\bar{y})+(x_3-\bar{x})(y_3-\bar{y})+...+(x_n-\bar{x})(y_n-\bar{y}))$$

を考えます。
これを、$x$ と $y$ の 共分散 といい、 $s_xy$ と表現します。

例えば、正の相関がある データについて考えてみます。
先ほどのデータについて思い出してください。

Data2

height2 = c(170,172,174,176,178,180,182,184,186,188,190)
weight2 = c(55, 56, 58, 59, 62, 65, 66, 66, 68, 74, 79)

height2の平均値は、180 、weight2の平均値は、64.36364 です。
(170, 55) の人は、平均値 (180, 64.363) をどちらも下回ります。
(172, 56) の人も、(174,58)の人も、(176, 59)の人も、(178, 62)の人も同様です。
(180, 65) の人は、体重のみ少し上回りますが、ほとんど平均といって構いません。
(182, 66) の人は、平均値 (180, 64.363) をどちらも上回ります。残りの人も同様です。
すなわち、(i番目の人の身長 - 身長の平均値)(i番目の人の体重 - 体重の平均値) は、 正×正 か、負×負 になり、ほとんどが正の数になります。
つまり、共分散 は 正 になります。

逆に考えれば、負の相関を持っているデータは 正×負 か、 負×正 になり、ほとんどが負の数になります。

つまり、共分散の正負は、相関関係の正負の目安になる といっても良さそうです。

次に、相関関係の強弱を考えるために、共分散 $s_xy$ を、 $s_x$ と $s_y$ の積、 $s_x s_y$ で割った値について考えてみます。
この値を、$x$ と $y$ の 相関係数 といい、 $r$ で表現します。

さて、それではRで (height, weight) と、 (height2, weight2) の相関関係を計算してみます。

correlationCoefficient.R

height  = c(181,179,167,190,166,155,162,168,164,170,181)
weight  = c(99, 58, 68, 85, 57, 46, 60, 59, 60, 71, 77)

height2 = c(170,172,174,176,178,180,182,184,186,188,190)
weight2 = c(55, 56, 58, 59, 62, 65, 66, 66, 68, 74, 79)

cor(height, weight) # -> 0.7789796
cor(height2, weight2) # -> 0.9682458

これだけ長い計算式ですが、Rでは cor関数が全て導出してくれます。
(height,weight) の相関係数は 0.7789796　、 (height2, weight2) の相関係数は 0.9682458 だということがわかりました。
では、ここで相関係数の重要な性質について説明します。

相関係数の性質

1. $-1 \leqq r \leqq 1$
2. $r=1$ である時、散布図は右上がりの直線になる
3. $r=-1$ である時、散布図は右下がりの直線になる
4. $r\fallingdotseq0$ である時、直線的な相関関係はない

つまり、(height2, weight2) の相関係数は (height, weight) より1に近いので、より強い正の相関関係があると数値的な根拠を持って言えるようになりました。

まとめ

元々は、ここから既存のコロナ感染者のデータから、何かの傾向を見つけ出すところまで取り組みたかったのですが、その後手が回らなくなってしまい、放置されていたため、この分だけでも公開しようと思い、コロナで一斉休校になっている間に書いていた記事を公開しました。
統計は数学的な知識がたくさん要求されるので、興味はありつつなかなか踏み出せない方々は、高校数学で扱う統計レベルから学んでいかれてもいいかもしれません。
私も、大学に進学したらもっと専門的な内容について学びたいと改めて思いました。