Kosugitti's BLOG

アンドロイドは正規分布の夢を見るか

この画面は、簡易表示です

R

習作;混合正規分布モデル(Gaussian Mixture Model) by rstan

Stanのマニュアルに載っている簡単なモデルなんだけど,最後の「各要素の所属確率」を算出するのになぜか手間取った。softmax関数を使うときの型の問題でした。ちぇ。

Stanコードはこの通り。generated quantitiesのところでsoftmax関数を使います。

data{
  int<lower =2> K; #number of clusters
  int<lower =1> N; #number of observations
  real X[N]; #observed data
}

parameters{
  vector[K] mu;
  vector<lower =0,upper=10>[K] sig2;
  simplex[K] theta;
}

transformed parameters{
  vector[K] ps[N];
  for(n in 1:N){
    for(k in 1:K){
      ps[n,k] = log(theta[k])+normal_lpdf(X[n]|mu[k],sig2[k]);
    }
  }
  
}

model{
  sig2 ~ cauchy(0,2.5);
  mu ~ normal(0,100);
  for(n in 1:N){
    target += log_sum_exp(ps[n]);
  }
}

generated quantities{
  simplex[K] u[N]; #class membership probability
  for(n in 1:N){
    u[n] = softmax(ps[n]);
  }
}

キックするコードは次の通り。

library(rstan)
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())

# 1 dimensional Mixed Gauss
#  K個の正規分布を混ぜる

K < - 3
mean1 <- 5
mean2 <- -5
mean3 <- 0
sd1 <- 1
sd2 <- 2
sd3 <- 3

N1 <- 500
N2 <- 300
N3 <- 200
N <- N1+N2+N3
X <-c(rnorm(N1,mean1,sd1),rnorm(N2,mean2,sd2),rnorm(N3,mean3,sd3))

### MCMC!
stanmodel <- stan_model("develop/gmm_single.stan",model_name="GMM")
standata <- list(K=K,N=N,X=X)
max.iter <- 3000
burn.in <- 1000
itr <- max.iter - burn.in
C <- 4


# vbは速いしラベルスイッチングとかないから嬉しい
fit_vb <- vb(stanmodel,data=standata)
print(fit_vb,pars=c("mu","sig2","theta"))



# 一つならラベルスイッチングは起きない
fit_sp <- sampling(stanmodel,data=standata,chain=1,iter=max.iter) 
# Rhatも優秀
print(fit_sp,pars=c("mu","sig2","theta"))
print(fit_sp,pars=c("u"))


# チェインを複数にするとラベルスイッチングが起きる 
fit_sp <- sampling(stanmodel,data=standata,chain=C,iter=max.iter,warmup=burn.in)
# Rhatが悪くなってもくじけちゃいけない
print(fit_sp,pars=c("mu","sig2","theta"))
# 目で見ると綺麗なラベルスイッチングが確認できるよ。
traceplot(fit_sp,pars="mu")


塾合宿のおもひで

お盆休みを満喫しております。

 

お休みに入る前,8/8-8/11の二泊三日で魁!ベイズ塾夏合宿2016に参加してきました。今年から新塾生も加入し,総勢19名というソコソコの集団サイズ。残念ながら,最近の若手は魁!男塾を知らないようで,そこかしこにある「塾のオマージュ」が伝わっていないようです。

二日間にわたって,StanとRを使いながら,ベイズ推定に色々な角度からアプローチ。日中は講義スタイル,夕食後はLTで自主的に研究発表をするというスタイルだったのだけど,さすが男気あふれる男塾,みんな自分で面白いことをやってますなあ。私はゼミ生のモデルの紹介と,状態空間モデルの話を披露しました。個人的にはアニメの公式ツイッターアカウントの,フォロアー数を自動的にとってこれるAPIの存在を知ったことがでかかったかな。いや,もちろん一番は国里先生の機械学習の話でしたけどね。

今後の活動について,どうしましょうかという話が最後にありましたけど,基本的には塾生のスキルアップ,塾からの情報発信だと思うので,塾長に何か面白いことを思いついてもらうかなぁ。最近はいろんなところでベイズの理解者が増えてきたとは思うので,目新しいことをするのが難しくなってきたかなとは思うね。



