はじめに、そしてRについて

 この講義で使用する教科書は、東大出版会の統計学入門です。万人にオススメできる標準的な教科書になってます。

 統計では理論もデータ解析能力(i.e., プログラミング能力)もどちらも重要です。講義中ではプログラミングについて解説する時間はなさそうなのですが、Rによるコードも紹介しようと思います。

 Rはフリーで利用できる統計解析ソフトです。現在はPythonの陰に隠れてあまり目立ちはしませんが、統計処理や結果の図示には好まれて利用されることが多いようです。特に、最近ではRStudioと呼ばわれる無料で利用できる統合環境が整ってきており、とても使い勝手が良いと思います。

 この講義資料は、RStudioを利用して作っています。RStudioでは文章作成にMarkdown文法が利用可能です。 ちなみにRStudioではPythonも使えます。英語に限定すれば(おそらく)、LaTeXも利用できます。 一度は触れてみることをおすすめします。

 例えば、同じ作業ディレクトリにある.csvファイルにはこのようにアクセスできます。 参考: Rで学ぶ確率統計学 一変量統計編のwebページ

・関数read_excelが含まれているreadxlのようなR packageは無料でダウンロード可能です。起動するたびにパッケージは読み込む必要があります。

# setwd("xxx") # ここに自分の作業ディレクトリを打ち込みます。
library(readxl) # R特有。便利なライブラリを読み込む。ここでは.xlsxの読み込み関数を含むパッケージをインストール
## Warning: package 'readxl' was built under R version 3.5.2
sampledata <- read_excel("sampledata.xlsx") # データ読み込み (データの代入は矢印でも=でもOKです)
math = sampledata$math  # データ読み込み
head(math, 3) # データ最初の3行表示
## [1] 90 65 60
summary(math) # データの代表的な統計量を表示
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   10.00   45.00   65.00   61.61   80.00  100.00
length(math) # データ長
## [1] 146
math[14] # 14番目のデータ
## [1] 80

データの可視化

 データを見てみましょう。「データ」とは、例えばクラス全体の数学の得点、身長、体重、など様々に計測したものを意味します。例えば、とあるクラスの数学のテストの得点を見てみましょう。

##   [1]  90  65  60  75  50  60  25  35  85 100  45  95  80  80  95  55  55  65
##  [19]  50  65  40  85  70  65  65  15  65  75  60  75  25  65  55  40  70  60
##  [37]  60  35  75  45  60  55  75  75  90  35  20  90 100  65  15  80  80  85
##  [55]  25  10  95  45  55  80  80  40  90  95  90  70  80  25  65  30  60  25
##  [73]  35  20  70 100  50  60  50  90  65  45  35  50  85  30  50  70  85  85
##  [91]  75  80  40  60  70  80  95  50  40  25  60  35  30  45  35  70  75  55
## [109]  40 100  50  55  90  95  60  70  80  70  45  55  75  45  65  30  25  85
## [127]  90  90  60  85  20  75  70  75  50  90  60  15  90  75  70  35  95  60
## [145]  45  55

さて、このデータを見るだけで何かわかったことはあるでしょうか？データの羅列を見ても、その特徴はほとんど何もわからない(気がする)。

  重要なことは、データを見ることです。まず見るべきは、可能であればヒストグラム。ヒストグラムとは、とある範囲内のデータがいくつ観測されたかを示すグラフのことです。観測数をそのまま表示しても、観測回数で割り算したものでもどちらでもOKです。

hist(math) # 度数

hist(math, prob=TRUE) # 相対度数

 なんとなく数字の羅列を見るよりはわかった気になるかと思います。51-60点の人が一番多い(農工大では悲惨な結果になってしまう得点レンジですね…)、評点で言うところのSが10人くらい、Aが20人くらい…など、なんとなく数字見るよりはわかりやすいかと思います。みなさんがデータを解析することになるとすると、実験のときや卒論のときだと思うのですが、とりあえずヒストグラム書いてみるというのは悪くない手段です。ただし、ヒストグラムがわかりやすい状況は見たいデータが1変数ないしは2変数のときに限定されます。

 また、Rは統計解析に加えて、グラフィクスにも定評があります。 特に、ggplotを利用した図表は可視性が高いという評判 (ただ文法がややこしいと思う人が多いらしい)。

library(ggplot2)
## Warning: package 'ggplot2' was built under R version 3.5.2
qplot(math, geom="bar", xlab="score", ylab="count number")

データの要約(平均)

 データを見た後には、データの特徴を捉える指標を計算しましょう。平均値\(\bar{x}\)は、データの数が\(N\)、\(i\)番目のテストの得点が\(x_i\)のとき(\(i = 1, ..., N\))、

\[\begin{equation} \bar{x} = \frac{1}{N}\sum_{i=1}^N x_i \end{equation}\]

により与えられます(\(N\)はデータ長)。Rでは平均の計算など簡単にこなすことが可能です(Rでなくとも大方のソフトウェアでも簡単)。

xbar = mean(math) # 平均値
xbar
## [1] 61.60959
# 平均値の計算の仕方その2
xbar = sum(math)/length(math)
xbar
## [1] 61.60959
## 愚直に平均を計算する ##
xbar = 0
N = length(math)
for(i in 1:N){
  xbar = xbar + math[i]
}
xbar = xbar/N
xbar
## [1] 61.60959

