Vehicle model

Modeling
1

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
2

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