Non converging model, which should be reparameterized

[edit: escape code]

I am modelling 4 variables A_i~N, which in sum give Z=sum(A_i) and Z follows a generalized normal distribution.
The data shows that the standard deviations D[A_i(Z)] are a function of Z and rho between A rho_ij(Z) as well. The code compiles and fits parameters however, stan throws a warning that the model should be reparametrized.

“WARNING:pystan:E-BFMI below 0.2 indicates you may need to reparameterize your model”

I am grateful for any advice and thank you in advance!

// MODEL 104

functions{
        real func_linear(real z, real zmin, real zmax, real rho_zmin,  real rho_zmax ){
            //  linear interpolation between zmin and zmax
            real y;
            y =  rho_zmin + (z- zmin) * (rho_zmax - rho_zmin) / (zmax - zmin)   ;
            return y;
        }

        real gen_normal_lpdf( vector x, real mu_, real alpha_, real beta_ ) {
            vector[rows(x)] lp;
            int N;

            N = rows(x);
            for(i in 1:N){
               // pdf according to wikipedia
               // https://en.wikipedia.org/wiki/Generalized_normal_distribution
               lp[ i ] = log(beta_) - (log(2*alpha_)+lgamma(1/beta_)) - pow( fabs(x[i]-mu_)/alpha_, beta_ );
            }
            // returns loglikelihood to maximise it (not neg. likelihood!)
            return sum( lp );
        }
}

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

    matrix[data_n, axes_n] A;   // data values of axe loads
    vector[data_n] Z;

    real < lower = 0, upper = 35 > prior_A_mean;          // mean of axe load for prior distribution
    real < lower = 0, upper = 10 > prior_A_std;           // std dev of axe load for prior distribution
    real < lower = 0, upper = 1 > prior_A_cov;           // CoV for the priors

    real  prior_rho_c12;           //
    real  prior_rho_c13;           //
    real  prior_rho_c14;           //
    real  prior_rho_c23;           //
    real  prior_rho_c24;           //
    real  prior_rho_c34;           //

    real  prior_rho_d12;           //
    real  prior_rho_d13;           //
    real  prior_rho_d14;           //
    real  prior_rho_d23;           //
    real  prior_rho_d24;           //
    real  prior_rho_d34;           //

    real < lower = 0 > Z_0;      // beta parameter for the Z variable
    real < lower = 0 > Z_1;      // beta parameter for the Z variable
    real < lower = 0 > Z_beta;      // beta parameter for the Z variable
}

transformed data {
    vector <lower=0, upper = 35> [data_n] A1;
    vector <lower=0, upper = 35> [data_n] A2;
    vector <lower=0, upper = 35> [data_n] A3;
    vector <lower=0, upper = 35> [data_n] A4;

    A1 = A[:, 1];
    A2 = A[:, 2];
    A3 = A[:, 3];
    A4 = A[:, 4];
}

parameters {
    real < lower = 0, upper = 30 > A1_mean;               //  mean of axe load
    real < lower = 0, upper = 30 > A2_mean;               //  mean of axe load
    real < lower = 0, upper = 30 > A3_mean;               //  mean of axe load
    real < lower = 0, upper = 30 > A4_mean;               //  mean of axe load

    real < lower = -1, upper = 1 > rho_c12;               //
    real < lower = -1, upper = 1 > rho_c13;               //
    real < lower = -1, upper = 1 > rho_c14;               //
    real < lower = -1, upper = 1 > rho_c23;               //
    real < lower = -1, upper = 1 > rho_c24;               //
    real < lower = -1, upper = 1 > rho_c34;               //

    real < lower = -1, upper = 1 > rho_d12;               //
    real < lower = -1, upper = 1 > rho_d13;               //
    real < lower = -1, upper = 1 > rho_d14;               //
    real < lower = -1, upper = 1 > rho_d23;               //
    real < lower = -1, upper = 1 > rho_d24;               //
    real < lower = -1, upper = 1 > rho_d34;               //

    real < lower = 0 > A1_std_c;               //
    real < lower = 0 > A2_std_c;               //

    real < lower = 0 > A1_std_d;               //
    real < lower = 0 > A2_std_d;               //

    real < lower = 0, upper = axes_n * 20 >     Z_alpha;  //  parameter for Z
}

