身長や体重の分布は正規分布か - ジョンとヨーコのイマジン日記

はじめに：身長の分布

身長の分布が正規分布っぽくなるというのはよく言われる。

ためしに、

学校保健統計調査－令和2年度（確定値）の結果の概要：文部科学省

から高校生男子17歳の身長の分布を持ってきて図示してみる。

保険統計調査のエクセルの表は1cm刻みで割合（千分率）になっている。

青い点線が正規分布で黒い実線は身長の分布である。

なかなか正規分布っぽい。

この青い点線の正規分布のパラメータは度数分布表から平均と標準偏差を定義通りに求め、

平均： $\bar x = \sum_i x_i p_i$

標準偏差： $s = \sum_i (x_i - \bar x)^2 p_i$

とした。ここでは階級 $x_i$ の確率（割合）を $p_i$ と表記した。

体重の分布

同じように、高校生男子17歳の体重の分布を持ってきて図示してみる。

正規分布らしくはまったくない。

通っぽい人なら仮想的なモデル人体みたいなものを考え、その身長が $s$ 倍されるとき、体積はおおむね $s^3$ に比例、体重はおおむね体積に比例する……つまり身長が正規分布なら体重は正規分布に従う確率変数を3乗した分布になるのでは！？と思うかもしれない（2乗3乗の法則 - Wikipedia）。

正規分布に従う確率変数を3乗したときの分布は求めようと思えば求められるので微分積分学などをやった人は計算してみてもいいかもしれないが、ここでは乱数を使って楽をしよう。

階級の値 $x_i$ を $p_i$ に比例する確率で10000回復元抽出して、それを1/3乗してみる。

3乗を1/3乗して、これで正規分布っぽくなる……かと思うとならない。

青い点線が正規分布で黒い実線は乱数で再現した体重の1/3乗の分布である。

たぶん人間の体は大きい人と小さい人で相似みたいにはなってないのかもしれない。

データ分析というのはこういうもので、紙の上だけで考えた理屈はこうやってあっさり退けられる。

身長の分布再考

もう一度身長の分布を見ると最頻値の171cm、170cmあたりで正規分布より少し急峻になっている。

これは偶然だろうか。

この調査の対象になった126900人だそうなので、高校3学年で割って大体40000人くらいのサンプルはいそうである。

先に求めた平均 $\bar x$ 、標準偏差 $s$ の正規分布から40000個乱数を取ってきて、1cm刻みの度数分布表を作ってみよう。

こうすることによって正規分布からのズレが偶然だとして、どれくらいの確率で起こるのかわかる。

こういったやり方をパラメトリック・ブートストラップという。

正規乱数から作った度数分布表を5000回実際の身長の分布に重ねたのが次の図である。

赤い線が実際の身長の分布、乱数から作ったほうがグレーである。

何本もあるグレーの線が赤い線より高くなることはそうそうない。

正規乱数から作った度数分布表の最大値の分布（1-経験分布関数）を次の図に示そう。

縦の点線が実際の身長の分布の171cmのところの確率で、このように値が集中するのは正規分布ではなかなかレアといえそうである。

身長は多分人間が目視で測るので170cmくらいに心理的に誘導されるみたいなことはあり得るかもしれない。

こんな具合に、正規分布に従うデータというのは自然にはまずないと思ったほうがいいと私は思う。

統計学で「正規分布で近似」というときは、多くの場合、データをパラメータ $\theta$ を持つ確率分布 $f(x;\theta)$ でモデル化したときに、 $\theta$ の推定量が正規分布に近づく、という話で、一般のデータ $x$ そのものが正規分布になるみたいなことを保証する定理はない。

この辺の区別はややこしくて、多くの人が当たる壁と思う。

R のコード

まとめて貼るぞ。

library(readxl)
library(ggplot2)
library(dplyr)
#https://www.e-stat.go.jp/stat-search/files?page=1&toukei=00400002&tstat=000001011648

