このページは、 テキストを参考にして完成するワークブックです。
サイトも参考に復習することをお勧めします。

1 これまでに学んだこと.

  • ヒストグラム:データの頻度をグラフ化したもの.
    (面積が確率を表すようにすると確率密度関数になる)
  • 累積相対度数折線:累積度数を確率グラフにしたもの.
    (縦軸が累積確率を表す).

このセクションでは代表値を学びますが、

   ヒストグラムと累積度数折線における代表値の意味.

を理解することが重要です。

2 代表値.

2.1 最頻値.

  • ヒストグラムの山の最も高くなるところ.
  • 確率密度関数のグラフが最大となる x.

確率的に最も起こりやすいデータ x であり、推測統計では最尤推定値と呼ばれ重要である。
最頻値はグラフを見ればすぐ分かる。

2.2 中央値.

メジアン(median)と呼ばれ、大小の順に並べて中央に位置するxである。
データが偶数個で中央のデータがない場合は前後の値の平均を中央値とする。
(グラフ上では平均は2つの点の中点であり直線近似で求めていることになる).

  • 確率密度関数の面積を左右二等分するx.
  • 累積相対度数折線で確率が0.5となるx.
    (累積度数が全個数の半分となる x と言ってもよい).

実データがあれば中央値は容易に求まるが、度数分布表のみから算出する計算は手間がかかり(p13)誤差もあるので深入りする必要はない。

2.3 平均.

まず p8 の例で、50人の平均点 3.12 を計算する方法を見てみましょう。

得点 人数
1 4
2 6
3 25
4 10
5 5

平均は50個を合計して50で割ります。
同じ値が繰り返されているので掛け算を用いて合計すると.
   \(m=\dfrac{1}{50}\left[1\cdot4+2\cdot6+3\cdot25+4\cdot10+5\cdot5\right]=3.12\).

実際の計算では、次のような表を作ります。

得点 人数 得点✖️︎人数
1 4 4
2 6 12
3 25 75
4 10 40
5 5 25

ここで分配法則を用いて、
   \(m=1\cdot\dfrac{4}{50}+ 2\cdot\dfrac{6}{50}+3\cdot\dfrac{25}{50}+4\cdot\dfrac{10}{50}+5\cdot\dfrac{5}{50}\)  とすると.

平均は、値と確率をかけてたしたものと見ることができます。

度数を\(f_i\) 総数を\(N\) とすると、
   \(m=x_1\cdot\dfrac{f_1}{N}+x_2\cdot \dfrac{f_2}{N}+\cdots+x_n\cdot\dfrac{f_n}{N}\).
p9の参考では、
   \(m=\dfrac{1}{N}\left[x_1\cdot f_1+x_1\cdot f_1+\cdots+x_n\cdot f_n\right]\).
となっていますが同じ式です。

平均は期待値と呼ばれます。
期待値は、度数を重みとした加重平均になっています。

注意点として、

  • てこの原理で、重心から離れたデータの影響が大きい.
    (重心から離れたデータを外れ値という).
  • 中央値は外れ値の影響を受けないので、代表値としてよく用いられる.

この2点は覚えておきましょう。

この例での最頻値は 3 、中央値も 3 となっています。

<力学的意味>.
(確率密度を重さと解釈したとき)
平均はグラフと横軸で挟まれる面積の重心を意味します。
数理的・力学的理由は以下の通りですが、立ち入らなくてよいです。
わすれなぐさ 様 によれば、物理における重心の定義は

   \(\displaystyle x_G = \frac 1M \sum x_i m_i= \sum x_i\cdot \frac {m_i}{M}\)

   \(\dfrac {m_i}{M}\) を確率と解釈すれば重心は期待値の式に一致します。

3 問題.

3.1 例題2.

data <- data.frame(x=seq(22,30,by=2),y=c(1,2,10,5,2))
colnames(data) <-c("階級値","度数")
data
階級値 度数
22 1
24 2
26 10
28 5
30 2

度数分布表

library(stringr)
dat <- data.frame(階級=str_c(seq(22,30,by=2)-1,"-",seq(22,30,by=2)+1),度数=c(1,2,10,5,2))
dat
階級 度数
21-23 1
23-25 2
25-27 10
27-29 5
29-31 2
library(ggplot2)
ggplot(data, aes(x=階級値, y=度数),color = 3) + geom_bar(stat="identity")+theme_gray(base_family = "HiraKakuPro-W3") + coord_fixed((ratio=1/2)) +
  scale_x_continuous(breaks=seq(0,30,2))+ scale_y_continuous(breaks=seq(21,31,2))+ geom_col(color = "blue",fill="grey92")

