Kosugitti's BLOG

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

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

2016 / 9月

40歳

若い頃はまあ早く死んでもいいというぐらいの勢いで,無謀なことをするわね。
年寄りはなんか逆に命に執着するような印象があるわね。

どちらも正しいのであれば,折り返しである40ぐらいに,なんか考え方の転換があるはずだわね。

最近思うのは,自分が生きて来た!という自負が,実は周りの人に生かされてた!という感覚。
自分がやって来たこと,やれたと思っていることも,周りの人がサポートしてくれて,やらせてくれていたんだと。周りの人が「小杉がしたいようにできるように」と,状況を用意してくれていたんだと。

妻はね。20代からずっと一緒で,そうしてくれているのは知っているというか。わかっていて甘えているところがあって,ごめんなさいとかありがとうとしか言いようがないんだけども。
子供ができることで,こいつらも自由にさせてやろう,という気持ちが湧いてくるわね。

同じ頃に,学生の面倒を見るようになって。一緒に遊んでくれるのもありがたいんだけど,遊んでもらえるなあとか,遊ばせてやらねばなあとか,そういう観点が自分にでてくるようになって。感謝されたいんじゃなくて,こっちが勝手に上から目線で,遊ばせてやってるキャラになりたくなるんだよね。

もちろん,「俺がこんなに面倒見てやってるんだぞ」感は一人になった時に思うんだけど,トータルで見ると,先人にしてもらったことの恩返しでね。反省するやら情けないやらで自己嫌悪を経て,大人なふりをしようとするんだよね。

そうすると,死んでられないというか,ちゃんとしてやらななぁ,って思う。「てぇ出したんなら,さいごまでやんな」ってトトロのババアも言うてた。

こうやって,大人になっていくんだな。

でも,まだ,そっちにいきたくないって気持ちがあるんだよ。
みんなに迷惑をかけて,面倒を見てもらいたいって思うんだよ。暴れたいんだよ。
そこまでの才能・能力をちゃんと育てて来たか?今からまだ伸びると思ってるのか?って自問自答してます,毎日。

俺はここまでなのか。いや,まだ・・・



習作;ラベルスイッチング問題をなんとかする関数

前回の続き。Stanコードに対数尤度を出す部分を追加。。。
(どう見てもカッコ悪い。スマートに書き直したいものだ。ていうか,そもそもこれで合ってるのかな?)


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{
  real log_lik[N];
  simplex[K] u[N]; #class membership probability

  for(n in 1:N){
    u[n] = softmax(ps[n]);
    log_lik[n] = log_sum_exp(ps[n]);
  }

}

で,鎖が複数になるとラベルスイッチング問題が生じるんだけど,これを解決するlabel.switching関数を使ってみる。

## どこに分類されたか
allocK < - rstan::extract(fit_sp,pars="u")$u
# 各チェインごとに割り当てられたクラス
zChain <- matrix(ncol=N,nrow=itr*C)
for(i in 1:(itr*C)){
  zChain[i,] <- apply(allocK[i,,],1,which.max)
}

## MCMCで得られたパラメタ
mcmc.params <- rstan::extract(fit_sp,pars=c("mu","sig2","theta"))

# mcmc * num.of.clusters * num.of.params
mcmc.pars <- array(data=NA,dim=c(itr*C,K,3))
mcmc.pars[,,1]<- mcmc.params$mu
mcmc.pars[,,2]<- mcmc.params$sig2
mcmc.pars[,,3]<- mcmc.params$theta

# この関数はSJW法を使う時に必要になってくる。completeオプションに関数を渡さないといけないので!
complete.normal.loglikelihood<-function(x,z,pars){
  #	x: data (size = n)
  #	z: allocation vector (size = n)
  #	pars: K\times J vector of normal mixture parameters:
  #		pars[k,1] = mean of the k-normal component
  #		pars[k,2] = variance of the k-normal component
  #		pars[k,3] = weight of the k-normal component
  #			k = 1,...,K
  g <- dim(pars)[1] #K (number of mixture components)
  n <- length(x)	#this denotes the sample size
  logl<- rep(0, n)	
  logpi <- log(pars[,3])
  mean <- pars[,1]
  sigma <- sqrt(pars[,2])
  logl<-logpi[z] + dnorm(x,mean = mean[z],sd = sigma[z],log = TRUE)
  return(sum(logl))
}


# パッケージのお出まし
library(label.switching)
# 全方法試す
set <- c("STEPHENS", "PRA", "ECR", "ECR-ITERATIVE-1", "ECR-ITERATIVE-2", "AIC", "SJW","DATA-BASED")