毒くらわば皿まで

rstan 2.10が来てますね。

2.9とちがって,式の代入が< -から=にすることが推奨されているようで,そうなっちゃうと今までのコードも全部書き換えないといけない。
とはいえこれから先はこうなるってんなら,早めに対応しておこうか,と母艦などのパッケージを思い切ってアップデート!

本艦のMacはすんなりアップデートでき,サンプルとして昔ながらの8schoolsを試したら,無事に走りました。もっともエラーとして,


DIAGNOSTIC(S) FROM PARSER:
Warning (non-fatal): assignment operator < - deprecated in the Stan language; use = instead.

って言われちゃいましたけどね。あー,やっぱりかー(当たり前だ)と。

で,書き直しなんだけど,よく見たら


theta[j] = mu + tau * eta[j];

だけでなく,モデルセクションも


model {
target += normal_lpdf(eta | 0, 1);
target += normal_lpdf(y | theta, sigma);
}

って変わってるのね。もとは


model {
eta ~ normal(0, 1);
y ~ normal(theta, sigma);
}

だったんですよ。

ともかくこれで走りました。コンパイルの時にちょっと警告が出たけどね。
これからは,コンパイル時のこの警告を潰していかないといけません。うう,やっと関数が出来上がったと思ったら・・・(苦笑)

[crayon-5bc6fba9b7ae2030158412/]

 

ついでに,とMacの上の仮想Windowsは,っとみたら,Rが3.2で止まってた。R3.3.1にあげて,rstanもアップデートする。

ここでRtoolsのバージョン違いがあってうまくいかなかった。最新版Rtools34はダメみたいで,Rtools33にしか対応してないみたい。でそこからが大変。コンパイル系のエラーが出まくって,rstanがうまく入りません!

仕方がないので,RもRtools33もRstudioもぜんぶアンインストール&再インストールしてから,install.packages("rstan")した。コンパイル時になんか警告が出たっぽいけど,これでやっとインストールできました。
ちなみにrstudioapiとかcolorなんとかってパッケージが入らないみたいで,要求されるままに手動で入力。ちょっと面倒だったな。今度学生に準備させるんだけど,事前に練習しておいてよかったよ・・・。

最後に Ubuntu。久しぶりに起動したので,システムそのもののアップデートをして,インストール。何だかうまくいかなかったんだけど,最終的にソースからコンパイルするとうまく走りました。Windowsより簡単ではあったな。

 

ともあれ,よいこのみんなは締め切りが迫っている中,環境を変えるのは良く無いってことを覚えておこうね!



Kandai.Rで発表してきました

6/18,第一回Kandai.Rで発表してきました。母校で行われる記念すべき第一回目に,初心者講習を引き受けさせていただきました。

Rとベイズは私の最後の恋です。Rを広めるためであれば,初心者講習をずっとやり続けることも厭わないのであります。上手か下手かは別にして・・・。

今回は「くれぐれも初心者向けに,やさしくね」と言われていたこともあって,非常に丁寧に進めたつもりです。一人でもファンが増えれれば嬉しいなあと思いながら。

Kandai R 入門者講習 from 考司 小杉

夜の情報交換会!では,妥当性としての「てっちゃんの手品」についてディスカッションができてよかったです。

久しぶりに大学以外の自由な雰囲気の中でのプレゼンテーションでした。楽しかったわ。



読書会を始めます

学内で,読書会を始めます。

[amazonjs asin=”4254122144″ locale=”JP” title=”はじめての 統計データ分析 ―ベイズ的〈ポストp値時代〉の統計学―”]
これをみんなで読んでいきます。

Q&Aを読み終わって,自分も何か始めないではいられなかった。半ば強引に,ゼミ生を誘い込みました。
私の考えでは,ベイズ統計学になると,もっと情報仮説が増えてくるはず。そうでなければならなかったんだ。そして,そういう仮説を思い付くには,帰無仮説を立てる訓練をしたように,統計的枠組みの中で仮説を発想する力が必要になる。

