パフォーマンス一覧にSASが無かったのでやってみた
Twitterを見ていたら流れてきた話題
Juliaのコードを勝手に最適化してみた - りんごがでている
スライド内のパフォーマンス一覧にSASが入っていないので、やってみました。*1
書いたコード
実行ログ
pi=3.14144144 NOTE: DATAステートメント処理(合計処理時間): 処理時間 5.17 秒 CPU時間 5.17 秒
何度か実行してみるとだいたい5.2秒くらいになるみたい。*2
まとめ
SAS雲丹(SAS University Edition)、仮想OS上で動作する造りなので正直パフォーマンスについてはあまり期待していなかったんですが、悪くないですね。
みんな使おうSAS雲丹!
Processingで正規分布の(手動)最尤推定ツール
はじめに
Processing Advent Calendar 2014 17日目の記事です。
最近は便利な統計解析ツールがたくさん出てきているので、最尤推定によってパラメータを求めるのはとても手軽にできるようになりましたが、確率分布の形をグリグリ動かしながら対数尤度が変化する様を見ることができたら面白いんじゃないかと思い、こんなツールを作ってみました。
(手動)最尤推定ツール
以下のデータが観測されたとします。
1.364 0.235 -0.846 -0.285 -1.646
l
で平均(mu)が+0.01
j
で平均(mu)が-0.01
k
で標準偏差(sigma)が+0.01
i
で標準偏差(sigma)が-0.01
と正規分布のパラメータが変化するので、対数尤度(Log likelihood)が一番大きくなるポイントを探ってみましょう。
(グラフ上をクリックしてからキーボード操作してみて下さい)
コード
コードは以下のようになります。
ArrayList myPoints; float[] myData = {1.364, 0.235, -0.846, -0.285, -1.646}; float mu = 0; float sigma = 1; float margin = 20; void setup(){ size(540, 350); smooth(); noStroke(); myPoints = new ArrayList(); for (int i=0; i<myData.length;i++){ myPoint mp = new myPoint(myData[i]); myPoints.add(mp); } } void draw(){ background(255); stroke(0); noFill(); rect(margin, margin, width - margin * 2, height - margin * 2); textSize(12); for(int i=-5; i<=5; i++){ float xl = map(i, -5, 5, margin, width - margin); float yl = height - 5; text(str(i), xl, yl); } text("0", 5, height - margin); text("0.5", 1, height / 2); text("1", 5, margin); for(int i=0; i<width*10; i++){ float mi1 = map(i, 0, width*10, -5, 5); float l = (1 / (sqrt(2 * PI) * sigma)) * exp(-1/(2 * pow(sigma, 2)) * pow(mi1 - mu, 2)); float ml = map(l, 0, 1, height - margin , 0 + margin); float mi2 = map(mi1, -5, 5, 0 + margin, width - margin); noStroke(); fill(100); ellipse(mi2, ml, 1, 1); } float lh = 0; for(int i=0; i<myPoints.size();i++){ myPoint mp = (myPoint) myPoints.get(i); mp.update(); lh += mp.lh; } fill(100); textSize(15); text("mu = " + nf(mu, 1, 2), 420, 50); text("sigma = " + nf(sigma, 1, 2), 400, 70); text("Log likelihood = " + nf(lh, 1, 2), 350, 90); } class myPoint{ float xPos, yPos, lh; myPoint(float x){ xPos = x; } void update(){ yPos = (1 / (sqrt(2 * PI) * sigma)) * exp(-1/(2 * pow(sigma, 2)) * pow(xPos - mu, 2)); float mx = map(xPos, -5, 5, 0 + margin, width - margin); float my = map(yPos, 0, 1, height - margin , 0 + margin); stroke(80, 180, 90); line(mx, my, mx, height - margin); noStroke(); fill(80, 180, 90); ellipse(mx, my, 4, 4); ellipse(mx, height - margin, 6, 6); lh = log(yPos); } } void keyPressed(){ if(key=='i'){ sigma -= 0.01; if(sigma<=0){ sigma = 0; } }else if(key=='k'){ sigma += 0.01; }else if(key=='j'){ mu -= 0.01; }else if(key=='l'){ mu += 0.01; } }
人の手によって求められたパラメータ、機械的に作られたものには無い温かみを感じますね。
— ぷるうぬす (@Prunus1350) November 6, 2014
SASで「データ解析のための統計モデリング入門」を読む 第2章
まえおき
大盛況のうちに幕を閉じた「データ解析のための統計モデリング入門」 読書会 - connpass(通称みどりぼん)。私は最終回にUst参加しただけなのですが、とても楽しそうだったので私も独りでもくもく読んでみることにしました。(通称ひとりぼん)
ただ読んでもなかなか頭に入らないので、SAS University Editionを使って図や数値の再現を行いながら読んでいます。
この記事はSASのプログラムとその実行結果を羅列していくものです。(コードはこちらPrunus1350/midoribon · GitHub)
実行結果の解釈についての解説は特に行わないので、内容についてはみどりぼんをご覧下さい。
データ解析のための統計モデリング入門――一般化線形モデル・階層ベイズモデル・MCMC (確率と情報の科学)
- 作者: 久保拓弥
- 出版社/メーカー: 岩波書店
- 発売日: 2012/05/19
- メディア: 単行本
- 購入: 16人 クリック: 163回
- この商品を含むブログ (21件) を見る
第2章 確率分布と統計モデルの最尤推定
2.1 例題:種子数の統計モデリング
データ読込
%let home_dir = /folders/myfolders/midoribon/chapter02/; data data; infile "&home_dir.data.csv" dsd missover firstobs = 2; attrib y length = 8 label = "種子数"; input y; run;
P15 dataの内容を表示
proc print data = data; run;
実行結果
OBS | y |
---|---|
1 | 2 |
2 | 2 |
3 | 4 |
4 | 6 |
5 | 4 |
以下略
P15 データ数の確認
data _null_; set data nobs = n; put n; stop; run;
実行ログ
50
P16 最小値・四分位点・標本平均・中央値・最大値
proc means data = data min q1 median mean q3 max; run;
実行結果
分析変数 : y 種子数 | |||||
---|---|---|---|---|---|
最小値 | 下側四分位点 | 中央値 | 平均 | 上側四分位点 | 最大値 |
0 | 2.0000000 | 3.0000000 | 3.5600000 | 5.0000000 | 7.0000000 |
SASとRとではQuantileの計算方法が違うため、みどりぼんの値が再現できない http://en.wikipedia.org/wiki/Quantile
【おまけ】5通りのQuantile計算方法を試す
%macro quantile_test; %do i = 1 %to 5; proc univariate data = data pctldef = &i. noprint; output out = temp min = min q1 = q1 median = median mean = mean q3 = q3 max = max; run; proc print data = temp; var min q1 median mean q3 max; run; %end; %mend quantile_test; %quantile_test;
実行結果
OBS | min | q1 | median | mean | q3 | max |
---|---|---|---|---|---|---|
1 | 0 | 2 | 3 | 3.56 | 4.5 | 7 |
OBS | min | q1 | median | mean | q3 | max |
---|---|---|---|---|---|---|
1 | 0 | 2 | 3 | 3.56 | 5 | 7 |
OBS | min | q1 | median | mean | q3 | max |
---|---|---|---|---|---|---|
1 | 0 | 2 | 3 | 3.56 | 5 | 7 |
OBS | min | q1 | median | mean | q3 | max |
---|---|---|---|---|---|---|
1 | 0 | 2 | 3 | 3.56 | 5 | 7 |
OBS | min | q1 | median | mean | q3 | max |
---|---|---|---|---|---|---|
1 | 0 | 2 | 3 | 3.56 | 5 | 7 |
みどりぼんと一致する結果は得られませんでした。
P16 度数集計
proc freq data = data; tables y / missing nocol norow nopercent; run;
実行結果
種子数 | ||
---|---|---|
y | 度数 | 累積 度数 |
0 | 1 | 1 |
1 | 3 | 4 |
2 | 11 | 15 |
3 | 12 | 27 |
4 | 10 | 37 |
5 | 5 | 42 |
6 | 4 | 46 |
7 | 4 | 50 |
P17 図2.2 例題の種子数データのヒストグラム
proc sgplot data = data; histogram y / scale = count binwidth = 1; run;
実行結果
P17 標本分散・標本標準偏差
proc means data = data noprint; output out = sample_variance(drop=_type_ _freq_) var = var std = std; run; data _null_; set sample_variance; put var = ; put std = ; sqrt_var = sqrt(var); put sqrt_var = ; run;
実行ログ
var=2.986122449 std=1.72804006 sqrt_var=1.72804006
2.2 データと確率分布の対応関係をながめる
P19 種子数がyであると観察される確率を生成
data pois; do y = 0 to 9; prob = pdf("poisson", y, 3.56); output; end; run;
P20 図2.3 平均3.56のポアソン分布の確率分布
proc print data = pois; run;
実行結果
OBS | y | prob |
---|---|---|
1 | 0 | 0.02844 |
2 | 1 | 0.10124 |
3 | 2 | 0.18021 |
4 | 3 | 0.21385 |
5 | 4 | 0.19033 |
6 | 5 | 0.13551 |
7 | 6 | 0.08040 |
8 | 7 | 0.04089 |
9 | 8 | 0.01820 |
10 | 9 | 0.00720 |
P20 図2.4 平均λ=3.56のポアソン分布
proc sgplot data = pois; series x = y y = prob / markers lineattrs = (pattern=dash); run;
実行結果
P21 図2.5 観測データと確率分布の対応
proc freq data = data; tables y / missing noprint out = data2(drop=percent); run; data data3; set data2; prob = pdf("poisson", y, 3.56) * 50; run; title "Histogram of data"; proc sgplot data = data3; vbar y / response = count; vline y / response = prob markers lineattrs = (pattern=dash); run; title;
実行結果
2.3 ポアソン分布とは何か?
P23 図2.6 さまざまな平均(λ)のポアソン分布
data pois2; do y = 0 to 20; prob1 = pdf("poisson", y, 3.5); prob2 = pdf("poisson", y, 7.7); prob3 = pdf("poisson", y, 15.1); output; end; run; proc sgplot data = pois2; series x = y y = prob1 / markers lineattrs = (pattern=dash); series x = y y = prob2 / markers markerattrs = (symbol=diamond) lineattrs = (pattern=dash); series x = y y = prob3 / markers markerattrs = (symbol=triangle) lineattrs = (pattern=dash); yaxis label = "prob"; run;
実行結果
2.4 ポアソン分布のパラメーターの最尤推定
%macro m2_7; * ポアソン分布のパラメータ ; %let lambda1 = 2.0; %let lambda2 = 2.4; %let lambda3 = 2.8; %let lambda4 = 3.2; %let lambda5 = 3.6; %let lambda6 = 4.0; %let lambda7 = 4.4; %let lambda8 = 4.8; %let lambda9 = 5.2; data data4; set data2; %do i = 1 %to 9; prob&i. = pdf("poisson", y, &&lambda&i..) * 50; %end; run; * パラメータ毎に対数尤度を計算 ; data data5; set data4; %do i = 1 %to 9; log_l&i. = (y * log(&&lambda&i..) - &&lambda&i.. - log(fact(y))) * count; %end; run; * 対数尤度をマクロ変数に格納 ; proc sql noprint; select %do i = 1 %to 8; sum(log_l&i.) format = 8.1, %end; sum(log_l9) format = 8.1 into %do i = 1 %to 8; :log_l&i., %end; :log_l9 from data5; quit; proc template; define statgraph poisson_graph; begingraph; layout lattice / rows = 3 columns = 3; %do i = 1 %to 9; layout overlay / xaxisopts = (display=(tickvalues)) yaxisopts = (linearopts=(viewmin=0 viewmax=15 tickvaluelist=(0 5 10 15)) display=(tickvalues)); layout gridded / border = false halign = right valign = top; entry halign = left "lambda=&&lambda&i.."; entry halign = left "log L=&&log_l&i.."; endlayout; barchart x = y y = count; seriesplot x = y y = prob&i. / display = all markerattrs = (symbol=circle) lineattrs = (pattern=dash); endlayout; %end; endlayout; endgraph; end; run; * P26 図2.7 ; proc sgrender data = data4 template = poisson_graph; run; %mend m2_7; %m2_7;
実行結果
P28 図2.8 対数尤度
%macro m2_8; proc datasets library = work nolist; delete log_likelihood; run; quit; %do i = 20 %to 50; %let lambda = %sysevalf(&i. / 10); data __m1; set data; log_l = y * log(&lambda.) - &lambda. - log(fact(y)); run; proc summary data = __m1; var log_l; output out = __m2(drop=_type_ _freq_) sum = ; run; data __m2; set __m2; lambda = &lambda.; run; proc append base = log_likelihood data = __m2; run; %end; proc genmod data = data; model y = / dist = poisson link = id; ods output parameterestimates = __m10; run; data _null_; set __m10; if parameter eq "Intercept" then call symputx("est_param", estimate); run; proc sql noprint; select max(log_l), min(log_l) into :maxy, :miny from log_likelihood; quit; data sganno; retain function "line" drawspace "datavalue" linepattern "dash"; x1 = &est_param.; x2 = &est_param.; y1 = &miny.; y2 = &maxy.; output; run; proc sgplot data = log_likelihood sganno = sganno; series x = lambda y = log_l; run; %mend m2_8; %m2_8;
実行結果
2.4.1 擬似乱数と最尤推定値のばらつき
%macro mle_poisson(lambda=, rep=); proc datasets library = work nolist; delete __parameter; run; quit; ods listing close; ods html close; %do i = 1 %to &rep.; data __m1; do i = 1 to 50; y = rand("poisson", &lambda.); output; end; run; proc genmod data = __m1; model y = / dist = poisson link = id; ods output parameterestimates = __m2; run; data __m3(keep=estimate); set __m2; where parameter eq "Intercept"; run; proc append base = __parameter data = __m3; run; %end; ods html; ods listing; %mend mle_poisson; %mle_poisson(lambda=3.56, rep=3000);
P30 図2.9
proc sgplot data = __parameter; histogram estimate / binwidth = 0.1 scale = count; run;
実行結果
U-NOTEまとめ部さんにまとめ記事を書いて頂きました
第3回 「はじめてのパターン認識」 読書会のまとめ記事をU-NOTEまとめ部さんに書いて頂きました。
http://u-note.me/note/47485572
ベイズの識別規則は、これを見て学べ!- 第3回「はじめてのパターン認識」読書会まとめ【CodeIQ提供】 #はじパタ http://t.co/cYRLhwYEYr
— U-NOTE編集部 (@u_note) July 31, 2013
本編発表担当の@_kobackyさん @millionsmileさん、
LT発表の@wwackyさん @yamakatuさん @who_you_meさん @kenchan0130_akiさん、
みなさんのお陰でステキな記事になりました。ありがとうございました。
8/1(木)の時点でU-NOTEデイリーランキング1位になってます。
いえ~い!
第4回 「はじめてのパターン認識」 読書会開催しました
第3回 「はじめてのパターン認識」 読書会開催しました
7/16(火)に株式会社リクルートキャリアにて、第3回 「はじめてのパターン認識」 読書会を開催しました。
次回は7/30(火)にニフティ株式会社セミナールームにて開催予定です。
読書会の後に懇親会を予定しておりますので、是非ご参加下さい。
第2回 「はじめてのパターン認識」 読書会開催しました
7/2(火)に株式会社リクルートキャリアにて、第2回 「はじめてのパターン認識」 読書会を開催しました。