身長と体重とJulia

histがどういう実装になってるか紐解いていく。マルチメソッドで、同名関数が２つあるけど、引数が少ない方を見る。(Histogram1Dの方ね。同様にHistogram2Dが有るけど、こちらは結果がヒートマップで表示される。何かの時に役立ちそう。)

すると、fitなんてのが核になっている事が分かった。

冒頭で、StatsBaseをusingしてるので、多分これだろうと消去法で当たりを付ける。だって他は、色の扱いっぽいものとREPLなんだもの。

使うに先立って、PkgからStatsBaseを登録して、usingしておく。Gnuplotで使ってるから共用させてよってのは、出来ない相談。地獄を避けるためだ。

julia> h = fit(Histogram, x, 0:100:1500)
Histogram{Int64,1,Tuple{StepRange{Int64,Int64}}}
edges:
  0:100:1500
weights: [15, 42, 20, 6, 3, 2, 1, 0, 1, 0, 1, 0, 2, 1, 0]
closed: left
isdensity: false

ビンゴでしたねぇ。一発で的中したけど、pythonみたいにパッケージのヒントが付いていないので、それらしいのが沢山あると、苦労しそう。

身長5cm単位の分布図 うん、正に正規分布だね。 身長の分布は本当に正規分布に従うのか！？ 統計に関する説明が有った。 身長偏差値チェッカー こんな便利なのもある。みんな統計情報って好きなのね。

上で紹介されてた、統計で見る日本 で、EXECLファイルが落とせるんだけど、おいらそんなの持っていないぞ。 まてよ、 libreofficeで開けるかな？

h26_hoken_toukei_02.xlsx サフィックスの最後にxが付いているのが味噌なんだな。ちゃんと開けて閲覧出来た。けど、欲しいデータに変換するのが超面倒。それに、あのボタン沢山のパネルには、卒倒させられるぞ。よくあんなものを（仕事とは言え）使ってたものだな。

ここは、えいやっと、

こういうデータを使おう。東京農大の保健室が、 身長・体重・BMIの平均値と標準偏差 なんてのを公開してくれていたので、有難く使わせてもらった。

データの統計量(Statistics) って事で解説されてた。 そして、一般式をJuliaに展開されてた。 正規分布

ようするに、この式って確率密度関数　な訳ね。

snd(x,m,s) = (1/sqrt(2*pi*s^2)) * exp(-((x-m)^2 / (2*s^2)))
nd(x)  = snd(x,0,1)
hnd(x) = snd(x,170,6)
wnd(x) = snd(x,63,10)

ベルカーブ（正規形）を計算する式を定義。mは平均値、sは標準偏差の積り。その一般式を使って、ndは標準形のカーブ。hndは、身長のカーブ。wndは体重のカーブを定義した。（haskellで言う所のカリー化みたい）

統計学では平均値の事を当たり前のように、μと表す。標準偏差はσだ。juliaは数学者にも平気で使って貰えるように、変数名とかにこれらのギリシャ文字を使える。

emacs-modeでも、 \mu とか入力してからTABを叩けば、μに変換出来る。ただ、残念な事にCUIモードで使うemacsでは、誤動作する。そう、ひと昔の某editorみたいに、byte単位でしか扱えないやつのごとくだ。カーソル位置の制御が正しくない。ユニコード文字が泣き別れになって、文字化けとか発生。ゆえに、涙を呑んで、ありきたりの文字を採用した。残念至極。

julia> x = 150:190
150:190
julia> @gp x [hnd(y) for y in x] "w l"

これが実際に身長のカーブを描く例。150から190cmの範囲に渡って、出現頻度をグラフにしてみた。 同様に体重のカーブを描いてみると、幅広になった。標準偏差が大きいと、幅が広がる（要するに、バラツキが大きい）のが伺える。

さて、次は、身長と体重をランダムに発生させて、そのデータからBMIを計算させればいいな。 そんな事出来るのか？　取り合えず using Random だな。

randn([rng=GLOBAL_RNG], [T=Float64], [dims...])

Generate a normally-distributed random number of type T with mean 0 and
standard deviation 1

randnってのが、偏った比率で乱数を発生してくれるとな。本当か検算してみる。

julia> q = randn(3000)
3000-element Array{Float64,1}:
 -0.5795832062760665
 -0.0897716967411955
   :
julia> using Statistics
julia> mean(q)
-0.023806058060213582
julia> std(q)
0.9946257193967464
julia> minimum(q)
-3.202317530254624
julia> maximum(q)
4.013185378506785

平均が０、標準偏差が１。範囲は-3から4ぐらい。

julia> h = hist(q)
Gnuplot.Histogram1D([-4.5, -3.5, -2.5, -1.5, -0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5], [0.0, 3.0, 64.0, 401.0, 1057.0, 1033.0, 376.0, 61.0, 4.0, 1.0, 0.0], 1.0)

julia> @gp h

一応目視確認。ベルカーブになってるね。でも幅の拡がりが±3ぐらいしかなくて、これじゃ、身長や体重の代用になりそうにないな。どうする？

q .+ 170

とかすれば、中心を並行移動出来るけど、拡がりはどうする。元データの標準偏差が１でこれが拡がりを表しているんだから、単純に掛け算しておけばいいのか。 先に拡げて、それを並行移動だな。

julia> r = q .* 6
julia> z = r .+ 170
julia> h = hist(z)
Gnuplot.Histogram1D([147.5, 152.5, 157.5, 162.5, 167.5, 172.5, 177.5, 182.5, 187.5, 192.5, 197.5], [0.0, 24.0, 124.0, 458.0, 919.0, 899.0, 439.0, 121.0, 13.0, 3.0, 0.0], 5.0)

うん、これでよさそうだ。一応検算って事で、meanとstdを取ってみたら、指定した通りになってた。

