廿TT

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

ガンマ・ポアソン分布回帰による margarine データの分析(特に根拠のない推定法)

今日の川柳

ディリクレ・多項分布回帰による margarine データの分析 - 廿TT と同じデータを扱います。

説明変数の行列を X (N行D列)で表し、家庭 n におけるマーガリン k の購入数を y_{n,k} とします。

以下のモデルを考えます。

B=\exp(-X\beta)
z_{n,k} \sim \mathrm{Gamma}(\alpha, B_{n,k})
y_{n,k} \sim \mathrm{Poisson}(z_{n,k})

ここでの指数関数は要素ごとの指数関数です。B_{n,k} は B の (n,k) 成分です。

多項分布を使うとマーガリンの購入頻度そのものの情報は係数に反映されません。購入頻度を所与としたときのブランド選択の分析になります。

ポアソン分布を使うとマーガリンの購入頻度を直接分析できます。

y_{n,k} が与えられたときの z_{n,k} 事後分布はガンマ分布になります。

z_{n,k}|y_{n,k}\sim\mathrm{Gamma}(\alpha+y_{n,k}, B_{n,k}+1)

パラメータ推定の手順は以下のようになります。

1)
 \sum_n \sum_k \log \mathrm{Gamma}(z_{n,k}|\alpha, B_{n,k})
の上記の事後分布の下での期待値を求める。

2)
1で求めた関数を \alpha\betaに関して最大化する。

1、2を繰り返す。

\alpha\betaに関する変分推論の導出が難しかったの最尤推定で代用しました。

この推定法には特に理論的根拠がありません。

シミュレーション

この方法でパラメータがちゃんと推定できるか試します。

set.seed(1234)
K <- 10
beta <- matrix(rnorm(5*K),5,K)
N <- 100
x <- matrix(rnorm(N*5),N,5)
z <- matrix(rgamma(N*K,2,exp(-x%*%beta)),N,K)
y <- matrix(rpois(N*K,z),N,K)

VIP <-function(y,x,...,maxit=1000,tol=1e-8,seed=12){
  N <- nrow(y)
  K <- ncol(y)
  D <- ncol(x)
  set.seed(seed)
  theta0 <- runif(1+D*K)
for(i in 1:maxit){
  alpha <- exp(theta0[1])
  beta <- matrix(theta0[-1],D,K)
  B <- exp(-x%*%beta)
  Ez <- (alpha+y)/(B+1)
  Elogz <- digamma(alpha+y)-log(B+1)
  ffunc <-function(theta){
    alpha <- exp(theta[1])
    beta <- matrix(theta[-1],D,K)
    B <- exp(-x%*%beta)
    sum(alpha*log(B)+(alpha-1)*Elogz-B*Ez-lgamma(alpha))
  }
  opt<-optim(theta0,ffunc,...,control = list(fnscale=-1),method = "BFGS")
  if(all(abs(opt$par-theta0)<tol)){
    break
  }
  theta0 <- opt$par
}
  print(i)
  return(opt)
}
opt <- VIP(y,x)
exp(opt$par[1])
# 2.052756

plot(beta,opt$par[-1],xlim = c(-3,3),ylim = c(-3,3))
abline(0,1,lty=2)

推定された \beta と真の \beta のプロットです。

f:id:abrahamcow:20180816063149p:plain

そこそこ求まっている。

データ分析

ガンマ分布の期待値、\alpha/B (要素ごとの割り算)のプロットです。

f:id:abrahamcow:20180816063658p:plain

当てはまりはいまいちかなー。

係数 \beta のプロットです。

f:id:abrahamcow:20180816063914p:plain

前回の結果(ディリクレ・多項分布回帰による margarine データの分析 - 廿TT)と比べるとホワイトカラーが10を選ぶ強度が高いけど、どうなんだろう。

data(margarine,package = "bayesm")
library(tidyverse)
choice_hh <-margarine$choicePrice %>% 
  select(hhid,choice) %>% 
  left_join(margarine$demos) %>% 
  mutate(choice=factor(choice))

choice_s <-group_by(choice_hh, hhid, choice, Income, Fam_Size, college, whtcollar, retired) %>%
  tally() %>% 
  spread(choice,n,fill=0) %>% 
  ungroup() %>% 
  mutate(Income=Income/10)

W <-as.matrix(choice_s[,-c(1:6)])
X <-cbind(Intercept=1,as.matrix(choice_s[,2:6]))

fit <- VIP(W,X,maxit = 10000)

Wdf <-as.data.frame(W) %>% 
  mutate(row=row_number()) %>% 
  gather(col,value,-row) %>% 
  mutate(type="observed") %>% 
  mutate(col=factor(as.integer(col)))

beta <-matrix(fit$par[-1],ncol(X),ncol(W))


Bdf <-as.data.frame(exp(X%*%beta)) %>% 
  set_names(1:10) %>% 
  mutate(row=row_number()) %>% 
  gather(col,value,-row) %>% 
  mutate(type="fitted") %>% 
  mutate(col=factor(as.integer(col)))

jointdf <-bind_rows(Wdf,Bdf)

ggplot(jointdf,aes(x=col,y=row,fill=value))+
  geom_tile()+
  scale_fill_continuous(low="white",high="orange")+
  facet_wrap(~type)

betadf <-beta[-1,] %>% 
  as.data.frame() %>%
  set_names(1:10) %>% 
  mutate(name=colnames(X[,-1])) %>% 
  gather(category,value,-name) %>% 
  mutate(category=factor(as.integer(category)))

ggplot(betadf,aes(x=category,y=name,fill=value))+
  geom_tile(colour="black")+
  scale_fill_gradient2()+
  theme_bw(base_size=16)