transformed parameters{
    matrix[axes_n, axes_n] rho [data_n] ;                                      //  correlation coefficients
    matrix < lower = -1000, upper = 1000 >[axes_n,axes_n] COVMAT [data_n]  ;  //  covariance matrix
    cholesky_factor_cov [axes_n] L_COVMAT[data_n] ;                         //  choleski of COVMAT

    vector <lower=0, upper = 35> [axes_n]  A_mean;          //  vector of mean values

    real < lower = 0, upper = axes_n * 35  >  Z_mean;       //  mean of total axe load
    real < lower = 0, upper = axes_n * 200 >  Z_var;        //  variance of total axe load

    vector < lower = 0, upper = 20 > [data_n] A1_std;   //  std dev of axe load
    vector < lower = 0, upper = 20 > [data_n] A2_std;   //  std dev of axe load
    vector < lower = 0, upper = 20 > [data_n] A3_std;   //  std dev of axe load
    vector < lower = 0, upper = 20 > [data_n] A4_std;   //  std dev of axe load

    real   p;       //  parameter of the quadratic equation
    real   q;       //  parameter of the quadratic equation

    A_mean[1]= A1_mean;     //  building the vector with mean values
    A_mean[2]= A2_mean;
    A_mean[3]= A3_mean;
    A_mean[4]= A4_mean;

    //  calculate the mean of the total load Z
    Z_mean = A1_mean +  A2_mean + A3_mean + A4_mean ;

    //  caculate Z_var as a function of Z_alpha
    Z_var = square( Z_alpha ) * tgamma( 3 / Z_beta ) / tgamma( 1 / Z_beta ) ;

	for (i in 1:data_n){
        //  caculate the rhos as a function of the total axe load Z data
        rho[i,1,1] = 1;
        rho[i,1,2] = func_linear( Z[i], Z_0, Z_1, rho_c12,  rho_d12 ) ;
        rho[i,1,3] = func_linear( Z[i], Z_0, Z_1, rho_c13,  rho_d13 ) ;
        rho[i,1,4] = func_linear( Z[i], Z_0, Z_1, rho_c14,  rho_d14 ) ;
        rho[i,2,1] = rho[i,1,2];
        rho[i,2,2] = 1;
        rho[i,2,3] = func_linear( Z[i], Z_0, Z_1, rho_c23,  rho_d23 ) ;
        rho[i,2,4] = func_linear( Z[i], Z_0, Z_1, rho_c24,  rho_d24 ) ;
        rho[i,3,1] = rho[i,1,3];
        rho[i,3,2] = rho[i,2,3];
        rho[i,3,3] = 1;
        rho[i,3,4] = func_linear( Z[i], Z_0, Z_1, rho_c34,  rho_d34 ) ;
        rho[i,4,1] = rho[i,1,4];
        rho[i,4,2] = rho[i,2,4];
        rho[i,4,3] = rho[i,3,4];
        rho[i,4,4] = 1;

        //  calculate STD as a function of Z
        A2_std[i] = func_linear( Z[i], Z_0, Z_1, A2_std_c, A2_std_d) ;
        A3_std[i] = func_linear( Z[i], Z_0, Z_1, A2_std_c, A2_std_d) ;  // use A2 values
        A4_std[i] = func_linear( Z[i], Z_0, Z_1, A1_std_c, A1_std_d) ;  // use A1 values

        //  we calculate A1_std as a function of Z_var and the covariance matrix
        p = 2*( A2_std[i] *rho[i,1,2] + A3_std[i]*rho[i,1,3] + A4_std[i]*rho[i,1,4] );
        q = 0;
        q += -Z_var;
        q += A2_std[i]*A2_std[i]*rho[i,2,2] + A2_std[i]*A3_std[i]*rho[i,2,3] + A2_std[i]*A4_std[i]*rho[i,2,4];
        q += A3_std[i]*A2_std[i]*rho[i,3,2] + A3_std[i]*A3_std[i]*rho[i,3,3] + A3_std[i]*A4_std[i]*rho[i,3,4];
        q += A4_std[i]*A2_std[i]*rho[i,4,2] + A4_std[i]*A3_std[i]*rho[i,4,3] + A4_std[i]*A4_std[i]*rho[i,4,4];

        A1_std[i] = -p/2 + sqrt( square( p/2) - q );

        COVMAT[i,1,1] = rho[i,1,1] * A1_std[i] * A1_std[i] + 1e-10;
        COVMAT[i,1,2] = rho[i,1,2] * A1_std[i] * A2_std[i];
        COVMAT[i,1,3] = rho[i,1,3] * A1_std[i] * A3_std[i];
        COVMAT[i,1,4] = rho[i,1,4] * A1_std[i] * A4_std[i];
        COVMAT[i,2,1] = rho[i,2,1] * A2_std[i] * A1_std[i];
        COVMAT[i,2,2] = rho[i,2,2] * A2_std[i] * A2_std[i] + 1e-10;
        COVMAT[i,2,3] = rho[i,2,3] * A2_std[i] * A3_std[i];
        COVMAT[i,2,4] = rho[i,2,4] * A2_std[i] * A4_std[i];
        COVMAT[i,3,1] = rho[i,3,1] * A3_std[i] * A1_std[i];
        COVMAT[i,3,2] = rho[i,3,2] * A3_std[i] * A2_std[i];
        COVMAT[i,3,3] = rho[i,3,3] * A3_std[i] * A3_std[i] + 1e-10;
        COVMAT[i,3,4] = rho[i,3,4] * A3_std[i] * A4_std[i];
        COVMAT[i,4,1] = rho[i,4,1] * A4_std[i] * A1_std[i];
        COVMAT[i,4,2] = rho[i,4,2] * A4_std[i] * A2_std[i];
        COVMAT[i,4,3] = rho[i,4,3] * A4_std[i] * A3_std[i];
        COVMAT[i,4,4] = rho[i,4,4] * A4_std[i] * A4_std[i] + 1e-10;

        L_COVMAT[i,:,:] = cholesky_decompose( COVMAT[i,:,:] );

    }   // end loop over data
}

