高校生のための統計学習教材(総務省統計局)
の解説です。
上記サイトの第1講をダウンロードし、テキストとしてください。テキスト
にそって統計ソフト「RStudio」の適用例を提示します。必要に応じて「なるほど統計学園」の参考ページを参照してください。
Rについては、橋本
洸哉 様 「Rで統計解析」 を参考にさせていただきました。
[参考]なるほど統計学園(初級)4_グラフの作り方(初級編)ヒストグラム
データの分布はヒストグラムで分かります。
分布の形状を数式に当てはめれば、確率的な予想が可能になります。(モデリングといいます)
どんな数式に当てはめるかはヒストグラムを見て判断し、当てはまっているかどうかを統計的に判断することができます。
例えば、あるバスケットチームの試合得点を集計すると左右対称な山形になり、
\(~~~~~f(x)=\dfrac{1}{\sqrt{2\pi}\sigma}e^{-(x-m)^2/2\sigma^2}~~~~~\)
とモデリングできるかもしれません。(正規分布といいます)
サッカーの試合得点はポアソン分布という山が左に寄った関数でモデリングできるという説もあります。
データをヒストグラムに描くことは、分析の第1歩といってもよいでしょう。
さて、RStudioを使ってデータのヒストグラムを描いてみましょう。 RStudioはフォルダを决め、全ての作業をそのフォルダ内で行います。
Rに装備されているデータセット iris を読み込みます。
irisはアヤメという花で、setosa、versicolor、virginica の3つの種類があります。
3種類のアヤメ50個体(合計150)について、がく片(sepal)の長さ、幅、および花弁(petal)の長さ、幅を測定したデータがIrisです。
列名 | 説明 |
---|---|
Sepal.Length | がく片の長さ [cm] |
Sepal.Width | がく片の幅 [cm] |
Petal.Length | 花弁の長さ [cm] |
Petal.Width | 花弁の幅 [cm] |
Species | 花の種類(setosa、versicolor、virginicaの3種) |
よほど花に興味のある人でなければ、がく片や花弁のイメージが沸かないかもしれません。
アヤメについては、
機械学習と情報技術
様 【Python】scikit-learnでIrisデータセットを手軽に用意する
を参照してください。
str(iris)によってデータの概要を表示します。Speciesは種類を表します。
Sepal.Width(がく片の幅)についてヒストグラムを描いてみましょう。
str(iris)
## 'data.frame': 150 obs. of 5 variables:
## $ Sepal.Length: num 5.1 4.9 4.7 4.6 5 5.4 4.6 5 4.4 4.9 ...
## $ Sepal.Width : num 3.5 3 3.2 3.1 3.6 3.9 3.4 3.4 2.9 3.1 ...
## $ Petal.Length: num 1.4 1.4 1.3 1.5 1.4 1.7 1.4 1.5 1.4 1.5 ...
## $ Petal.Width : num 0.2 0.2 0.2 0.2 0.2 0.4 0.3 0.2 0.2 0.1 ...
## $ Species : Factor w/ 3 levels "setosa","versicolor",..: 1 1 1 1 1 1 1 1 1 1 ...
hist(iris$Sepal.Width)
breaks=seq(最小値,最大値,階級数の幅) の設定が可能です。right=FALSE は未満の意味です。
# 最小値,最大値,階級数の幅)の設定
hist(iris$Sepal.Width,
right=FALSE,
breaks=seq( 1,5,0.5 ))
x軸、y軸の範囲、ラベル、タイトル名をつけることもできます。
hist(iris$Sepal.Width,
right=FALSE,
breaks=seq( 1,5,0.5 ),
xlim=c( 0, 5 ),ylim=c( 0, 70),
xlab="Sepal.Width",ylab="frequency",main="histgram" )
グラフの棒をビンといいます。(テキストp1)
ビンの数の決め方には、スタージェスの方法、スコットの方法、平方根選択がありますが、(テキストp3)
ビンの数によって
データの見え方が異なることがあるので注意が必要です。(テキストp6)
高校生のための統計学習教材(総務省統計局)
ではPythonのコードも記述されています。
RStudioとPython両方に親しんでおくことは重要です。
RStudioの中でPythonが利用できます。
Pythonのコードとグラフを示しますので参考にしてください。
iris = sns.load_dataset('iris')
print(iris.head())
## sepal_length sepal_width petal_length petal_width species
## 0 5.1 3.5 1.4 0.2 setosa
## 1 4.9 3.0 1.4 0.2 setosa
## 2 4.7 3.2 1.3 0.2 setosa
## 3 4.6 3.1 1.5 0.2 setosa
## 4 5.0 3.6 1.4 0.2 setosa
sns.histplot(iris["sepal_width"], alpha=0.5,kde=True)#kdeは近似曲線の表示
この近似曲線を数式で表すことにより、推測に役立てようというのがデータサイエンスです。
[参考]なるほど統計学園(初級)4_グラフの作り方(初級編)箱ひげ図
箱ひげ図は、ヒストグラムを単純化した分布図です。
irisの”setosa”を抽出したデータをdataとし、そのSepal.Widthの箱ひげ図を描くと、
data <- iris[iris$Specie=="setosa",]
boxplot(data$Sepal.Width)
○は、外れ値を表します。
箱ひげ図は、複数の群の分布を並べて比較するときに便利です。
irisのSepal.Widthを種類別に比較してみましょう。
boxplot(Sepal.Width~Species, iris, ylab="Sepal.Width") #x軸には変数Speciesを指定、さらに引数ylabでy軸のラベルを指定
ヒストグラムの広がりを数値で表したものが分散です。平均との差の2乗の平均で計算されます。
データが、\(X=x_1,x_2,\cdots ,x_n\) であれば、分散V(X)は、平均をmとして
\(V(X)=\dfrac{1}{n} \left[(x_1-m)^2+(x_2-m)^2+\cdots+(x_n-m)^2\right]\)
例えばデータが \(X=3,3,3\)
であれば、平均 \(m=3\) 、散らばっていませんから
\(V(X)=\dfrac{1}{3}
\left[(3-3)^2+(3-3)^2+(3-3)^2\right]=0\)
\(Y=2,3,4\) であれば、平均 \(m=3\)で
\(V(Y)=\dfrac{1}{3}
\left[(2-3)^2+(3-3)^2+(4-3)^2\right]=\dfrac{2}{3}\)
\(Z=6,7,8\) であれば、平均 \(m=7\)ですが、散らばりは同じで
\(V(Z)=\dfrac{1}{3}
\left[(2-3)^2+(3-3)^2+(4-3)^2\right]=\dfrac{2}{3}\)
計算式を簡潔に表すには\(\sum\)
(シグマ)という 記号が必要です。なるほど統計学園の説明も\(\sum\)を用いています。
\(\displaystyle \sum_1^n\)
は、1番目からn番目までたしあわせるという意味です。
この記号を用いれば、
\(\displaystyle V(X)=\dfrac{1}{n}
\sum_1^n(x_i-m)^2\) 簡潔に表現することができます。
統計の世界では、nの代わりに(n-1)で割って
\(\displaystyle V(X)=\dfrac{1}{n-1}
\sum_1^n(x_i-m)^2\)
を分散として用います。(不偏分散といいます)
理由は、実際のデータ分析では全てのデータを入手することは不可能で分散が小さめに出てしまいます。
そこで、(n-1)で割り分散を大きめに見積もることで、より正確な分散の推定値としています。
下に示したRStudioでの分散計算も、(n-1)で割る不偏分散が用いられています。
A=c(3,3,3)
str(A)#データAの内容
## num [1:3] 3 3 3
var(A)#データAの分散
## [1] 0
hist(A,right=FALSE,breaks=seq( 2,5,1 ))
B=c(1,2,3)
str(B)#データBの内容
## num [1:3] 1 2 3
var(B)#データBの分散
## [1] 1
hist(B,right=FALSE,breaks=seq( 1,4,1 ))
分散は平均との差の平方の平均で、平方しているため値が大きくなりすぎます。
このため、分散の平方根である標準偏差が実際のイメージに合うでしょう。
実際のデータで調べてみると、
平均を中心にして
(平均)±
3(標準偏差) の範囲にヒストグラムの山はおさまります。
平均が 3 標準偏差が 0.5 のとき、ガウスの正規分布
\(f\left( x\right) =\dfrac{1}{\sqrt{2\pi
}\cdot 0.5}\exp \left[
-\dfrac{1}{2}\left(\dfrac{x-3}{0.5}\right)^2\right]\)
をグラフで描画すると、
curve(dnorm(x, 3, 0.5),1.5,4.5)
平均 3 を中心に、標準偏差0.5の3倍の幅にグラフがあることが分かります。
\(\dfrac{x-3}{0.5}\) という式に注意を向けた人もいるでしょう。
平均0,標準偏差1の正規分布を標準正規分布といい、
\(\phi \left( z\right)
=\dfrac{1}{\sqrt{2\pi }}\exp \left(
-\dfrac{z^{2}}{2}\right)\) で表されます。
\(z=\dfrac{x-3}{0.5}\)
はグラフをx軸方向に0.5倍に縮小し、x軸方向に3平行移動する変換です。
(標準化といいます)
curve(dnorm(x, 0, 1),xlim=c(-2, 4.5),ylim=c(0, 0.8), col=3)
curve(dnorm(x, 3, 0.5),xlim=c(1.5, 4.5), add=T)
緑のグラフが標準正規分布で、x軸方向に0.5倍に縮小しx軸方向に3平行移動することにより、黒のグラフに移っていることが分かります。
irisの”virginica”について、横軸にPetal.Length、縦軸にSepal.Lengthをとってグラフに点を打ってみましょう。
data <- iris[iris$Specie=="virginica",]
plot(Sepal.Length ~ Petal.Length, data) #がく片の長さをyに、花びらの長さをxに指定
正比例とまではいきませんが、花びら(Petal)の長いものは、がく片(Sepal)も長い傾向にあるようです。
散布図の直線度合いを表す指標として考えられたのが、相関係数です。
\(r_{xy} =
\dfrac{\frac{1}{n}\displaystyle\sum_{i=1}^{n}(x_i-\overline{x})(y_i-\overline{y})}{\sqrt{\frac{1}{n}\displaystyle\sum_{i=1}^{n}(x_i-\overline{x})^2}\sqrt{\frac{1}{n}\displaystyle\sum_{i=1}^{n}(y_i-\overline{y})^2}}\)
分母はxの標準偏差とyの標準偏差の積になっており、分子は共分散と言います。
実際の計算では次の式を使うことができます。
\(r_{xy} =
\frac{\displaystyle\sum_{i=1}^{n}(x_i-\overline{x})(y_i-\overline{y})}{\sqrt{\displaystyle\sum_{i=1}^{n}(x_i-\overline{x})^2}\sqrt{\displaystyle\sum_{i=1}^{n}(y_i-\overline{y})^2}}\)
ここでは、計算に立ち入ることはやめておきましょう。
右上がりの直線状態が強ければ1に近く、右下がりの直線状態が強ければ−1に近い。直線状態がなければ0に近いとイメージしておけばよいでしょう。
なるほど統計学園(上級)「複数の変数の関係性を見る」のサイトでは、
散布図に当てはまる近似直線を求める方法が説明されています。
この近似直線を回帰直線といいますが、x,y
を標準化したときの傾きは相関係数になります。
\(\dfrac{y-\overline{y}}{s_y}=r_{xy}\left(\dfrac{x-\overline{x}}{s_x}\right)\)
この式を整理すれば回帰直線が求まりますが、ここでは立ち入らなくてよいでしょう。
参考までに回帰直線の傾きは
\(\displaystyle
r_{xy}\cdot\dfrac{s_y}{s_x}=\dfrac{s_{xy}}{s_xs_y}\cdot\dfrac{s_y}{s_x}=\dfrac{s_{xy}}{s_x^2}=\frac{\sum
^{n}_{i=1}\left( x_{i}-\overline{x}\right) \left(
y_{i}-\overline{y}\right) }{ \sum ^{n}_{i=1}\left(
x_{i}-\overline{x}\right) ^{2}}\)
で計算されます。
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## `geom_smooth()` using formula = 'y ~ x'
回帰直線の当てはまりを表す指標は決定係数\(R^2\)
で、単回帰分析では相関係数の2乗で表されます。
\(0<R^2<1\)
直線にピタリと当てはまっていれば、相関係数=1、決定係数=1
当てはまりが悪ければ、決定係数は1から離れ0に近づきます。