fixed coupla_defination.stan of zero_constrain #69

yamikarajput546 · 2022-01-29T06:46:10Z

Please check the PR and merge the changes

spinkney

Let's just make correlation_angles_constrain_working to correlation_angles_constrain.stan

spinkney · 2022-01-29T11:45:18Z

examples/linear_algebra/stan/correlation_angles_constrain_example.stan

+ 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;
+}


Add these with #include

yes i've update these changes to constrain_example.stan and make the correlation_angles_constrain_working to correlation_angles_constrain.stan

spinkney · 2022-01-29T11:47:48Z

examples/linear_algebra/stan/correlation_angles_constrain_working.stan

+ vector lower_elements(matrix M){
+    int n = rows(M);
+    int k = n * (n - 1) / 2;
+    int counter = 1;
+    vector[k] 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;
+}
+
+ // set the upper triangle to 0
+ // only looking at strictly lower tri part
+ vector sparse_cholesky_lp (vector angle_raw, int[] csr_rows, int[] csr_cols, int Z, int N){
+  vector[Z + N] sparse_chol; // Z + N is the number of non-zero values in lower tri plus
+                             // the diagonal since the angles do not include the diagonal
+  int R = size(csr_rows);
+  int C = size(csr_cols);
+  int S[R, C] = append_array(csr_rows, csr_cols);
+  int count = 1;
+
+  sparse_chol[count] = 1;
+
+  // traversing in col-major order
+  // S[2, 1]           skips S[2, 2]
+  // S[3, 1], S[3, 2]  skips S[3, 3]
+  // etc
+  for (r in S) {
+    int this_rows_column_num = 0; 
+    for (c in r) {
+      count += 1;
+      this_rows_column_num += 1;
+
+      if(this_rows_column_num == 1)
+        sparse_chol[count] = cos(angle_raw[count]); // constrain first column
+
+
+    }
+
+    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;
+}


use #include to add these

included and update the changes to working.stan

spinkney · 2022-01-29T11:48:22Z

functions/copula/copula_definition.stan

+/** \defgroup copula Copula Functions
+  *
+  * A copula is a multivariate cumulative distribution function for which the 
+  * marginal probability distribution of each variable is uniform on the interval \f$[0, 1]\f$.
+  * Copulas are used to describe the dependence between random variables.
+  *
+  **/
+
+


Remove this change

yeah i removed it

spinkney · 2022-01-29T11:49:26Z

functions/linear_algebra/correlation_angles.stan

@@ -12,7 +12,7 @@
   */
 matrix angle2chol (matrix angle_mat){
  int N = rows(angle_mat);
-  matrix[N, N] inv_mat = identity_matrix(N);
+  matrix[N, N] inv_mat;


Let's add the doxygen style format.

Also include yourself under the copyright!

spinkney · 2022-01-29T11:50:31Z

functions/linear_algebra/triangular.stan

-  }
-
-  return out;
+  vector lower_elements(matrix M, int tri_size){


Does the previous version not work? I have some vague remembrance of a potential issue here.

spinkney

Let's just make correlation_angles_constrain_working to correlation_angles_constrain.stan

spinkney · 2022-01-29T11:51:06Z

@yamikarajput546 thanks for these! Welcome and look forward to working with you

yamikarajput546 · 2022-01-29T13:37:23Z

i've update the changes but i don't know why again the conflicts ?
@spinkney please review it

…46/helpful_stan_functions into yamikarajput546-fix_zero_constrain

spinkney

I merged the other pr 🎉! And I fixed the merge conflicts with this one.

For this one:

add those functions to separate .stanfunctions files and document
we need the doxygen comments at the top and a description of the math. Look at https://github.com/spinkney/helpful_stan_functions/blob/main/functions/copula/normal_copula.stanfunctions how to do the LaTeX and documentation.
rename the file correlation_angles_constrain_work.stan to just correlation_angles_constrain.stanfunctions

spinkney · 2022-01-30T11:10:16Z

examples/linear_algebra/stan/correlation_angles_constrain.stan

+ vector lower_elements(matrix M){
+    int n = rows(M);
+    int k = n * (n - 1) / 2;
+    int counter = 1;
+    vector[k] 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;
+}
+
+ // set the upper triangle to 0
+ // only looking at strictly lower tri part
+ vector sparse_cholesky_lp (vector angle_raw, int[] csr_rows, int[] csr_cols, int Z, int N){
+  vector[Z + N] sparse_chol; // Z + N is the number of non-zero values in lower tri plus
+                             // the diagonal since the angles do not include the diagonal
+  int R = size(csr_rows);
+  int C = size(csr_cols);
+  int S[R, C] = append_array(csr_rows, csr_cols);
+  int count = 1;
+
+  sparse_chol[count] = 1;
+
+  // traversing in col-major order
+  // S[2, 1]           skips S[2, 2]
+  // S[3, 1], S[3, 2]  skips S[3, 3]
+  // etc
+  for (r in S) {
+    int this_rows_column_num = 0; 
+    for (c in r) {
+      count += 1;
+      this_rows_column_num += 1;
+
+      if(this_rows_column_num == 1)
+        sparse_chol[count] = cos(angle_raw[count]); // constrain first column
+
+
+    }
+
+    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;
+}


this is more work but will be nice to have. Can you make each of these into their own file and then reference here?

no issue @spinkney i'll surely update all changes :)

spinkney · 2022-01-30T11:11:00Z

examples/linear_algebra/stan/correlation_angles_constrain_example.stan

+ 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;


again, add the files and reference here

spinkney · 2022-01-30T11:11:21Z

examples/linear_algebra/stan/correlation_angles_constrain_working.stan

+ vector lower_elements(matrix M){
+    int n = rows(M);
+    int k = n * (n - 1) / 2;
+    int counter = 1;
+    vector[k] 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;
+}
+
+ // set the upper triangle to 0
+ // only looking at strictly lower tri part
+ vector sparse_cholesky_lp (vector angle_raw, int[] csr_rows, int[] csr_cols, int Z, int N){
+  vector[Z + N] sparse_chol; // Z + N is the number of non-zero values in lower tri plus
+                             // the diagonal since the angles do not include the diagonal
+  int R = size(csr_rows);
+  int C = size(csr_cols);
+  int S[R, C] = append_array(csr_rows, csr_cols);
+  int count = 1;
+
+  sparse_chol[count] = 1;
+
+  // traversing in col-major order
+  // S[2, 1]           skips S[2, 2]
+  // S[3, 1], S[3, 2]  skips S[3, 3]
+  // etc
+  for (r in S) {
+    int this_rows_column_num = 0; 
+    for (c in r) {
+      count += 1;
+      this_rows_column_num += 1;
+
+      if(this_rows_column_num == 1)
+        sparse_chol[count] = cos(angle_raw[count]); // constrain first column
+
+
+    }
+
+    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;


add files and reference

yamikarajput546 · 2022-01-30T11:35:48Z

I merged the other pr tada! And I fixed the merge conflicts with this one.

For this one:

add those functions to separate .stanfunctions files and document

we need the doxygen comments at the top and a description of the math. Look at https://github.com/spinkney/helpful_stan_functions/blob/main/functions/copula/normal_copula.stanfunctions how to do the LaTeX and documentation.

rename the file correlation_angles_constrain_work.stan to just correlation_angles_constrain.stanfunctions

Thanks for that @spinkney i'll update as you mentioned like renamed and pushed the changes but as the workflow has failed so will try the doxygen comments

spinkney · 2022-02-01T11:45:51Z

I merged the other pr tada! And I fixed the merge conflicts with this one.
For this one:

add those functions to separate .stanfunctions files and document

we need the doxygen comments at the top and a description of the math. Look at https://github.com/spinkney/helpful_stan_functions/blob/main/functions/copula/normal_copula.stanfunctions how to do the LaTeX and documentation.

rename the file correlation_angles_constrain_work.stan to just correlation_angles_constrain.stanfunctions

Thanks for that @spinkney i'll update as you mentioned like renamed and pushed the changes but as the workflow has failed so will try the doxygen comments

Just tag me once you want a review. Thanks again!

spinkney and others added 11 commits February 11, 2021 13:32

added

8a3f550

works but had to comment out the jacobian adjustments

5277812

remove unnecessary jitter

0f1d02a

added reject statement

ebb52af

saving unsaved changes

7ce1ff1

fixed jacobian

3a433bc

Merge branch 'main' into zero_constrain

2e60ea3

working on sparse matrix for zeros

790c2ab

Merge branch 'main' into zero_constrain

fd76b68

saving progress

211c60a

fixed coupla_defination.stan of zero_constrain

2084f0b

spinkney reviewed Jan 29, 2022

View reviewed changes

yamikarajput546 added 2 commits January 29, 2022 18:50

resolve conflicts

5d7b582

fix changes

f4bec14

Merge branch 'fix_zero_constrain' of https://github.com/yamikarajput5…

29c40d7

…46/helpful_stan_functions into yamikarajput546-fix_zero_constrain

spinkney reviewed Jan 30, 2022

View reviewed changes

fix changes

ccdd0dd

spinkney force-pushed the main branch from e0fa60d to e70486d Compare March 12, 2022 13:19

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

fixed coupla_defination.stan of zero_constrain #69

fixed coupla_defination.stan of zero_constrain #69

yamikarajput546 commented Jan 29, 2022

spinkney left a comment

spinkney Jan 29, 2022

yamikarajput546 Jan 29, 2022

spinkney Jan 29, 2022

yamikarajput546 Jan 29, 2022

spinkney Jan 29, 2022

yamikarajput546 Jan 29, 2022

spinkney Jan 29, 2022

spinkney Jan 29, 2022

spinkney left a comment

spinkney commented Jan 29, 2022

yamikarajput546 commented Jan 29, 2022

spinkney left a comment •

edited

Loading

spinkney Jan 30, 2022

yamikarajput546 Jan 30, 2022

spinkney Jan 30, 2022

spinkney Jan 30, 2022

yamikarajput546 commented Jan 30, 2022 •

edited

Loading

spinkney commented Feb 1, 2022

fixed coupla_defination.stan of zero_constrain #69

Are you sure you want to change the base?

fixed coupla_defination.stan of zero_constrain #69

Conversation

yamikarajput546 commented Jan 29, 2022

spinkney left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

spinkney left a comment

Choose a reason for hiding this comment

spinkney commented Jan 29, 2022

yamikarajput546 commented Jan 29, 2022

spinkney left a comment • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

yamikarajput546 commented Jan 30, 2022 • edited Loading

spinkney commented Feb 1, 2022

spinkney left a comment •

edited

Loading

yamikarajput546 commented Jan 30, 2022 •

edited

Loading