model {

    //  rho parameter priors
    target += normal_lpdf( rho_c12 | prior_rho_c12, 1 );
    target += normal_lpdf( rho_c13 | prior_rho_c13, 1 );
    target += normal_lpdf( rho_c14 | prior_rho_c14, 1 );
    target += normal_lpdf( rho_c23 | prior_rho_c23, 1 );
    target += normal_lpdf( rho_c24 | prior_rho_c24, 1 );
    target += normal_lpdf( rho_c34 | prior_rho_c34, 1 );

    target += normal_lpdf( rho_d12 | prior_rho_d12, 1 );
    target += normal_lpdf( rho_d13 | prior_rho_d13, 1 );
    target += normal_lpdf( rho_d14 | prior_rho_d14, 1 );
    target += normal_lpdf( rho_d23 | prior_rho_d23, 1 );
    target += normal_lpdf( rho_d24 | prior_rho_d24, 1 );
    target += normal_lpdf( rho_d34 | prior_rho_d34, 1 );

    target += normal_lpdf( A1_std_c | prior_A_std, prior_A_std * prior_A_cov );
    target += normal_lpdf( A2_std_c | prior_A_std, prior_A_std * prior_A_cov );
    target += normal_lpdf( A1_std_d | prior_A_std, prior_A_std * prior_A_cov );
    target += normal_lpdf( A2_std_d | prior_A_std, prior_A_std * prior_A_cov );

    //  load std priors
    for (i in 1: data_n) {
        target += normal_lpdf( A1_std[i] | prior_A_std, prior_A_std * prior_A_cov );
        target += normal_lpdf( A2_std[i] | prior_A_std, prior_A_std * prior_A_cov );
        target += normal_lpdf( A3_std[i] | prior_A_std, prior_A_std * prior_A_cov );
        target += normal_lpdf( A4_std[i] | prior_A_std, prior_A_std * prior_A_cov );
    }

    //  load mean priors
    target += normal_lpdf( A1_mean | prior_A_mean, prior_A_mean * prior_A_cov );
    target += normal_lpdf( A2_mean | prior_A_mean, prior_A_mean * prior_A_cov );
    target += normal_lpdf( A3_mean | prior_A_mean, prior_A_mean * prior_A_cov );
    target += normal_lpdf( A4_mean | prior_A_mean, prior_A_mean * prior_A_cov );

    //  Total load priors
    target += normal_lpdf( Z_mean | prior_A_mean * axes_n, prior_A_cov * prior_A_mean * axes_n );
    target += normal_lpdf( Z_var  | axes_n * square( prior_A_std ), prior_A_cov * axes_n * square( prior_A_std ) );
    target += normal_lpdf( Z_alpha | sqrt( axes_n ) * prior_A_std * tgamma(1/Z_beta) / tgamma(3/Z_beta),  prior_A_cov*sqrt( axes_n ) * prior_A_std * tgamma(1/Z_beta) / tgamma(3/Z_beta) );

    //  Likelihood Total axe load
    target += gen_normal_lpdf( Z | Z_mean, Z_alpha, Z_beta );

    for (i in 1: data_n) {
        A[i,:] ~ multi_normal_cholesky_lpdf(A_mean, L_COVMAT[i,:,:] );
    }
}

I replied to your other post on the parameterization. Usually this kind of thing is a problem if a model’s not identified in the likelihood, but only through the priors.

1 Like