# 対数尤度
log_liks <- rstan::extract(fit_sp,pars="log_lik")$log_lik
library(loo)
loo(log_liks)
waic(log_liks)

#ピボットは対数尤度が一番大きいやつにしてみた(なんでもいいと思う)
# MCMCサンプルごとに対数尤度が一番大きいやつを算出(apply部)
# 表の形で集計,大きい順に並べかえる(sort(table)部分)
# 表のラベルを取り出して(names部分),数字にする(as.numeric)
pivot <- as.numeric(names(sort(table(apply(log_liks,1,which.max)),decreasing=T)[1]))
# さあ実行
ls <- label.switching(method = set,
                      zpivot = zChain[pivot,],
                      prapivot = mcmc.pars[pivot, , ],                      
                      z = zChain,
                      K = 3, 
                      p = allocK,
                      mcmc=mcmc.pars,
                      complete=complete.normal.loglikelihood,
                      data = X)
# 各推定法で結果が一致したかどうか
ls$similarity
# 推定された所属クラスタ番号
ls$clusters

練習ここまで。



習作;混合正規分布モデル(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")


学会の続きに

末っ子が入院しまして。

先週末からお熱を出してて,月曜日に小児科に行ったところ,とりあえず喉が赤いので抗生物質を出されて。まあ三日ぐらいして熱が下がったので,はあやれやれ,と金曜日には幼稚園に行けたのですが(と言っても,金曜日は久しぶりの登園で緊張したのか,園が近づくにつれて足取りが重くなり,行かない〜行かない〜と泣いたりしたんだけど)。

金曜の夜から出張に出てたので,わからなかったんだけど,土曜日・日曜日にまた症状が悪化したらしく,月曜日(祭日)の朝には休日救急病棟へ。診察の結果「肺が真っ白だねえ」と言われたというほどの,ザ・肺炎で,入院することに。病室で再会した娘はもう,点滴はされてるわ,バイタルを測る機械のケーブルが出てるわ,酸素吸入マスクしてるわで,ぐったりベットに寝ている姿を見ている方がつらいほど。

とはいえ,上の子二人は翌日から普通に学校があるわけです。3人の子育てってのはワンオペではどうしても回らんシーンが出てくるけど,今回はなかなかヘヴィな方です(スケジュールの調整がね)。寝泊まりの付き添い付き,ってことは家族が分断して生活をしないといけないことだからねえ。

幸い,裁量労働制なので,今日は朝だけ仕事をしに行って,夜通し付き添ってくれてた妻と交代。1日入院すると熱も下がったし,ずいぶん元気回復した感じ。機械の数値もいいように思います。だからまあ,そんなに長引かない・・・かな,と期待してます。楽しみにしていた,園最後の運動会はおやすみ確定なのでかわいそうだけど。明日は仕事が埋まっているけど,どうするかなー。

台風はこちらに関係ない感じですが,とまれ,なかなかイベントの多い秋です。



JSSP2016参戦記

この週末(2016.09.17-18)にあった日本社会心理学会大会に参加してきました。会場は母校。前回この学会大会が母校であったときは実行委員だったなー,とか思ってたら,あれからもう10年経ってたのね。

学部の建物も建て替わって,H号館という名前になっている。社会学部棟とH号館はつながっていて,似て非なるもので。しかも前の建物のレイアウトのイメージに合わせよう,とすると「あれ,これはどっち向きなの。今俺はあっち向いているはず?」とかえって混乱しちゃいます。

この学会はいろいろなセッションが並行で走るから,発表の内容によってはあっち行ったりこっち行ったり,とバタバタしている間に午前中が終わり,正午からの打ち合わせにはちょっと遅れて行くことに。無事打ち合わせが終わったと思ったら,今度は自分のポスターを貼る番で,在籍責任時間の前にちょっと別の打ち合わせをして,発表会場に戻ったら連名発表者が説明してくれてて・・・と目まぐるしいお働き。さすがにポスター発表後,ちょっとぐったりしたので,久しぶりにお会いしたfjt先生とコーヒーブレイクしましたよ。

懇親会はお酒を飲むだけ,というわけにもいかない。というのもIISという国際宇宙ステーションだかサーバのソフトだかと勘違いしてしまいそうな,アイデアの他流試合。ご指名いただいたので,懇親会の会場で議論することに。会場の端っこにでもポスターがあるのかな,と思いきや,ポスター会場を分断するかのように中央にずらりとポスターが並ぶという・・・。

師匠が小走りに壇上に向かって,説明してくださったワインも飲みたかったので,その辺をウロウロしているとIISに行くのが少し遅くなりまして,まあ潤滑油も入っていたのでいろいろ思いついたこと(普段は言わないようなこと!)を喋ってしまいました。ごめんなさい。

懇親会の後は西宮北口の,昔よく行った店で仲間と二次会。明日もあるから,とか言いながら看板の時間まで飲んでました。

ところが翌日も朝から真面目にセッションに参加しちゃうんだな,これが。学会特有の足腰&肝臓のトレーニングで少し疲労感はありましたが,この時期頑張らないと学者としてどうなん,ということで午前中はたっぷり会場。お昼は仲間とランチして,そのあとは研究グループと打ち合わせする時間に。ほんと,普段とは比べ物にならないぐらい,研究という愉悦に浸りました。

実行委員の皆様,ありがとうございました。大変楽しく充実した二日間を送ることができました。書籍コーナーでも拙著をいろいろ宣伝させていただきました。ご購入いただいた皆様,ありがとうございました。

研究に関する印象として,自分の発表は浮いてしまうんじゃないかと思っていたけど,結構興味を持ってくれる人がいたり,幾つかの発表の中でベイズ的だなとか,データ生成モデルの観点が広がってきているな,というようにも思いました。懇親会でも「ちょっと興味があるんだけど手を出しにくくて・・・」みたいな話も聞きました。

Rとベイズは自分の最後の恋と決めているのですが,素敵な恋を皆さんとも分かち合いたいので,研究会などに呼んでいただければ,伝道師としてすぐに参上しますよ。Enjoy Bayes! May the Bayes be with you!

 

学会二日目の終わりにはゼミの同窓会(その前に時間があったので0次会)までたっぷりと,体と脳みそをフル回転させたのでした。自分ところの院生をほったらかしてたのはほんまごめんなさい。

でも,ほんと,楽しかったでーす。

[amazonjs asin=”4762829439″ locale=”JP” title=”研究論文を読み解くための多変量解析入門 応用篇: SEMから生存分析まで”]

[amazonjs asin=”4762829404″ locale=”JP” title=”研究論文を読み解くための多変量解析入門 基礎篇: 重回帰分析からメタ分析まで”]



BSJ2016参加報告

行動計量学会第44回大会@札幌学院大学に,フル参加してきました!

例年,初日午後のチュートリアルセミナーはスケジュールの都合で参加できなかったりするんだけど,今年はベイズの話だしいけるぞ!とおもっていたら台風10号さんが大接近。羽田まで飛ぶかな?と心配になって便を変えてもらえないか,生協や事務に確認したところ,「じゃあ今日の17時の便で前入りしなさい」とのこと。これが決まったのが13時頃。慌てて帰宅して着替えを詰め,大学で資料を印刷して前入り。幸いホテルは駅近で,予定していたところの前日分追加で取れたのでずっと滞在できるし,ホテルの地下が居酒屋でルームーキーを持って行ったら20%オフだったし,チュートリアルセミナーは翌日午後からだから,体は疲れたものの結果オーライかな。

チュートリアルセミナーでしっかり勉強し,夜は馴染みのメンツで寿司を食べる。二日目は非対称セッションで非対称飲み会に参加,三日目は懇親会+懐かしの顔も加えて二次会。とにかく胃と肝臓と,もちろん脳みそフル回転の学会活動でした。

こうやって書くと「夜は飲んでばっかりじゃないか」と思われるかもしれませんが,そこで研究のアイデアを話し合ったり,普段会えない人からアイデアをもらったりして,とてもいい刺激になるのですよ。さすがに最終日は体力が尽きかけたけど,気合いを振り絞って飲みに行ったもんね。

学会はいつも余裕のあるスケジューリングで,プログラムが並列すぎないように,詰め込みすぎないようになっているから,空いているところで月末締め切りの仕事を片付けたりすることもできて,ほんとにこの学会はいい学会です。

二日目の自分の発表では,院生の頃のような「完璧にわかりきってないのに勢いで発表しちゃった」感が丸出しの話で,いろいろな人に突っ込まれても「この本にこう書いてあるからやってみたんすよ」ぐらいのお返事しかできなかったし,「こういう研究がありますよね」とかいわれても「へえそうなんですか」みたいな返事になって,いやあ恥ずかしいやら疲れるやらでした。でも,

のびしろですね!

頑張りますよ。刺激をもらった分,頑張り甲斐があるというか,発奮してるね。これが本当に学会に参加するいい効果なんだと思います。

やるぞ。やるぞ。いろいろやるぞ(仕事ももらってしまったし)。あー楽しかった。

 

追伸 次来るときはジンギスカンとビールもやりたいな。いや,海鮮と日本酒もよかったんだけどね。じゃがいも焼酎もね。

追伸2 現在帰りの飛行機が遅延中。行きも帰りもいろいろある旅だな,今回は。




top