20150103_ブログ用2

japanr.stan

data
{
    int<lower=0> N; //サンプルサイズ(ユニークな個人)
    int<lower=0> NT; //サンプルサイズ(全体)
    int<lower=0, upper=1> y[NT]; //目的変数（離反）
    int id[NT]; //ユニークな個人の数
    vector[2] x[NT]; //切片と価格
    vector[10] w[NT]; //環境変数
    matrix[N,4] z; //属性
}

parameters
{
    row_vector[2] alpha[N]; //次元はマーケティング変数の個数+1(切片)
    matrix[4,2] theta; //属性の数,（マーケティング変数の数+1）
    cov_matrix[2] V_alpha; //マーケティング変数の個数+1

    row_vector[10] beta[N]; //次元は環境変数の数
    matrix[4,10] lambda; //属性の数,環境変数の数
    cov_matrix[10] V_beta; //環境変数の数
}

model
{
//初期値の設定
    vector[8] mu1; //（属性の数）*（マーケティング変数の数+1）
    vector[8] s1_v;
    matrix[8,8] s1; //行数列数は（属性の数）*（マーケティング変数の数+1）
    vector[2] s_1_v;
    matrix[2,2] s_1; //次元はマーケティング変数の個数+1(切片)

    vector[40] mu2; //（属性の数）*（環境変数の数）
    vector[40] s2_v;
    matrix[40,40] s2; //行数列数は（属性の数）*（環境変数の数）
    vector[10] s_2_v;
    matrix[10,10] s_2; //（環境変数の数）

mu1[1] <- -1;
mu1[2] <- 0.5;
mu1[3] <- 0.5;
mu1[4] <- 0.5;
mu1[5] <- 0.5;
mu1[6] <- 0.5;
mu1[7] <- 0.5;
mu1[8] <- 0.5;

for (i in 1:8)
    s1_v[i] <- 100;

s1 <- diag_matrix(s1_v);

for (i in 1:2)
    s_1_v[i] <- 100;

s_1 <- diag_matrix(s_1_v);

mu2[1] <- -1;
mu2[2] <- 0.5;
mu2[3] <- 0.5;
mu2[4] <- 0.5;
mu2[5] <- 0.5;
mu2[6] <- 0.5;
mu2[7] <- 0.5;
mu2[8] <- 0.5;
mu2[9] <- 0.5;
mu2[10] <- 0.5;
mu2[11] <- 0.5;
mu2[12] <- 0.5;
mu2[13] <- 0.5;
mu2[14] <- 0.5;
mu2[15] <- 0.5;
mu2[16] <- -1;
mu2[17] <- 0.5;
mu2[18] <- 0.5;
mu2[19] <- 0.5;
mu2[20] <- 0.5;
mu2[21] <- 0.5;
mu2[22] <- 0.5;
mu2[23] <- 0.5;
mu2[24] <- 0.5;
mu2[25] <- 0.5;
mu2[26] <- 0.5;
mu2[27] <- 0.5;
mu2[28] <- 0.5;
mu2[29] <- 0.5;
mu2[30] <- 0.5;
mu2[31] <- -1;
mu2[32] <- 0.5;
mu2[33] <- 0.5;
mu2[34] <- 0.5;
mu2[35] <- 0.5;
mu2[36] <- 0.5;
mu2[37] <- 0.5;
mu2[38] <- 0.5;
mu2[39] <- 0.5;
mu2[40] <- 0.5;


for (i in 1:40)
    s2_v[i] <- 100;

s2 <- diag_matrix(s2_v);

for (i in 1:10)
    s_2_v[i] <- 100;

s_2 <- diag_matrix(s_2_v);

//尤度
    for (i in 1:NT)
        y[i] ~ bernoulli(inv_logit(dot_product(x[i],alpha[id[i]]) + dot_product(w[i],beta[id[i]])));

//パラメーターの事前分布（マーケティング変数）
    for (i in 1:N)
        alpha[i] ~ multi_normal(z[i]*theta,V_alpha);
    to_vector(theta) ~ multi_normal(mu1,s1);
    V_alpha ~ inv_wishart(10,s_1);

//パラメーターの事前分布（環境変数）
    for (i in 1: N)
        beta[i] ~ multi_normal(z[i]*lambda,V_beta);
    to_vector(lambda) ~ multi_normal(mu2,s2);
    V_beta ~ inv_wishart(20,s_2);
}