最頻値が26であることがわかる。

3.1.1 累積度数.

累積度数を求める関数は cumsum です。

data <- data.frame(x=seq(22,30,by=2),y=cumsum(c(1,2,10,5,2)))
colnames(data) <-c("階級値","累積度数")
data
階級値 累積度数
22 1
24 3
26 13
28 18
30 20
library(ggplot2)
ggplot(data, aes(x=階級値, y=累積度数),color = 3) + geom_bar(stat="identity")+theme_gray(base_family = "HiraKakuPro-W3") + coord_fixed((ratio=1/2)) +
  scale_x_continuous(breaks=seq(21,31,2))+ ylim(c(0, 20))+ geom_col(color = "lightblue",fill="grey92")+geom_line(color = "blue")+ geom_line(data = data, aes(x=階級値+1, y=累積度数),col="red")+geom_hline(yintercept = 10, size = 0.9, linetype = 3, color = "red")+geom_vline(xintercept = 26.4, size = 0.9, linetype = 3, color = "red")
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

この折線だと30で累積20になるので、31で累積20になるよう修正する。
棒グラフの中央を繋いでいるグラフを、右の端点を繋ぐように修正する。
科学センスを目指して 様「ヒストグラムの中央値の求め方」 によれば、
  「横軸に各区間の最大値、縦軸に累積度数を取ってグラフを描きます」

累積相対度数折線を用いて中央値を求めよう

横軸が各区間の最大値となるよう x を1増やします。

data1 <- data.frame(x=c(23,25,27,29,31),y=cumsum(c(1,2,10,5,2)))
plot(data1,type="l",xlim=c(23,31))
axis(side=1, at=seq(21,31)) 
abline(v = 26.4,col="lightblue")
abline(h = 10,col="lightblue")

(25,3), (27,13)を結ぶ直線とy=10 の交点が中央値

\(y-3=\dfrac{13-3}{27-25}(x-25)\) に y=10を代入して
   \(10-3=\dfrac{13-3}{27-25}(x-25)\).
    \(x-25=\dfrac{27-25}{13-3}(10-3)=1.4~~~~\therefore ~~x=26.4\).

累積度数がちょうど半分になるところが中央値であり、26.4 であることが分かります。
中央値のある区間の累積関数を直線で近似して求めています。
テキストp13に計算法がありますが、直線の方程式の応用です。

簡便法として、「中央値の所属する階級の階級値を中央値とする」立場もあるとテキストp13の注に記されています。
もの立場では、中央値は 26 になります。

x=seq(22,30,by=2)
X <- rep(x,c(1,2,10,5,2))
X
##  [1] 22 24 24 26 26 26 26 26 26 26 26 26 26 28 28 28 28 28 30 30
median(X)
## [1] 26

度数分布表から求める中央値は直線近似を用いた推定値であり、誤差があります。元になるデータが分からないのですから、あまり気にする必要はありません。

3.2 演習2.

data <- data.frame(x=seq(10,60,by=10),y=c(10,3,2,4,6,15))
colnames(data) <-c("階級値","度数")
data
階級値 度数
10 10
20 3
30 2
40 4
50 6
60 15

度数分布表

library(stringr)
dat <- data.frame(階級=str_c(seq(10,60,by=10)-5,"-",seq(10,60,by=10)+5),度数=c(10,3,2,4,6,15))
dat
階級 度数
5-15 10
15-25 3
25-35 2
35-45 4
45-55 6
55-65 15
library(ggplot2)
ggplot(data, aes(x=階級値, y=度数),color = 3) + geom_bar(stat="identity")+theme_gray(base_family = "HiraKakuPro-W3") + coord_fixed((ratio=1)) +
  scale_x_continuous(breaks=seq(5,65,5))+ ylim(c(0, 40))+ geom_col(color = "blue",fill="grey92")

最頻値が60であることがわかる。

3.2.1 累積度数.

累積度数を求める関数は cumsum です。

