fixed coupla_defination.stan of zero_constrain

examples/linear_algebra/R/correlation_angles_example.R

@@ -1,12 +1,15 @@
 library(cmdstanr)
 fp <- file.path("./examples/linear_algebra/stan/correlation_angles_example.stan")
-mod <- cmdstan_model(fp, include_paths = "./functions/linear_algebra")
+mod <- cmdstan_model(fp, include_paths = "./functions/linear_algebra", force_recompile = T)
 
 T <- matrix(c(0.8, -0.9, -0.9,
               -0.9, 1.1, 0.3,
               -0.9, 0.4, 0.9), 3, 3)
 
-
+T <- matrix(c(1, -0.9, -0.9,
+              -0.9, 1, 0.3,
+              -0.9, 0.4, 1), 3, 3)
+chol(T)
 mod_out <- mod$optimize(
   data = list(N = 3,
               R = T)
@@ -17,12 +20,10 @@ mod_out <- mod$sample(
               R = T),
   chains = 2,
   seed = 23421,
+
   parallel_chains = 2,
   iter_warmup = 600,
   iter_sampling = 600
 )
 
-mod_out$summary("R_hat")
-mat <- matrix(mod_out$summary("R_hat")$mean, 3, 3)
-mat
-chol(mat)
+mod_out$summary()

examples/linear_algebra/R/correlation_angles_example2.R

@@ -0,0 +1,120 @@
+library(cmdstanr)
+fp <- file.path("./examples/linear_algebra/stan/correlation_angles_constrain_example.stan")
+mod <- cmdstan_model(fp, include_paths = "./functions/linear_algebra", force_recompile = T)
+
+library(rstan)
+expose_stan_functions(fp)
+
+library(StanHeaders)
+stanFunction("multiply_lower_tri_self_transpose", x = L)
+
+
+N <- 3
+test_mat <- matrix(c(0.8, -0.9, -0.9,
+              -0.9, 1.1, 0.3,
+              -0.9, 0.4, 0.9), 3, 3)
+
+T <- matrix(c(1, -0.9, -0.9,
+              -0.9, 1, 0.3,
+              -0.9, 0.4, 1), 3, 3)
+
+
+test_mat <- matrix(c(1.0000,   0,         0,         0,   -0.9360,
+                     0,    1.0000,   -0.5500,   -0.3645,   -0.5300,
+                     0,   -0.5500,    1.0000,   0,    0.0875,
+                     0,   -0.3645,    0,    1.0000,    0.4557,
+                     -0.9360,   -0.5300,    0.0875,    0.4557,    1.0000),
+                   byrow = T, 5, 5)
+
+where_zero <- (test_mat[upper.tri(test_mat)] == 0) * 1
+sum(where_zero) * 1
+
+fit <- stan(file = fp, 
+            data = list(N = nrow(test_mat),
+                        R = test_mat,
+                        is_symmetric = 1,
+                        Z = sum(where_zero),
+                        K = ( nrow(test_mat) * (nrow(test_mat) - 1 )) / 2,
+                        where_zero = where_zero ),
+            chains = 1)
+
+K <- 5
+theta_test <- build_angle_mat(where_zero, runif(6, 0, pi), K)
+L = zero_constrain(theta_test, K)
+tcrossprod(L)
+
+angle_raw <-  runif(6)
+N <- length(where_zero);
+mat <- matrix(0, N, N)
+count <- 1;
+raw_count <- 1;
+
+for (k in 2:K){
+  mat[k, k] = 0;
+  for (j in 1:k - 1) {
+    if (where_zero[raw_count] != 1) {
+      mat[k, j] = angle_raw[count];
+      count <= count +  1;
+    }
+    raw_count = raw_count + 1;
+  }
+}
+
+count <- 0
+
+for (k in 2:K){
+  for (j in 1:(k - 1)) {
+    count <- count + 1
+    print(count)
+  }}
+
+
+
+mod_out <- mod$sample(
+  data = list(N = nrow(test_mat),
+              R = test_mat,
+              is_symmetric = 1,
+              Z = sum(where_zero),
+              K = ( nrow(test_mat) * (nrow(test_mat) - 1 )) / 2,
+              where_zero = where_zero ),
+  chains = 2,
+  #init = 0,
+  #seed = 23421,
+  #adapt_delta = 0.8,
+  #max_treedepth = 10,
+  parallel_chains = 4,
+  iter_warmup = 400,
+  iter_sampling = 400
+)
+
+N <- nrow(test_mat)
+
+round(matrix(mod_out$summary("R_out")$mean,nrow(test_mat), nrow(test_mat)), 3)
+
+
+chol(test_mat)
+
+mod_out <- mod$optimize(
+  data = list(N = nrow(test_mat),
+                      R = test_mat,
+                      is_symmetric = 1,
+                      Z = sum(where_zero),
+                      K = ( nrow(test_mat) * (nrow(test_mat) - 1 )) / 2,
+                      where_zero = where_zero ))
+
+round(matrix(mod_out$summary("R_out")$estimate,nrow(test_mat), nrow(test_mat)), 3)
+
+mod_out <- mod$sample(
+  data = list(N = nrow(test_mat),
+              R = test_mat,
+              is_symmetric = 1),
+  chains = 2,
+  seed = 23421,
+  adapt_delta = 0.9,
+  max_treedepth = 10,
+  parallel_chains = 2,
+  iter_warmup = 200,
+  iter_sampling = 200
+)
+mod_out$summary("R_out")
+round(matrix(mod_out$summary("R_out")$mean, N, N), 3)