この本にあるように,発想を自由にできるようになれば,そのあとのアプローチは伝統的仮説検定の枠組みよりもずっと簡単なはず。
この読書会ではRプログラムを走らせたりすることをメインにはしない。ゆっくり,データとともに本を読んで,考える時間を持ちたいです。
一時間と短いですが,毎週ちょっとずつでも,できるだけ前に進んでいきたい。

読書会フライヤー



RstudioとパッケージとGithubと

大事な午前中がまるまる潰れてしまった。

Rstudioでパッケージを作るのが楽になったのはわかる。RstudioでGithubのバージョンコントロールができることも知ってる。ちょっとできる。Githubで文章を加筆修正したりもできる。

その組み合わせがうまくいかぬ。

1.RstudioでFile>New Project>バージョンコントロール>Git

とすると,GitとRstudioの連携はうまくいく。が,パッケージを作るプロジェクトの雛形は用意してくれない。ここでpackage.skeleton()が必要なのか?

2.RstudioでFile>New Project>パッケージを作る>Gitにもチェック

とすると,create git repositoryのチェックボックスをオンにするも,レポジトリが作られるわけでもなく,gitがうまくいかない。パッケージの雛形は用意してくれる。

3.本日最大の問題点

mkmeans関数というのを作り,mkmeansパッケージというのを作ろうとしたら,いちおうエラーなくbuildするところまではできるのだけど,build>library(mkmeans)をしてもmkmeans関数が見つからないって言われる。関数名とパッケージ名が同じなのがまずかったのか,と思って変えてみたけどうまくいかない。どこに問題が・・・?

 

自分の悩みを整理するために日記にしてみた。順番にやっつけていくか。まだ負けてないよ,奮闘記録ですよちくしょう。

 

ちなみに,成果はGistに置いた。

解決編

3.について

exportの指定をしておかないと,外から触れない=見えない,ということ。内部だけで使いたい関数は隠しておけということね。外に出力するわけではなくて,解放するって感じか。この指定をしておかないと,サンプルコードすら走らない。ちなみに@exportsとスペルミス?してもアウトでした。

1と2について

パッケージよりもGithubとの連携が欲しいので,1.のやり方(githubと連携してからスケルトン)のほうがいいと思う。2のやり方はgithubとの連携ができないもん。

余談

roxygen2ってすげえ。すげえ便利。すげえ楽。

できた。

https://github.com/kosugitti/imKmeans

まとめ

  1. Githubにレポジトリを作る
  2. Rstudioでバージョンコントロール・フォルダをつくる。Githubと紐付ける。
  3. package.skeltom()で雛形を作る
  4. 作って上げる。修正する。enjoy Git.


Lesson10.Bayesian Modeling 2

#2016.01.07現在,Rpubsが不調なため(?),こちらにアップします。

#2016.01.22 Rpubsが復旧したので,そちらにリンクを貼るようにします。なお,前に載せてあった記事にはstanコードにバグがありましたので,修正しました。

Lesson10. Bayesian Modeling 2 http://rpubs.com/kosugitti/145116



因子軸の回転の可視化

「因子軸を回転って結局どういうイメージですか」

という質問をされたのですね。
探索的因子分析の後は,解釈をしやすくするために軸を回転しますよ,と説明しているのだけど,そのイメージがわかりにくいということなので,Rでプロットしてみることにしました。

psychパッケージのデータセットbfiを例にあげます。これは五次元のデータなんですけど,仮に2因子だとしましょう。
(後で出てきますが3次元空間の回転行列から角度を求めるのは行列計算的にちょっと面倒なのです。)

因子数を2にし,回転なし,バリマックス回転,プロマックス回転の結果をオブジェクトに代入しておきます。

[crayon-5bc6fba9b9898611861288/]

[crayon-5bc6fba9b98a4435825614/]

[crayon-5bc6fba9b98a9343991777/]

直交回転の場合