dat_h <-read_excel("./Downloads/r3_hoken_tokei_02.xlsx",skip=4)
h17 <- as.numeric(dat_h$`17歳`)
h17<- ifelse(is.na(h17),0,h17)
len <- length(h17)
height <- as.numeric(gsub("cm|～","",dat_h$区分))

df_height <- data.frame(height=height[-c(1,(len-1):len)],
                        prob=h17[-c(1,(len-1):len)]/1000)

m_h <- with(df_height, sum(prob*height))
s_h <- with(df_height, sqrt(sum(prob*(height-m_h)^2)))

ggplot(df_height, aes(x=height, y=prob))+
  geom_line()+
  geom_line(data = NULL, aes(y=diff(pnorm(c(min(height)-1,height),m_h,s_h))), colour="royalblue", linetype=2, linewidth=1)+
  theme_bw(16)
ggsave("hist_h.png")

######

dat_w <-read_excel("./Downloads/r3_hoken_tokei_03.xlsx",skip=4)
w17 <- as.numeric(dat_w$`17歳`)
w17 <- ifelse(is.na(w17),0,w17)
len <- length(w17)
weight <- as.numeric(gsub("kg|～","",dat_w$区分))
weight <- weight[-1]
df_weight <- data.frame(weight=weight[-1],
                        prob=w17[-1]/1000)
m_w <- with(df_weight, sum(prob*weight))
s_w <- with(df_weight, sqrt(sum(prob*(weight-m_w)^2)))

ggplot(df_weight, aes(x=weight, y=prob))+
  geom_line()+
  geom_line(data = NULL, aes(y=diff(pnorm(c(min(weight)-1,weight),m_w,s_w))), colour="royalblue", linetype=2, linewidth=1)+
  theme_bw(16)
ggsave("hist_w.png")

sim_w <- sample(df_weight$weight, size = 10000, prob=df_weight$prob, replace = TRUE)


ggplot(data.frame(stat=sim_w^(1/3)), aes(stat))+
  geom_density()+
  stat_function(fun = function(x)dnorm(x, mean=mean(sim_w^(1/3)), sd=sd(sim_w^(1/3))), colour="royalblue", linetype=2, linewidth=1)+
  theme_bw(16)+
  labs(x="stat (1/3)",y="density")
ggsave("hist_w3.png")

ggplot(data.frame(stat=sim_w^(1/2)), aes(stat))+
  geom_density()+
  stat_function(fun = function(x)dnorm(x, mean=mean(sim_w^(1/2)), sd=sd(sim_w^(1/2))), colour="royalblue", linetype=2, linewidth=1)+
  theme_bw(16)+
  labs(x="stat (1/2)",y="density")
ggsave("hist_w2.png")

ggplot(data.frame(stat=log(sim_w)), aes(stat))+
  geom_density()+
  stat_function(fun = function(x)dnorm(x, mean=mean(log(sim_w)), sd=sd(log(sim_w))), colour="royalblue", linetype=2, linewidth=1)+
  theme_bw(16)+
  labs(x="stat (log)",y="density")
ggsave("hist_wlog.png")

######
sample_h <- function(i, n, m_h, s_h, height){
  set.seed(i)
  z <- rnorm(n,m_h,s_h)
  r <- range(z)
  b <- seq(floor(r[1]),ceiling(r[2]))
  h0 <- hist(z, breaks = b, plot = FALSE)
  data.frame(height=b[-1],prob=h0$counts/n,id=i)
}

simout <- lapply(1:5000, function(i){sample_h(i,40000, m_h, s_h, height)})
dfout <- bind_rows(simout)

modep <- group_by(dfout,id) %>% 
  summarise(maxp=max(prob))

ggplot(dfout, aes(x=height,y=prob))+
  geom_line(aes(group=id), alpha=0.002)+
  geom_line(data = df_height, colour="firebrick")+
  theme_bw()
ggsave("hist_sample.png")

ggplot(modep, aes(x=maxp))+
  geom_step(aes(y=1-after_stat(y)), stat="ecdf",pad=FALSE)+
  geom_vline(xintercept = max(prob_h), linetype=2)+
  labs(x="stat",y="prob")+
  theme_bw(16)
ggsave("ccdf_h.png")