examples/linear_algebra/stan/correlation_angles_constrain_example.stan

@@ -0,0 +1,124 @@
+functions {
+ vector lower_elements(matrix M, int tri_size){
+    int n = rows(M);
+    int counter = 1;
+    vector[tri_size] out;
+
+    for (i in 2:n){
+      for (j in 1:(i - 1)) {
+        out[counter] = M[i, j];
+        counter += 1;
+      }
+    }
+    return out;
+  }
+
+   vector lower_elements_constrain(matrix M, int A){
+    int n = rows(M);
+    int counter = 1;
+    vector[A] out;
+
+    for (i in 2:n){
+      for (j in 1:(i - 1)) {
+        if(M[i, j] > 0){
+         out[counter] = M[i, j];
+         counter += 1;
+        }
+      }
+    }
+    return out;
+  }
+
+ matrix build_angle_mat (vector where_zero, vector angle_raw, int K) {
+  int N = num_elements(where_zero);
+  matrix[K, K] mat;
+  int count = 1;
+  int raw_count = 0;
+
+  mat[1, 1] = 0.0;
+  for (k in 2:K){
+    mat[k, k] = 0.0;
+    for (j in 1:k - 1) {
+      raw_count += 1;
+      if (where_zero[raw_count] != 1) {
+         mat[k, j] = angle_raw[count];
+         count += 1;
+      }
+    }
+  }
+
+  return mat;
+}
+
+ matrix zero_constrain (matrix angle_raw, int N){
+  matrix[N, N] inv_mat;
+  matrix[N, N] angle = angle_raw;
+
+  inv_mat[1, 1] = 1;
+
+  for (i in 2:N) {
+    // constrain first column
+    // if C = BB^T then for the first column
+    // c_{i, 1} = b_{i, 1} = 0 so 
+    // if c_{i, 1} = 0 then cos(\theta_{i, 1}) = 0
+    // acos(0) = \theta_{i, 1} = pi / 2
+     if (is_nan(angle_raw[i, 1]) == 1) {
+       inv_mat[i, 1] = 0;
+       angle[i, 1] = pi() / 2;
+     } else inv_mat[i, 1] = cos(angle[i, 1]);
+
+    if (i > 2) {
+      for (j in 2:(i - 1)) {
+        real prod_sines = prod(sin(angle[i, 1:j - 1]));
+        real cos_theta;
+        if (is_nan(angle_raw[i, j]) == 1) {
+          cos_theta = -(dot_product(inv_mat[j, 1:j - 1], inv_mat[i, 1:j - 1]) ) / ( inv_mat[j, j] * prod_sines );
+            if ( cos_theta < -1 || cos_theta > 1 ) reject("cos_theta is ", cos_theta, " and must be in [-1, 1]"); // cos_theta = 0;
+           // else if( cos_theta > 1) cos_theta = 1;
+            angle[i, j] = acos( cos_theta );
+          }
+        inv_mat[i, j] = cos(angle[i, j]) * prod_sines;
+      }
+    }
+    inv_mat[i, i] = prod(sin(angle[i, 1:(i - 1)]));
+  }
+  return inv_mat;
+}
+}
+data {
+  int N;
+  matrix[N, N] R;
+  int<lower=0, upper=1> is_symmetric;
+  int Z; // number of zeros
+  int K;
+  vector[K] where_zero;
+}
+transformed data {
+ // int K = N * (N - 1) / 2;
+  row_vector[N] m;
+
+  for (n in 1:N)
+    m[n] = N - n + 1;
+}
+parameters {
+  vector<lower = 0, upper = pi()>[K - Z] theta_free;
+} 
+transformed parameters {
+  matrix[N, N] theta_mat = build_angle_mat(where_zero, theta_free, N);
+  matrix[N, N] L = zero_constrain(theta_mat, N);
+  matrix[N, N] R_hat = multiply_lower_tri_self_transpose(L);
+}
+model {
+  // log_absd_jacs 
+  // sin(theta) is always > 0 since theta in (0, pi)
+  // because of this the diagonal of L is always > 0 also
+   target += 2 * sum(log(sin(lower_elements_constrain(theta_mat, K - Z))));            // angle2chol
+   target += N * log(2) + m * diagonal(log(L));   // cholesky tcrossprod
+   if (is_symmetric == 1)
+    lower_elements(R, K) ~ normal(lower_elements(R_hat, K), 0.001);
+    else to_vector(R) ~ normal(to_vector(R_hat), 0.001);
+}
+generated quantities {
+   matrix[N, N] theta_out = build_angle_mat(where_zero, theta_free, N);
+   matrix[N, N] R_out = multiply_lower_tri_self_transpose( zero_constrain(theta_out, N) );
+}