さてここから加工。まず因子負荷行列をmatrix型で抜き出し,データフレーム型に変換。psychパッケージのfaオブジェクトはいろいろ付いているからね。無駄なものをそぎ落とします。

[crayon-5bc6fba9b98ae861872833/]

因子名をわかりやすくFx,Fyにしましょう(2次元だから)

[crayon-5bc6fba9b98b3860636906/]

回転なしの因子負荷量をつかって因子負荷行列のプロットを描いてみます。

[crayon-5bc6fba9b98b8309232026/]

 

unnamed-chunk-4-1

 

回転後の因子負荷行列の散布図も描いてみます。

[crayon-5bc6fba9b98bf514800543/]

unnamed-chunk-5-1

うむ,確かにこれでは分かりにくい。

さて,回転は行列をかけることで行っています。座標を変換するんですね。Wikiのページにもあるように,回転行列とはcos,sinをつかった行列です。つ回転行列

この回転行列はrot.matの中に入っています。

[crayon-5bc6fba9b98cb955437407/]

1行1列目は$\cos \theta$の値でしたね。この値から逆に角度を求めます。逆cosはacosです。

[crayon-5bc6fba9b98d4268002424/]

[crayon-5bc6fba9b98de136246751/]

0.71ラジアンぐらいですね

図にしましょう。まず回転前のプロット。

[crayon-5bc6fba9b98e9835446187/]

回転軸を引きます。角度を傾きにするため,tanをとっています

[crayon-5bc6fba9b98f2235470056/]

第2次元は第1次元と直角に交わるので,0.5piを足してtanをとっています。

[crayon-5bc6fba9b98fc984350747/]

unnamed-chunk-10-1

たしかに,元の青い軸のプロットよりも赤い軸のプロットの方が,点と軸の距離が近いようです。強調した,というのがちょっとは実感できるでしょうか。

斜交回転

promax回転の場合。最初の処理は同じです。

[crayon-5bc6fba9b9907162344607/]

回転行列を抜き出します。

[crayon-5bc6fba9b9910145265290/]

角度を抜き出します。

[crayon-5bc6fba9b991a663925068/]

さて,斜交回転の場合は軸が相関します。相関行列はPhiに入っています。相関係数は(\cos \theta) の結果なので,逆算(acos)して第一軸と第二軸の角度を求めます。

[crayon-5bc6fba9b9923219818076/]

これでプロット。

[crayon-5bc6fba9b992c918590725/]

unnamed-chunk-15-1

こんなもんでどうでしょう。

本当のこというと,適当な2軸を選んで描画し,軸をぐいっと動かして相関係数を見せるようなインタラクティブなR画面を作りたかったけど,それは技量と時間の限界ということで,誰か助けてエロい人。



初めてのスクレイピング

スクレイピングとぞいふもの

へー,webからデータを取ってくることをスクレイピングっていうんだ,というレベルから始めたのですがね。

毎年,データ演習の授業では色々なデータを使うのだけど,分布が偏ってるのもあって,適度に関係性が見られて,カテゴリカル変数もあって,それでいて身近な例で,というような丁度いいデータってのが欲しいもんなのですよ。

なので私はいつも,プロ野球データFreak さんのデータをコピペして,データを作っていました。
が,これが簡単にできるパッケージがあるらしい。のでやってみた。

[crayon-5bc6fba9ba8e2767913905/]

[crayon-5bc6fba9ba8f3831471965/]

同様に読み込んでデータ化していけばいいんだけど,残念ながらデータによっては数字じゃないところがあったり,最後にも変数名が入っていたりして,全部chr型になってたりする。

面倒だけど,これを修正する作業を行いました。

ピッチャーの年俸データを整える

文字列操作パッケージstringrの解説を参考に,文字列操作をしていく。

[crayon-5bc6fba9ba8ff969209149/]

[crayon-5bc6fba9ba909782000374/]

バッターの年俸データを整える

同じような作業。

[crayon-5bc6fba9ba912187580660/]

[crayon-5bc6fba9ba91c812424233/]

成績のデータを整える。

最後の一行に変なのが入っているせいで全部キャラクターになっちゃってる。

