I encountered the similar kind of numerical issues in this post, i.e., see codes below, basically my covariance matrix was reported to be insymmetric during run execution because of rounding.
vector theta_joint_pred_rng(matrix x, matrix x_pred, real tau, real alpha, real l,
matrix PHI, matrix PHI_pred, vector sqrt_SPD, vector theta_obs){
int n = rows(PHI);
int n_pred = rows(PHI_pred);
int M_nD = cols(PHI);
vector[n_pred] theta_pred;
{
matrix[n,M_nD] PHI_sqrt_SPD = PHI .* rep_matrix(sqrt_SPD',n);
matrix[n_pred,M_nD] PHI_pred_sqrt_SPD = PHI_pred .* rep_matrix(sqrt_SPD',n_pred);
matrix[n_pred,n] K21 = PHI_pred_sqrt_SPD * PHI_sqrt_SPD';
matrix[n,n] K11_inv= add_diag(-PHI * inverse_spd(add_diag(PHI'*PHI, tau^2/sqrt_SPD^2))*PHI', 1)/tau^2;
matrix[n_pred,n_pred] K22 = alpha^2 * add_diag(PHI_pred_sqrt_SPD * PHI_pred_sqrt_SPD', tau^2);
vector[n_pred] theta_pred_mu = K21 * K11_inv * theta_obs;
matrix[n_pred,n_pred] theta_pred_cov = K22 - K21 * K11_inv * K21';
theta_pred = multi_normal_rng(theta_pred_mu, theta_pred_cov);
}
return theta_pred;
}
I have two questions.
- I wonder whether the easiest way to resolve this issue is to use, for example:
theta_pred_cov = symmetrize_from_lower_tri(theta_pred_cov);
- I’ve tried to define these matrices as covariance matrix and was wondering whether that might help enforce the symmetry. However, not sure why, I couldn’t compile my codes when I tried that. See codes below:
cov_matrix[n] K11_inv;
K11_inv = add_diag(-PHI * inverse_spd(add_diag(PHI'*PHI, tau^2/sqrt_SPD^2))*PHI', 1)/tau^2;
I got the following error while trying to compare these changed codes. I wonder whether anyone knows why?
Variable declaration, statement or "}" expected.