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,:,:] );
}
}