データの要約(分散)

 平均値のみでは、例えば全員が61.60959点のときと、さきほど解析したような場合の区別がつきません。 全員が50点のテストか、0点から100点の人がいて平均が50点のテストとでは違うものでしょう。 そこで、データのばらつきも特徴の一つになりえそうです。

 データのばらつきを表現する際には、分散\(\mathrm{Var}(x)\)、または標準偏差\(\mathrm{SD}(x)\)を利用することが多いです。これらは、

\[\begin{equation} \mathrm{Var}(x) = \frac{1}{N}\sum_{i=1}^N(x_i - \bar{x})^2 \end{equation}\]

\[\begin{equation} \mathrm{SD}(x) = \sqrt{\mathrm{Var}(x)} \end{equation}\]

により与えられます。平均と単位が揃う指標は、分散ではなく標準偏差であることに注目してください。つまり、平均と足し算引き算できるのは、標準偏差です。

Var = var(math) # 分散
Var
## [1] 506.5293
# 分散の計算の仕方その2
Var = sum((math-xbar)**2)/(length(math)-1)
Var
## [1] 506.5293

 ちなみに、ほぼすべてのソフトウェア上ではなぜか、分散は\(N\)ではなく、\(N-1\)で割り算されています。 その理由は意外と面白くそしてやや複雑であるため後々説明しますが、\(N-1\)で割らないと真の分散値に近くならないという性質があります。

 分散は、下記のようにばらつきを意味します。

library(ggplot2)
data_plot = data.frame(
  x1 = rnorm(1000,0,1), # 平均0、分散1(標準偏差1)の正規乱数(ガウス分布という分布から生成した確率変数)
  x2 = rnorm(1000,0,4)  # 平均0、分散16(標準偏差4)の正規乱数(ガウス分布という分布から生成した確率変数)
)
head(data_plot)
##           x1         x2
## 1  1.6217864 -0.7477662
## 2  0.9213522 -4.0740587
## 3  1.1787546 -4.6402474
## 4  2.0293072  0.7341869
## 5 -0.3649124 -0.7147746
## 6 -0.3818758  5.7979682
k = ceiling(1 + log2(1000)) # スタージェスの公式(ヒストグラムのビンの数を決める方法、使ってない)

# 画像として保存
ggplot(data_plot, aes(x=x1)) + geom_histogram(color="blue", bins = k, position="identity", alpha=0.0, fill = "white", binwidth = 0.2) + geom_histogram(aes(x=x2), color="red", bins=k, position="identity", alpha=0.0, fill = "white", binwidth = 0.2)

赤枠のデータの方が分散が大きく、ばらついていることがわかる。

分位点と箱ひげ図

 他にもデータを要約する方法は様々あります。例えば中央値も有力な指標です。 (2k+1)個(kは正の整数)のデータを観測したとき、小さい順ないしは大きい順に(k+1)番目のデータを中央値といいます。 中央値は外れ値と呼ばれる明らかに変な値があってもあまり変化しない頑健性を有します。それに対して、平均、分散は外れ値には弱い性質をもちます。 中央値のわずか欠点と言っても差し支えない点は、解析計算がわずかにややこしいことです。平均値、分散の方が解析しやすいため、好まれることが多いように思います。

## [1] 65
## [1] 78

データが偶数個のときは、k番目、(k+1)番目のデータの平均が中央値。

 その他の指標を挙げます。小さい順に、25%の点を第一四分位点(Q1)、75%の点を第三四分位点(Q3)と呼びます。 Q3-Q1は、四分位偏差と呼ばれる(inter-quantile range, IQR)。 これらを示した図を箱ひげ図と呼びます。

中央の黒い太線が中央値、中央の近くの四角の上線がQ3, 下線がQ1、線の一番上は最大値、一番下は最小値。 このように、箱ひげ図とデータを同時にプロットすると、色々わかりやすい。

演習課題

1: データの数が\(N = 100\)、\(i\)番目のデータの値が\(x_i = i\)のとき(\(i = 1, ..., 100\)) [例えば5番目のデータの値は5、12番目のデータの値は12ということです]、下記の問に答えよ。なるべく手書き(LaTeXやMarkdown使ってもらっても構わない)で計算せよ。そして、できる人はR、RStudio、Julia、もしくはPythonなどお好みのソフトウェアにより、その手書きの計算が正しいことを確認せよ。


2: 手元にある

\(\boldsymbol{x}_{\rm orig} = (1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11)\)

という11個のデータがある。ある日子供がPCをいじってしまい、

\(\boldsymbol{x} = (1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 100, 100)\)

と値を追加してしまった。このとき、\(\boldsymbol{x}_{\rm orig}, \boldsymbol{x}\)の分散と中央値を計算せよ。どちらが大きく変化したであろうか。


3: データの数が\(N\)(自然数)、\(i\)番目のデータの値が\(x_i\)(\(i = 1, ..., N\))、その平均を\(\bar{x}\)、分散を\(\mathrm{Var}(x)\)とするとき、以下の等式を証明せよ。

\[\begin{equation} \mathrm{Var}(x) = \frac{1}{N}\sum_{i=1}^Nx_i^2 - \bar{x}^2 \end{equation}\]

分散は2乗の平均から平均の2乗を引いたものになる。これが必ず0より大きくなることも確認せよ。