この方法を思いつたのは、正規形のやつの平均が０で標準偏差が１って点。両数字とも、特別な値だ。０は加算の単位元、１は乗算の単位元。

何も考えずに加算すれば、平均値がずれる。何も考えずに乗算すれば、標準偏差を簡単に変更できる。若し、元データ（乱数）が、それ以外になってたら、それを考慮して演算しなきゃいけない。そういう意味では、０と１ってのは、良く吟味された値だと思うよ。数学者の怠慢に乾杯（誉め言葉です、念の為）

後は、これを体重にも適用すれば、シュミレーション用のデータが完成。

bmi(h,w) = w / (h/100)^2

reqrn(m,s,n=1) = (randn(n) .* s) .+ m

reqrnは、平均値と標準偏差と個数（省略時は１）を指定して、乱数を発生させる関数だ。

 ht = reqrn(170,6,3000)
 wt = reqrn(63,10,3000)
 b = [bmi(ht[i],wt[i]) for i in 1:3000]
 h = hist(b)
Gnuplot.Histogram1D([9.0, 11.0, 13.0, 15.0, 17.0, 19.0, 21.0, 23.0, 25.0, 27.0, 29.0, 31.0, 33.0, 35.0], [0.0, 10.0, 39.0, 118.0, 294.0, 508.0, 600.0, 583.0, 416.0, 256.0, 118.0, 38.0, 20.0, 0.0], 2.0)

出来た。

julia> mean(b)
21.873364531202803

julia> std(b)
3.801185631721734

julia> maximum(b)
33.99864944391385

julia> minimum(b)
10.285352614719926

実データと比較してみると、まあ再現出来ているな。

上でBMIのシュミレーションをしたけど、極限値はどうなるだろう？

極大になるには、分子が体重が大きくて、分母の身長が小さければ良い。でぶでチビと言う取合せ。極小はその逆だから、瘦せでノッポさん。

julia> bmi(170-6*3, 63+10*3)
40.25277008310249

julia> bmi(170+6*3, 63-10*3)
9.336803983703033

平均から３σ離れた人は、1000人中に１人ぐらいしか現れない特異な人。身長も体重も何の相関が無いとして極端と極端を組み合わせると、WHO真っ青、厚生省の役人も飛んでくるような値になりました。

bmiの計算式を見てると、体重を定数と見れば、身長の２乗に反比例してる。体重を固定して、身長を変化させた時、どんなグラフになるか書かせてみるか。

Gnuplotで一つのグラフに複数の折れ線を入れる方法 Plot two curves: として、例が載ってた。これを応用すればいいんだな。

x = 150:190
@gp "set grid"
@gp :-  x [bmi(h,65) for h in x] "w l"
@gp :-  x [bmi(h,85) for h in x] "w l"
@gp :-  x [bmi(h,45) for h in x] "w l"

これで、身長を150から190まで振った時で、３種の体重の線を書けた。

次は３Dグラフだな。これだと、全体像が分かる。

julia> x = 150:5:190
150:5:190

julia> y = 45:5:85
45:5:85

julia> z = [bmi(h,w) for h in x, w in y]
9×9 Array{Float64,2}:
 20.0     22.2222  24.4444  26.6667  …  31.1111  33.3333  35.5556  37.7778
 18.7305  20.8117  22.8928  24.974      29.1363  31.2175  33.2986  35.3798
   :
 13.1483  14.6092  16.0701  17.531      20.4529  21.9138  23.3747  24.8356
 12.4654  13.8504  15.2355  16.6205     19.3906  20.7756  22.1607  23.5457

julia> @gsp x y z "w l"

x軸に5cm刻みの身長を割り当て。ｙ軸は5Kg刻みの体重を割り当て。ｚ軸方向はbmi値を計算させた。後は、@gspで、サーフェイスグラフを描くだけ。例によってマウスでぐりぐりとグラフを回せるぞ。

ちょっとおまけ。シュミレーション用のデータとして正規分布の乱数を 注意して 使った。現実に即してね。何も考えないと、シュミレーションなら（普通の）乱数でいいじゃんとなる。

rand(150:190, 3000)とかやると整数の乱数が返ってきちゃう。それは、あんまりだろう。小数点付きで欲しい（ちょっと知恵が回った状態）。さあ、どうする？

julia> q = (rand(3000) .* (190 - 150)) .+ 150
3000-element Array{Float64,1}:
 171.3243907166107
 155.6558037805561
  :

本当にまんべんなく散らばっているか？　例によって、検定

julia> mean(q)
170.06559037144146
julia> std(q)
11.486356303816363
julia> minimum(q)
150.00071148280264
julia> maximum(q)
189.99112901304898

よさそうだけど、まだ信用ならん。やっぱり度数分布を取ってみなければ安心出来ない。

julia> counts(q, 150:190)
ERROR: MethodError: no method matching counts(::Array{Float64,1}, ::UnitRange{Int64})
Closest candidates are:
  counts(::AbstractArray{T,N} where N where T<:Integer, ::UnitRange{T} where T<:Integer) at /home/sakae/.julia/packages/StatsBase/unDUx/src/counts.jl:79
 :

入力はIntしか受け付けないんだな。どうやってFloatをIntにする。ぐぐれ。

julia> qi = floor.(Int, q)
julia> @gp 150:190 counts(qi, 150:190) "w l"

多少のバラツキはあるけど、出現頻度はほぼ同じになった。明らかにベルカーブとは違う。

julia> rand(150.0:0.01:190.0, 5)
5-element Array{Float64,1}:
 181.29
 156.6
 165.44
 189.07
 175.62

ああ、自前でやらなくても、レンジを上手に使えばいいんだよ。これだから原始人は困るよ。