Stanで多次元正規分布

• 多次元データが取れる場合と、多次元だけれど、その一部次元のデータしか取れない場合っていうのは、どういう違いがあるのか…
• multinormal2.stan

data {
  int N_subjects;
  int N_items;
  matrix[N_subjects,N_items] y;
}

parameters {
  vector[N_items] mu;
  cov_matrix[N_items] Sigma;

}

transformed parameters {
}

model {
  matrix[N_items,N_items] identity;
  mu ~ normal(0,100);
  identity <- diag_matrix(rep_vector(1.0,N_items));
  Sigma ~ inv_wishart(N_items,identity);
  for(i in 1:N_subjects)
    y[i] ~ multi_normal(mu,Sigma);
}

• multinormal4.stan

data {
  int N_subjects[3];
  int N_items;
  matrix[N_subjects[1],2] y12;
  matrix[N_subjects[2],2] y23;
  matrix[N_subjects[3],2] y31;
}

parameters {
  vector[N_items] mu;
  cov_matrix[N_items] Sigma;

}

transformed parameters {
  cov_matrix[2] Sigma12;
  cov_matrix[2] Sigma23;
  cov_matrix[2] Sigma31;
  vector[2] mu12;
  vector[2] mu23;
  vector[2] mu31;
  
  Sigma12[1,1] <- Sigma[1,1];
  Sigma12[1,2] <- Sigma[1,2];
  Sigma12[2,1] <- Sigma[2,1];
  Sigma12[2,2] <- Sigma[2,2];
  
  Sigma23[1,1] <- Sigma[2,2];
  Sigma23[1,2] <- Sigma[2,3];
  Sigma23[2,1] <- Sigma[3,2];
  Sigma23[2,2] <- Sigma[3,3];
  
  Sigma31[1,1] <- Sigma[3,3];
  Sigma31[1,2] <- Sigma[3,1];
  Sigma31[2,1] <- Sigma[1,3];
  Sigma31[2,2] <- Sigma[1,1];
  
  mu12[1] <- mu[1];
  mu12[2] <- mu[2];
  mu23[1] <- mu[2];
  mu23[2] <- mu[3];
  mu31[1] <- mu[3];
  mu31[2] <- mu[1];
}

model {
  matrix[N_items,N_items] identity;

  mu ~ normal(0,100);
  identity <- diag_matrix(rep_vector(1.0,N_items));
  Sigma ~ inv_wishart(N_items,identity);
  
  for(i in 1:N_subjects[1]){
    y12[i] ~ multi_normal(mu12,Sigma12);
  }
  for(i in 1:N_subjects[2]){
    y23[i] ~ multi_normal(mu23,Sigma23);
  }
  for(i in 1:N_subjects[3]){
    y31[i] ~ multi_normal(mu31,Sigma31);
  }
}

library(mvtnorm)
library(rstan)
library(MCMCpack)
N_subjects <- rep(100,3)
N_items <- 3

S <- riwish(N_items,diag(rep(1,N_items)))

mu <- rnorm(N_items)
tmp <- rmvnorm(N_subjects[1],mean=mu,sigma=S)
y12 <- tmp[,1:2]
tmp <- rmvnorm(N_subjects[2],mean=mu,sigma=S)
y23 <- tmp[,2:3]
tmp <- rmvnorm(N_subjects[3],mean=mu,sigma=S)
y31 <- tmp[,c(3,1)]

dat1 <- list(N_subjects=N_subjects,N_items=N_items,y12=y12,y23=y23,y31=y31)

out1 <- stan(file="multinormal4.stan",data=dat1)

ymp <- rmvnorm(sum(N_subjects)/3*3,mean=mu,sigma=S)
dat2 <- list(N_subjects=N_subjects[3],N_items=N_items,y=tmp)
out2 <- stan(file="multinormal2.stan",data=dat2)