読者です 読者をやめる 読者になる 読者になる

廿TT

譬如水怙牛過窓櫺 頭角四蹄都過了 因甚麼尾巴過不得

R {mixtools} で混合多項分布のパラメータ推定をやってみた

R

場面設定

  1. 全く同意できない
  2. 同意できない
  3. 同意できる
  4. 非常に同意できる

のような四択形式のアンケート20問を100人に対して行ったとする.

問題ごとに質問項目への回答数を集計してプロットしたのが以下の図である.

f:id:abrahamcow:20161102200359p:plain

アンケートの反応パターンから, 質問項目を分類したい.

こういうときに混合多項分布が使えるのではないか.

シミュレーション

mixtools の multmixEM 関数を使うと混合多項分布のパラメータ推定が行われる.

posterior probabilities が最大になるクラスを, そのアンケートの属するクラスとする.

シミュレーションで仮定した真のクラスと, 推定されたクラスを比べると, 正しく推定が行われていることがわかる.

アンケートの属するクラスによって色を変えてプロットしたのが次の図である.

f:id:abrahamcow:20161102201254p:plain

多項分布の確率パラメータも, シミュレーションで仮定した真値とおおむね近い値が求まっている.

f:id:abrahamcow:20161102201329p:plain

以下に R のコードを示す.

library(tidyr)
library(cowplot)
library(mixtools)
library(gridExtra)

prob1 <- c(0.1,0.2,0.3,0.4)
prob2 <- c(0.4,0.3,0.2,0.1)
set.seed(1)
clabel <- ifelse(rbinom(20,1,0.5)==1,1,2)
x <-t(sapply(clabel,
             function(x){if(x==1){rmultinom(1,100,p=prob1)}else{rmultinom(1,100,p=prob2)}}))
dat <-data.frame(x,clabel,id=1:20)
dat4plot <-gather(dat,item,value,-c(id,clabel))

p1 <-ggplot(dat4plot,aes(x=item,y=value,fill=factor(clabel)))+
  geom_bar(stat = "identity")+
  facet_wrap(~id)+
  theme(legend.position="none")+
  ggtitle("Correct answer")
fit1 <-multmixEM(x)
estlabel <-apply(fit1$posterior,1,which.max)
dat2 <-data.frame(x,estlabel,id=1:20)
dat4plot2 <-gather(dat2,item,value,-c(id,estlabel))

p2 <-ggplot(dat4plot2,aes(x=item,y=value,fill=factor(estlabel)))+
  geom_bar(stat = "identity")+
  facet_wrap(~id)+
  theme(legend.position="none")+
  ggtitle("Estimated")

grid.arrange(p1,p2,ncol=2)

comp1 <-data.frame(True =prob2, Estimated=fit1$theta[1,]) %>% 
  gather()
comp2 <-data.frame(True =prob1, Estimated=fit1$theta[2,]) %>% 
  gather()
p3 <-ggplot(comp1,aes(x=rep(1:4,2),y=value,fill=key))+
  geom_bar(stat = "identity",position = "dodge")+
  xlab("")+
  ggtitle("prob 1")
  
p4 <-ggplot(comp2,aes(x=rep(1:4,2),y=value,fill=key))+
  geom_bar(stat = "identity",position = "dodge")+
  xlab("")+
  ggtitle("prob 2")

grid.arrange(p3,p4)

参考文献:
https://cran.r-project.org/web/packages/mixtools/vignettes/mixtools.pdf

関連エントリ: