Chapter5 練習問題 Part2

(3)5.3.2項の続きで、アルバイトが好きかどうかごとにYを集計せよ

元データはこんなもの。

data-attendance-3.txt

PersonID,A,Score,Weather,Y
1,0,69,B,1
1,0,69,A,1
1,0,69,C,1
1,0,69,A,1
1,0,69,B,1
1,0,69,B,1
1,0,69,C,0
1,0,69,B,1
1,0,69,A,1

コードは以下のとおり。

#データ読み込み
d2 <- read.csv(file="data-attendance-3.txt")
#aggregate(x=list(<表示列名>=<変数名>$<集計列>),by=list(<表示列名>=<変数名>$<キー列>,<表示列名>=<変数名>$<キー列>,……),FUN=sum)

aggregate(x=list(Y=d2$Y),by=list(A=d2$A),FUN=table) #自分の回答

aggregate(Y~A,data=d2,FUN=table) #回答はこっち

結果はどちらも一致。

  A Y.0 Y.1
1 0 288 994
2 1 386 728

(4)5.3.3項では曇りと雨の影響を固定した。晴れを０とし、曇りと雨の影響をパラメーター化したモデルも考えられる。Stanで推定せよ。

まずモデル式を本から書き写すと

これに基づくStanプログラムは

model5-5.stan

data{
  int I;
  int<lower=0, upper=1> A[I];
  real<lower=0, upper=1> Score[I];
  real<lower=0, upper=1> W[I];
  int<lower=0, upper=1> Y[I];
}

parameters{
  real b[4];

}

transformed parameters{
  real q[I];
  for(i in 1:I)
    q[i]=inv_logit(b[1]+b[2]*A[i]+b[3]*Score[i]+b[4]*W[i]);
}

model{
  for(i in 1:I)
    Y[i]~bernoulli(q[i]);
}

これを修正して天気をパラメーターに入れる。

model5-5-2.stan

data{
  int I;
  int<lower=0, upper=1> A[I];
  real<lower=0, upper=1> Score[I];
  int<lower=1, upper=3> WID[I];
  int<lower=0, upper=1> Y[I];
}

parameters{
  real b[3];
  real bw2;
  real bw3;
}

transformed parameters{
  real q[I];
  real bw[3];
  bw[1]=0; //晴れは0
  bw[2]=bw2;
  bw[3]=bw3;
  for(i in 1:I)
    q[i]=inv_logit(b[1]+b[2]*A[i]+b[3]*Score[i]+bw[WID[i]]);
}

model{
  for(i in 1:I)
    Y[i]~bernoulli(q[i]);
}

実行コードは以下のとおり

conv <- c(1,2,3)
names(conv)<- c('A','B','C')
data <- list(I=nrow(d2),A=d2$A, Score=d2$Score/200, WID=conv[d2$Weather], Y=d2$Y)
fit4 <- stan(file="model5-5-2.stan", data=data, seed=1234)

トレースプロットは・・・

Stanコード間違って何度かやり直したけど、ようやくちゃんと収束しました。

続きはまた後日。