Vectorization of multi_normal, multi_normal_prec and multi_normal_cholesky + tests

Author: randommm

basic_estimators/normal_multi.data.R

@@ -0,0 +1,13 @@
+N <- 10
+J <- 5
+sigma <-
+structure(c(1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 
+0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 
+0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 
+0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 
+0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1), .Dim = c(10L, 
+10L))
+y <- structure(c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 
+16, 17, 1, 19, 20, 21, 22, 23, 24, 25, 26, 3, 28, 29, 30, 31, 
+32, 33, 3, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 
+48, 49, 50), .Dim = c(5L, 10L))

basic_estimators/normal_multi.stan

@@ -0,0 +1,52 @@
+data {
+  int<lower=0> N; //response variable size
+  int<lower=0> J; //number of observations
+  vector[N] y[J];
+  matrix[N, N] sigma;
+}
+transformed data {
+  row_vector[N] ry[J];
+  matrix[N, N] inv_sigma;
+  for (n in 1:N)
+    for (j in 1:J)
+      ry[j][n] <- y[j][n];
+  inv_sigma <- inverse_spd(sigma);
+}
+parameters {
+  vector[N] beta[J];
+  row_vector[N] rbeta[J];
+}
+model {
+  for (n in 1:J)
+    ry[n] ~ multi_normal(beta[n], sigma);
+  y ~ multi_normal(beta, sigma);
+  ry ~ multi_normal(beta[1], sigma);
+  y[1] ~ multi_normal(beta, sigma);
+  for (n in 1:J)
+    y[n] ~ multi_normal(rbeta[n], sigma);
+  y ~ multi_normal(rbeta, sigma);
+  y ~ multi_normal(rbeta[1], sigma);
+  ry[1] ~ multi_normal(rbeta, sigma);
+  
+  for (n in 1:J)
+    ry[n] ~ multi_normal_cholesky(beta[n], sigma);
+  y ~ multi_normal_cholesky(beta, sigma);
+  ry ~ multi_normal_cholesky(beta[1], sigma);
+  y[1] ~ multi_normal_cholesky(beta, sigma);
+  for (n in 1:J)
+    y[n] ~ multi_normal_cholesky(rbeta[n], sigma);
+  y ~ multi_normal_cholesky(rbeta, sigma);
+  y ~ multi_normal_cholesky(rbeta[1], sigma);
+  ry[1] ~ multi_normal_cholesky(rbeta, sigma);
+  
+  for (n in 1:J)
+    ry[n] ~ multi_normal_prec(beta[n], inv_sigma);
+  y ~ multi_normal_prec(beta, inv_sigma);
+  ry ~ multi_normal_prec(beta[1], inv_sigma);
+  y[1] ~ multi_normal_prec(beta, inv_sigma);
+  for (n in 1:J)
+    y[n] ~ multi_normal_prec(rbeta[n], inv_sigma);
+  y ~ multi_normal_prec(rbeta, inv_sigma);
+  y ~ multi_normal_prec(rbeta[1], inv_sigma);
+  ry[1] ~ multi_normal_prec(rbeta, inv_sigma);
+}