diff --git a/functions/copula/centered_gaussian_copula.stanfunctions b/functions/copula/centered_gaussian_copula.stanfunctions
index 8151e6f..550f4d9 100644
--- a/functions/copula/centered_gaussian_copula.stanfunctions
+++ b/functions/copula/centered_gaussian_copula.stanfunctions
@@ -31,12 +31,12 @@ array[] matrix normal_marginal(matrix y, matrix mu_glm, vector sigma) {
   array[2] matrix[N, J] rtn;
   // Initialise the jacobian adjustments to zero, as vectorised lpdf will be used
   rtn[2] = rep_matrix(0, N, J);
-
-  for (j in 1:J) {
-      rtn[1][ : , j] = (y[: , j] - mu_glm[ : , j]) / sigma[j];
-      rtn[2][1, j] = normal_lpdf(y[ : , j] | mu_glm[ : , j], sigma[j]);
+  
+  for (j in 1 : J) {
+    rtn[1][ : , j] = (y[ : , j] - mu_glm[ : , j]) / sigma[j];
+    rtn[2][1, j] = normal_lpdf(y[ : , j] | mu_glm[ : , j], sigma[j]);
   }
-
+  
   return rtn;
 }
 
@@ -63,11 +63,11 @@ array[] matrix bernoulli_marginal(array[,] int y, matrix mu_glm, matrix u_raw) {
   int J = cols(mu_glm);
   matrix[N, J] matrix_y = to_matrix(y);
   matrix[N, J] mu_glm_logit = 1 - inv_logit(mu_glm);
-
+  
   matrix[N, J] Lbound = matrix_y .* mu_glm_logit;
   matrix[N, J] UmL = fabs(matrix_y - mu_glm_logit);
-
-  return { inv_Phi(Lbound + UmL .* u_raw), log(UmL) };
+  
+  return {inv_Phi(Lbound + UmL .* u_raw), log(UmL)};
 }
 
 /**
@@ -89,15 +89,15 @@ array[] matrix bernoulli_marginal(array[,] int y, matrix mu_glm, matrix u_raw) {
  * @return 2D array of matrices containing the random variables
  *         and jacobian adjustments
  */
-array[] matrix binomial_marginal(array[,] int num, array[,] int den,
-                                 matrix mu_glm, matrix u_raw) {
+array[] matrix binomial_marginal(array[,] int num, array[,] int den, matrix mu_glm,
+                                 matrix u_raw) {
   int N = rows(mu_glm);
   int J = cols(mu_glm);
   matrix[N, J] mu_glm_logit = inv_logit(mu_glm);
   array[2] matrix[N, J] rtn;
-
-  for (j in 1:J) {
-    for (n in 1:N) {
+  
+  for (j in 1 : J) {
+    for (n in 1 : N) {
       real Ubound = binomial_cdf(num[n, j] | den[n, j], mu_glm_logit[n, j]);
       real Lbound = 0;
       if (num[n, j] > 0) {
@@ -108,7 +108,7 @@ array[] matrix binomial_marginal(array[,] int num, array[,] int den,
       rtn[2][n, j] = log(UmL);
     }
   }
-
+  
   return rtn;
 }
 
@@ -135,9 +135,9 @@ array[] matrix poisson_marginal(array[,] int y, matrix mu_glm, matrix u_raw) {
   int J = cols(mu_glm);
   matrix[N, J] mu_glm_exp = exp(mu_glm);
   array[2] matrix[N, J] rtn;
-
-  for (j in 1:J) {
-    for (n in 1:N) {
+  
+  for (j in 1 : J) {
+    for (n in 1 : N) {
       real Ubound = poisson_cdf(y[n, j] | mu_glm_exp[n, j]);
       real Lbound = 0;
       if (y[n, j] > 0) {
@@ -148,11 +148,10 @@ array[] matrix poisson_marginal(array[,] int y, matrix mu_glm, matrix u_raw) {
       rtn[2][n, j] = log(UmL);
     }
   }
-
+  
   return rtn;
 }
 
-
 /**
  * Mixed Copula Log-Probability Function
  *
@@ -167,23 +166,23 @@ real centered_gaussian_copula_cholesky_lpdf(array[,] matrix marginals, matrix L)
   int N = rows(marginals[1][1]);
   int J = rows(L);
   matrix[N, J] U;
-
+  
   // Iterate through marginal arrays, concatenating the outcome matrices by column
   // and aggregating the log-likelihoods (from continuous marginals) and jacobian
   // adjustments (from discrete marginals)
   real adj = 0;
   int pos = 1;
-  for (m in 1:size(marginals)) {
+  for (m in 1 : size(marginals)) {
     int curr_cols = cols(marginals[m][1]);
-
-    U[ : , pos:(pos + curr_cols - 1)] = marginals[m][1];
-
+    
+    U[ : , pos : (pos + curr_cols - 1)] = marginals[m][1];
+    
     adj += sum(marginals[m][2]);
     pos += curr_cols;
   }
-
+  
   // Return the sum of the log-probability for copula outcomes and jacobian adjustments
   return multi_normal_cholesky_copula_lpdf(U | L) + adj;
 }
 
-/** @} */
+/** @} */
\ No newline at end of file
diff --git a/functions/copula/normal_copula.stanfunctions b/functions/copula/normal_copula.stanfunctions
index 4a2a481..31bd81e 100644
--- a/functions/copula/normal_copula.stanfunctions
+++ b/functions/copula/normal_copula.stanfunctions
@@ -30,7 +30,7 @@
   */
 real normal_copula_lpdf(real u, real v, real rho) {
   real rho_sq = square(rho);
-
+  
   return (0.5 * rho * (-2. * u * v + square(u) * rho + square(v) * rho)) / (-1. + rho_sq)
          - 0.5 * log1m(rho_sq);
 }
@@ -52,12 +52,12 @@ real normal_copula_lpdf(real u, real v, real rho) {
 real normal_copula_vector_lpdf(vector u, vector v, real rho) {
   int N = num_elements(u);
   real rho_sq = square(rho);
-
+  
   real a1 = 0.5 * rho;
   real a2 = rho_sq - 1;
   real a3 = 0.5 * log1m(rho_sq);
   real x = -2 * u' * v + rho * (dot_self(u) + dot_self(v));
-
+  
   return a1 * x / a2 - N * a3;
 }
 
@@ -97,7 +97,7 @@ real multi_normal_copula_lpdf(matrix u, matrix L) {
   int N = cols(u);
   real inv_sqrt_det_log = sum(log(diagonal(L)));
   matrix[K, N] x = mdivide_left_tri_low(L, u);
-
+  
   return -N * inv_sqrt_det_log - 0.5 * sum(columns_dot_self(x) - columns_dot_self(u));
 }
 
@@ -122,10 +122,10 @@ real bivariate_normal_copula_cdf(vector u, real rho) {
   real alpha_u = a * (pv / pu - rho);
   real alpha_v = a * (pu / pv - rho);
   real d = 0;
-
-  if (u[1] < 0.5 && u[2] >= 0.5 || u[1] >= 0.5 && u[2] < 0.5)
+  
+  if (u[1] < 0.5 && u[2] >= 0.5 || u[1] >= 0.5 && u[2] < 0.5) 
     d = 0.5;
-
+  
   return avg_uv - owens_t(pu, alpha_u) - owens_t(pv, alpha_v) - d;
 }
 
@@ -165,7 +165,6 @@ real bivariate_normal_copula_cdf(vector u, real rho) {
   * @return log density of the distribution
   */
 
-
 /**
   * Multivariate Normal Copula lpdf (Cholesky parameterisation)
   *
@@ -176,7 +175,7 @@ real bivariate_normal_copula_cdf(vector u, real rho) {
 real multi_normal_cholesky_copula_lpdf(matrix U, matrix L) {
   int N = rows(U);
   int J = cols(U);
-  matrix[J,J] Gammainv = chol2inv(L);
+  matrix[J, J] Gammainv = chol2inv(L);
   return -N * sum(log(diagonal(L))) - 0.5 * sum(add_diag(Gammainv, -1.0) .* crossprod(U));
 }
-/** @} */
+/** @} */
\ No newline at end of file