正規分布X,Yの和X+Yについてjulia言語で考えてみた。
はじめに
Twitterから
この,会話文の中で,
とありました。「まあ,違うよな...」と思いながら,「どんな分布になるんだろう?」と思ったので,julia言語
を使って分布を書いてみることにしました。
日本人の身長の分布について
問題文で,太郎さんと花子さんの身長と,日本人の男女別の身長の平均と標準偏差が与えられています。そして,身長の分布は正規分布に従うそうです。
- 太郎くんの身長 168.6cm
- 花子さんの身長 185.4cm
太郎くんは私と同じくらいです。花子さんは高いですね。
平均 | 標準偏差 | |
---|---|---|
日本人の男子 | 171.2cm | 5.8cm |
日本人の女子 | 158.2cm | 5.2cm |
まずは問題を解いてみます。
julia言語
を用いて,問題を解いてみます。
第1問
ポストでは正規分布表も与えられているのですが,ここはコマンドで求めてみます。Distributions.jl
のパッケージを利用します。コマンドccdf(f,x)
は分布fでx以上の確率を返します。
using Distributions # 分布パッケージ
ccdf(Normal(158.2,5.2),168.6) # ccdf(f,x) f の分布で x 以上の確率を返す。
0.02275013194817915
太郎さんより身長が高い日本人女子は2.27%。100人中2人くらいしかいませんね。
ccdf(Normal(171.2,5.8),185.4)
0.007177085392942717
花子さんより身長が高い日本人男子は0.71%。100人中1人もいません。
次のような結論になります。
日本人女子と日本人男子が同数とすると,太郎さんより身長の高い日本人女子の方が,花子さんより身長の高い日本人男子よりも多い。
第2問
このことを考える前に,正規分布について成り立つ性質があります。
この性質は正規分布の再生性と呼ばれ,多くのサイトで説明があります。
「今回はこの性質なのかな?」と考えてみると,ちょっと違うことがわかります。まず,集団Aの平均は集団BとCの平均をそれぞれ加えたものではありません。おそらく男女同数であれば,Aの平均は2で割る必要があると思われます。
今回の分布は2つの正規分布を1:1で混ぜた混合ガウス分布(GMM) となります。確率密度関数それぞれ,0.5倍して加えます。
これらのグラフをStatsPlots.jl
パッケージを使って書いてみます。
using StatsPlots
X = Normal(171.2,5.8) # 日本人男子の分布(正規分布)
Y = Normal(158.2,5.2) # 日本人女子の分布(正規分布)
Z = Normal(158.2+171.2,√(5.8^2+5.2^2)) # X+Yの分布
G = MixtureModel([X,Y],[.5,.5])# 混合ガウス分布 X:Y = 1:1
plot(X,label="JPN men")
plot!(Y,label="JPN women")
plot!(Z,label="X+Y")
plot!(x->pdf(G,x),label="GMM")
第3問
混合ガウス分布を利用して,求めてみます。
ccdf(G,185.4)*100
0.3588584898003955
花子さんは上位0.358%であることがわかります。
これは男女別に調べてもわかります。
a = ccdf(X,185.4)*100
0.7177085392942717
b = ccdf(Y,185.4)*100
8.440306529671548e-6
(a+b)/2
0.3588584898004007
同じ結果ですね。
正規分布表を用いて手計算で解いてみる。
第1問
太郎さんより身長が高い日本人女子は2.28%。
花子さんより身長が高い日本人男子は0.71%。
第2問
一般的に正しくない。分布は一般に混合ガウス分布となり,正規分布とは異なる。
第3問
正規分布表は3.09までしかないですね。
よって,
花子さんは上位0.355%であることがわかります。
データで確認してみる。
- 男子1,000,000人,女子1,000,000人を正規乱数を用いてヒストグラムを作成。
- 正規乱数で作成した男子1,000,000人,女子1,000,000人を混ぜて,GMMのヒストグラムを作成。
- 正規乱数で作成した男子1000人(X),女子1000人(Y)より,X+Y(データ数1,000,000)のヒストグラムを作成。
using StatsPlots,Distributions
function graph(n)
p = map(x->(5.8*x+171.2),randn(10^(2*n))) # 日本男子10^n人
q = map(x->(5.2*x+158.2),randn(10^(2*n))) # 日本女子10^n人
l=[] # GMM
for i=1:Int(10^(2*n))
append!(l,p[i])
append!(l,q[i])
end
k=[] # X+Y
for i=1:10^n,j=1:10^n
append!(k,p[i]+q[j])
end
# ヒストグラム
histogram(p,alpha=0.1,label="mem")
histogram!(q,alpha=0.1,label="womem")
histogram!(l,alpha=1,label="GMM")
histogram!(k,alpha=0.1,label="X+Y")
end
graph(3)
X+Yの分布を除いてみてみます。
やはり,正規分布ではないですね。
Discussion