# Vehicle model

Dear colleagues,

from the code you will see that I am a new to stan (pystan) and I run into problems when I want to fit the model (compiling works).

I very often get a message that the matrix for the Choleski decomposition is not positive definite and sure it should be positive definite.

Message: “Informational Message: The current Metropolis proposal is about to be rejected because of the following issue: Exception: cholesky_decompose: Matrix m is not positive definite (in ‘model_03_multivariate_ln.stan’ at line 68)”

I warmly welcome any helping advice.
With kind regards.
Oliver

Here is the model: It is about a vehicle with 4 axe. With axe load data I want to create a stan model.

``````//  MODEL 03
//	4 axes vehicle

data {
int < lower = 1 > data_n;   // number of values
int < lower = 1 > axes_n;   // number of axes

vector[data_n] A1;     // data values of axe loads
vector[data_n] A2;     // data values of axe loads
vector[data_n] A3;     // data values of axe loads
vector[data_n] A4;     // data values of axe loads

// some constraints
real < lower = 0 > prior_mu_mu_Ax;      //  prior for mu of mu_Ax (mean of axe load)
real < lower = 0 > prior_sigma_mu_Ax1;  //  prior for the sigma of mu Ax (diagonal elements)
real < lower = 0 > prior_sigma_mu_Ax2;  //  prior for the sigma of mu Ax ( non diagonal elements)

real < lower = 0 > upper_sigma_Ax1_sqrt; // constraint
real < lower = 0 > upper_variance_value; // constraint
}

parameters {
// rho
real<lower=-1, upper=1> rho_12;
real<lower=-1, upper=1> rho_13;
real<lower=-1, upper=1> rho_14;
real<lower=-1, upper=1> rho_23;
real<lower=-1, upper=1> rho_24;
real<lower=-1, upper=1> rho_34;

vector[axes_n] mu_Ax;

// This is to restrict the covariance matrix for Ax
real<lower=0 , upper=upper_sigma_Ax1_sqrt> sigma_Ax1_sqrt;
real<lower=0 , upper=upper_sigma_Ax1_sqrt> sigma_Ax2_sqrt;
real<lower=0 , upper=upper_sigma_Ax1_sqrt> sigma_Ax3_sqrt;
real<lower=0 , upper=upper_sigma_Ax1_sqrt> sigma_Ax4_sqrt;
}

transformed parameters{
vector [ axes_n] Ax[ data_n]; // data values loads all Axis
vector [ axes_n] mu_mu_Ax; // mean of prior of axes loads

matrix <lower=-upper_variance_value, upper=upper_variance_value> [ axes_n, axes_n ] sigma_Ax;     // covariance of the axle loads
matrix <lower=-upper_variance_value, upper=upper_variance_value> [ axes_n, axes_n ] sigma_mu_Ax;  // covariance of mean of the  axle loads
matrix [ axes_n, axes_n ] L_Ax;                              // Cholesky of the axle loads
matrix [ axes_n, axes_n ] L_mu_Ax;                              // Cholesky of mean of axle loads

sigma_Ax[1,1] = sigma_Ax1_sqrt;
sigma_Ax[1,2] = sqrt(sigma_Ax1_sqrt) * sqrt(sigma_Ax2_sqrt) * rho_12;
sigma_Ax[1,3] = sqrt(sigma_Ax1_sqrt) * sqrt(sigma_Ax3_sqrt) * rho_13;
sigma_Ax[1,4] = sqrt(sigma_Ax1_sqrt) * sqrt(sigma_Ax4_sqrt) * rho_14;

sigma_Ax[2,1] = sigma_Ax[1,2];
sigma_Ax[2,2] = sigma_Ax2_sqrt;
sigma_Ax[2,3] = sqrt(sigma_Ax2_sqrt) * sqrt(sigma_Ax3_sqrt) * rho_23;
sigma_Ax[2,4] = sqrt(sigma_Ax2_sqrt) * sqrt(sigma_Ax4_sqrt) * rho_24;

sigma_Ax[3,1] = sigma_Ax[1,3];
sigma_Ax[3,2] = sigma_Ax[2,3];
sigma_Ax[3,3] = sigma_Ax3_sqrt;
sigma_Ax[3,4] = sqrt(sigma_Ax3_sqrt) * sqrt(sigma_Ax4_sqrt) * rho_34;

sigma_Ax[4,1] = sigma_Ax[1,4];
sigma_Ax[4,2] = sigma_Ax[2,4];
sigma_Ax[4,3] = sigma_Ax[3,4];
sigma_Ax[4,4] = sigma_Ax4_sqrt;

L_Ax = cholesky_decompose( sigma_Ax );

for (i in 1 : data_n){
Ax[i,1] = A1[i];
Ax[i,2] = A2[i];
Ax[i,3] = A3[i];
Ax[i,4] = A4[i];
}

mu_mu_Ax[1] = prior_mu_mu_Ax;
mu_mu_Ax[2] = prior_mu_mu_Ax;
mu_mu_Ax[3] = prior_mu_mu_Ax;
mu_mu_Ax[4] = prior_mu_mu_Ax;

sigma_mu_Ax[1,1] = pow( prior_sigma_mu_Ax1  ,2 ) ;
sigma_mu_Ax[1,2] = pow( prior_sigma_mu_Ax2 ,2 );
sigma_mu_Ax[1,3] = pow( prior_sigma_mu_Ax2 ,2 );
sigma_mu_Ax[1,4] = pow( prior_sigma_mu_Ax2 ,2 );

sigma_mu_Ax[2,1] = pow( prior_sigma_mu_Ax2 ,2 );
sigma_mu_Ax[2,2] = pow( prior_sigma_mu_Ax1  ,2 );
sigma_mu_Ax[2,3] = pow( prior_sigma_mu_Ax2 ,2 );
sigma_mu_Ax[2,4] = pow( prior_sigma_mu_Ax2 ,2 );

sigma_mu_Ax[3,1] = pow( prior_sigma_mu_Ax2 ,2 );
sigma_mu_Ax[3,2] = pow( prior_sigma_mu_Ax2 ,2 );
sigma_mu_Ax[3,3] = pow( prior_sigma_mu_Ax1  ,2 );
sigma_mu_Ax[3,4] = pow( prior_sigma_mu_Ax2 ,2 );

sigma_mu_Ax[4,1] = pow( prior_sigma_mu_Ax2 ,2 );
sigma_mu_Ax[4,2] = pow( prior_sigma_mu_Ax2 ,2 );
sigma_mu_Ax[4,3] = pow( prior_sigma_mu_Ax2 ,2 );
sigma_mu_Ax[4,4] = pow( prior_sigma_mu_Ax1  ,2 );

L_mu_Ax = cholesky_decompose( sigma_mu_Ax );

}

model {
// prior
mu_Ax ~ multi_normal_cholesky_lpdf(mu_mu_Ax, L_mu_Ax);

//
Ax ~ multi_normal_cholesky_lpdf( mu_Ax, L_Ax );
}
``````
1 Like

Unfortunately I cannot delve deeply into the model, but what sometimes happens is that the matrix that is analytically positive-semidefinite by construction is not positive semi-definite in practice due to numerical inaccuracies of the floating point mathematic. A frequently used hack is to add a small number (say 1e-10) to the diagonal of the matrix and if your code is otherwise correct, it usually makes those warnings go away without causing any real difference to the posterior.

Does that help?

Best of luck with your model!

1 Like