投稿

11月, 2017の投稿を表示しています

R で正規分布の確率密度函数を描いてみよう

イメージ
できるだけ無駄の少ない PDF 図が欲しいときには、R が良いですね。 R で正規分布の確率密度函数を描いてみましょう。 x できあがりは下図です。 qnorm(0.025) で、正規分布時の 97.5 パーセンタイルが求まります。 縦の線2本は、2.5 パーセンタイルと、97.5 パーセンタイルです。 y=0 を追加して描いてあります。

t 検定は n がどのくらい大きくなれば、正規分布に近づくのか?

イメージ
前回の julia で t 分布と正規分布の確率密度函数 PDF を描いてみよう を見ていると、n=20 くらいで、かなり t 分布は正規分布に近いような感じがします。 実用上、気になるのは、両側検定の時に検定率 5% となる 97.5 パーセンタイルの t 値が、正規分布の時の 1.9599639845400583 にどのくらい近いか? ということではないでしょうか。 そこで、julia + pyplot で、n=3 から 100 までの自由度を x 軸に、それぞれの 97.5 パーセンタイルを y 軸に描いて、正規分布の時の 1.9599639845400583 を点線で描いてみましょう。 using Distributions using PyCall using PyPlot @pyimport seaborn as sns x = collect(3:100) y=Array[] y= quantile(TDist(3), 0.975) for i in 4:100 y = vcat(y, quantile(TDist(i), 0.975)) end plot(x,y) n_y = quantile(Normal(),0.975) hlines(n_y, 0, 100 , linestyles="dashed") ylim(0, 3.5) できあがりは、下図です。 n=100 のとき、97.5 パーセンタイルは、1.98397 です。0.02 程がどれほどの差かは、時と場合によるのでしょうか。n=20 くらいでは、かなり遠く感じますね。 quantile(TDist(3), 0.975) で、自由度 3 の時の 97.5 パーセンタイルが求まります。

julia で t 分布と正規分布の確率密度函数 PDF を描いてみよう

イメージ
julia で t 分布と正規分布の確率密度函数 probability density function、PDF を描いてみましょう。 確率密度函数は、Distributions パッケージの中にあるので、Distributitons を Pkg.add("Distributitons") で、インストールしましょう。 正規分布も t 分布も Univariable Distributions の中の Continuous Distributions というところで定義されています。 正規分布は、Normal() 、自由度 n の t 分布は、TDist(n) で良いようです。 using Distributions using PyCall using PyPlot @pyimport seaborn as sns x=collect(-4:0.1:4) fig, ax = subplots() for i in [3,5,10,20] y = pdf.(TDist(i), x) plot(x, y, alpha=0.6, label="n = $i") end y = pdf.(Normal(), x) plot(x, y, alpha=0.6, label="Normal Distribution") ax[:legend]() julia 0.6 と、PyCall.pyversion v"2.7.13" だと、これで、下記のような図ができます。 n が増えると、正規分布に近付くの図です。 x は、[-4:0.1:4] では動かずに、collect()の中に入れると動きます。 x=linspace(-4:4,801) でも動きます。 fig, ax = subplots() という1行を入れると、label を貯めて、ax[:legend]() で1ヶ所に描画してくれます。 ax[:plot](...) という描き方では、julia 0.6 では上手く動きませんでした。plot(...) と単純な方が良いようです。 ちなみに、seaborn は趣味なので、不要ですね。 たくさん重ねて描きたい時には、pypl...

Vim で、英単語数を調べる方法

:%s/\i\+/&/gn これで、重いアプリケーションを立ち上げずに済みます。 vimで単語数や行数を数える方法 より。 なお、ハイライトを消すには、 :noh

Rで、サンプルデータ作製からヒートマップ表示とクラスタリング表示 (Afiinity propagation含む) まで

イメージ
Rで、サンプルデータ作製して、ヒートマップ表示と、k-medoids (pam) と Afiinity propagation含むいくつかのクラスタリング表示を作る機会がありました。せっかくなので、記載します。今回ようなデータの処理と目的では、k-means のような方法が適していました。DBSCAN などは他をご参考をいただけますようお願いします。ひたすら駆け足で、忘備録程度の内容です。 対象は、50人分、10種類の観測値のサンプルデータです。 目的は、ヒトごと、10種類の観測値ごとのそれぞれのクラスターを見付けることです。 平均値とSD の異なる正規分布で、50人分、10種類の観測値のサンプルデータを作ります。 A1=cbind( matrix(rnorm(50, -1, 1),nc=5), matrix(rnorm(50,0,1),nc=5) ) A2=cbind( matrix(rnorm(30,-0.5,1),nc=3), matrix(rnorm(50,0,2),nc=5), matrix(rnorm(20,-3,1),nc=2) ) A3=cbind( matrix(rnorm(20,0.5,2),nc=2), matrix(rnorm(30,-2,1),nc=3), matrix(rnorm(20,0.2,0.2),nc=2), matrix(rnorm(30,1,1),nc=3) ) A4=cbind( matrix(rnorm(10,2,0.8),nc=1), matrix(rnorm(40,1,0.7),nc=4), matrix(rnorm(10,-0.5,0.1),nc=1), matrix(rnorm(40,0,1),nc=4) ) A5=cbind( matrix(rnorm(10,-0.1,0.2),nc=1), matrix(rnorm(40,-2,1),nc=4), matrix(rnorm(10,2,1),nc=1), matrix(rnorm(40,0.4,1),nc=4) ) M_Example = rbind(A1,A2,A3,A4,A5) データを50人の個人をランダムに並べ変えます。 クラスタライングの為には...