[crayon-5bc6fba9ba925022725428/]

成績と年棒を合わせる

ホームラン一本おいくら万円とか知りたいじゃないですか。
名前でマージします。

[crayon-5bc6fba9ba92f775018765/]

これで出来上がり。あとはwrite.tableでもなんでもどうぞー。



キングオブコントの評価軸

先週末に行われたキングオブコント2015。リアルタイムで観て,昨日また録画したものを見直して。
いやー,なかなかいい戦いをしていますね。

ところで今回は評価方法が変わって,5人の審査員が100点ずつ持ち寄って500点満点で評価するというスタイルになりました。
私も見ながら採点して行ったんだけど,振り返ると「世界観の完成度」「突発的な面白さ」「タイミングの良さ」「オリジナリティ」など,自分の中でも幾つかの評価次元があって,0-100の一次元に落とすのがとても難しいなと感じました。

で,審査員たちはどういう評価をしているんだろう?と思いますよね。
こういう話(「キングオブコント2015」採点データ分析で見た意外な真実。さまぁ〜ず三村が鍵を握っていた)がホットエントリーになってて,楽しく読ませてもらったんですが,残念,平均と標準偏差出したところまで。

では私の出番ですね(謎

まずデータを入れます。

[crayon-5bc6fba9bb505194290200/]

審査員も口々に(最初の演者を評して)「これが基準になっちゃうからね」とか「少し高めで始めちゃったから(大竹)」と言ってたように,個人の中での相対的な基準でしかないと思うので,これを相関行列に変えて,相関の高さ=近い,と距離行列にしたMDSをしてみます。

まずは相関行列。

[crayon-5bc6fba9bb515604981219/]

[crayon-5bc6fba9bb51d970353508/]

これを見ると,最低でも0.46ぐらいの相関はあるということだから,やはり皆さんお笑い芸人として共通するところを持っている,と判断していいんじゃないかな(目安として0相関が.3より大きくなると関係がある=独立ではない,とするのが社会科学では一般的かと)。
最大で三村さんと日村さんの0.84。コンビ同士の相関がバナナマンで0.64,さまぁ〜ずで0.54と中程度の相関があるのも妙に納得。

さて地図を描こう。

MDS-1

 

これを見ると,日村さんが意外と中心的というか,当たり障りのない評価をしていて,三村–設楽軸と松本–大竹軸があるんですなぁ。
コメントから考えると,三村さんは当日の笑いの多さ,設楽さんはここに至るまでの努力や練習量を評価するという違いがあったのかも。
松本さんと大竹さんは,キャラクター的には似ているのかと思ったけど,評価次元という意味ではオリジナリティ重視vs基本重視,という違いがあったのかしらん。

いや,正直ちょっと評価に迷います。
ほな評定者の分類だけじゃなくって,演者のスタイルと掛け合わせてみたらわかりやすくなるんでしょうか。やってみましょう。
行と列の両方を弄るのなら,双対尺度法,まあスコアが連続変数だから対応分析っていうべきかな。

[crayon-5bc6fba9bb527722727965/]

[crayon-5bc6fba9bb52f616075042/]

[crayon-5bc6fba9bb537861303987/]

 

unnamed-chunk-1-1

最初の二次元で80%の説明ができているから,これはこれでいいとして。
こうしてみると,やっぱり三村–設楽軸ははっきりしてて,設楽さんの「シュール・若手」好み,三村さんの「ベタ・渋・古参」好みが見られたかな。日村さんと大竹さんが意外と近くなるのね。

ボキャブラ天国の軸でいうと,横軸は シブい–インパクト ,縦軸が 知的–バカ で,松本さんのシブ知,日村さん・大竹さんのバカパク,設楽さんのインパク知,三村さんのバカシブ好み,ってところかなぁ。
で,第一次元の寄与率が62%あることから,やはりインパクトがあるかどうか(インパクトの反対はシブいと言うよりスタンダードかどうか,かな)が大きくて,

皆さんも自分の採点結果をここに入れて,楽しんでみてください。




top