data <- data.frame(x=seq(10,60,by=10),y=cumsum(c(10,3,2,4,6,15)))
colnames(data) <-c("階級値","累積度数")
data
階級値 累積度数
10 10
20 13
30 15
40 19
50 25
60 40
library(ggplot2)
ggplot(data, aes(x=階級値, y=累積度数),color = 3) + geom_bar(stat="identity")+theme_gray(base_family = "HiraKakuPro-W3") + coord_fixed((ratio=1)) +
  scale_x_continuous(breaks=seq(5,65,5))+ ylim(c(0, 40))+ geom_col(color = "lightblue",fill="grey92") + geom_line(color = "blue")+geom_line(color = "blue")+ geom_line(data = data, aes(x=階級値+5, y=累積度数),col="red")+geom_hline(yintercept = 20, size = 0.9, linetype = 3, color = "red")+geom_vline(xintercept = 46.7, size = 0.9, linetype = 3, color = "red")

階級の右端をつないだ折線を作図し、中央値を図から推定しましょう

累積度数折線は次のグラフになり、グラフから中央値が 46.7 と推定されます。

data1 <- data.frame(x=c(15,25,35,45,55,65),y=cumsum(c(10,3,2,4,6,15)))
plot(data1,type="l",xlim=c(15,65),ylim=c(0,40))
axis(side=1, at=seq(10,65,by=5)) 
abline(v = 46.7,col="lightblue")
abline(h = 20,col="lightblue")

4 研究.

ヒストグラムの面積を2等分する線として、中央値を求めてみましょう。

4.1 演習2.

data <- data.frame(x=seq(10,60,by=10),y=c(10,3,2,4,6,15))
colnames(data) <-c("階級値","度数")

ggplot(data, aes(x=階級値, y=度数),color = 3) + geom_bar(stat="identity")+theme_gray(base_family = "HiraKakuPro-W3") + coord_fixed((ratio=1)) +
  scale_x_continuous(breaks=seq(5,65,5))+ ylim(c(0, 40))+ geom_col(color = "grey85",fill="lightblue")+  geom_vline(xintercept = 45+10/6,
             linetype = "dotted",
             color = "red",
             size = 0.5)

45 と 55 を 1:5 に内分する点は.
   \(x=45+\dfrac{1}{6}\times 10=46.66\).

4.2 連続型確率分布へ.

数理モデルでは、グラフの面積が確率を表すようにすることで、データ数の異なるデータを統一的に扱います。
代表値の図形的意味を理解しておくと、連続型分布の理解が容易になります。

4.2.1 サッカーの得点.

サッカーの1試合での得点は小さい値で、大きな得点になる確率は小さいので、ポアソン分布に従うと言われています。

<参考サイト>.

小中英嗣 様 「ポアソン分布を利用したサッカーシミュレーション」.

平均得点を3.5 として、1000個のデータを出力しヒストグラムを描きます.

x <- rpois(1000, 3.5)
head(x)
## [1] 1 0 2 2 0 6
median(x)
## [1] 3
table(x)
## x
##   0   1   2   3   4   5   6   7   8   9  10  11 
##  35 108 184 216 180 141  75  38  14   7   1   1
names(which.max(table(x)))
## [1] "3"
h <- hist(x,prob=T)
abline(v = median(x),col="red")
abline(v = mean(x),col="green")
abline(v = names(which.max(table(x))),col="blue")

ヒストグラムに偏りがあるので.   最頻値 = 中央値 < 平均値  となっています。

4.3 数理統計モデル.

4.3.1 ガンマ分布.

Covid19においては数理統計モデルによる解析が注目されました。
国立感染症研究所によると、
 発症間隔とは一次感染者の発症時刻から二次感染者の発症時刻までの時間間隔.
ですが、 新型コロナウイルス感染症対策専門家会議(第 17 回)資料1では、

「発病間隔データの確率密度関数を基に、平均 4.8 日、標準偏差 2.3 日のガンマ分布に従うものと想定した」.
とあります。

この「ガンマ分布、平均、標準偏差」という3つの情報が、Covid19の発症間隔についての詳しい情報を伝えてくれます。

@kidaufo 様「ガンマ分布のパラメータ指定方法まとめ」 によると、

\(\theta = \frac{v}{\mu}\)  \(k = \frac{\mu^2}{v}\)  より、

  • shape parameter \(k=4.35\).
  • scale parameter  \(\theta=1.10\).
x <- rgamma(1000, shape=4.35, scale = 1.1)
median(x)
## [1] 4.553499
mean(x)
## [1] 4.907656
hist(x,prob=T,breaks="sturges")
abline(v = median(x),col="red")
abline(v = mean(x),col="green")
curve(dgamma(x, shape=4.35, scale = 1.1),add=T,col="blue")

最頻値はグラフが最大となる x です。

  最頻値 < 中央値 < 平均    が読み取れます。