Why my estimations not converged?

Modeling
,
1

Hi, I am new to stan. I am trying to use rstan to estimate parameters in a multivariate normal model like this
y_{ij}^{(1)}=X*\beta^{(1)}+Z*b_i^{(1)}+e_{ij}^{(1)};
y_{ij}^{(2)}=X*\beta^{(2)}+Z*b_i^{(2)}+e_{ij}^{(2)};
…
with \ell=10 responses, and b\sim Normal(0,\Psi), e_i^{(\ell)}\sim Normal(0,\sigma^{(\ell)})
I have to estimate \beta^{(\ell)}, \Psi, \sigma^{(\ell)} , However, I have problems in converging by using the following stan code.

functions{
  real likelihood_lpdf(real[,] y, real[,,] Z, real[,,] X,int N,real[] quadf,real w){

 real lprob; 
 lprob = 0;
 for(i in 1:N){
       lprob=lprob-(0.5)*log(2*pi())+(0.5)*log(w)-0.5*quadf[i];
}
 return lprob;
 }
}

data{
 int N; //number of observations
 int n; //number of branches
 int resp;//number of responses=l
 int p; //dim of X
 int q; //dim of Z
 matrix[(resp*n),(resp*n)] I1;
 matrix[(resp*q),(resp*q)] I2;
  vector[(resp*q)] bmu;
 real y[(n*resp),N]; //outcomes array
 real X[N,(n*resp),(p*resp)]; //branch-varying variables array
 real Z[N,(n*resp),(q*resp)]; //option-varying variables array
}

parameters{
 vector[(resp*p)] nbeta; //
 matrix[(resp*q),N] B;
 real<lower=0> sigma[resp];
 cov_matrix[(resp*q)] sigmab;
}

transformed parameters{
  real mu[N,(n*resp)];
  real res[N,(n*resp)];
  matrix[(n*resp),(n*resp)] sigmae;
  real w;
  real quadf[N];
  
  sigmae=I1;
  for(j in 0:(resp-1))
  {
    for(k in 1:n)
    {
          sigmae[(j*n+k),(j*n+k)]=sigma[(j+1)];
    }
  }
  w=determinant(inverse(sigmae));
  for(i in 1:N)
  {
    for(l in 1:(n*resp))
    {
      mu[i,l]=dot_product(to_row_vector(X[i,l,]),nbeta)+dot_product(to_row_vector(Z[i,l,]),B[,i]);
    res[i,l]=y[l,i]-mu[i,l];
    }
            quadf[i]=quad_form(inverse(sigmae),to_vector(res[i,]));
  }
}
model{
//specify priors
 for(i in 1:(p*resp)) 
 {
   nbeta[i] ~ normal(1,10);
 }
 for(i in 1:resp)
 {
    sigma[i]~normal(2,10);
 }
 sigmab~inv_wishart((resp*q), I2);
 for(i in 1:N){
    B[,i]~multi_normal(bmu,sigmab);
 }
      y~likelihood(Z,X,N,quadf,w);
}

}

I1<-diag((respn))
I2<-diag((respq))
bmu<-rep(0,(resp*q))

2

I am sorry but the model is quite long and complex and I am not able to understand it. What is the use of likelihood_lpdf? Why are you calculating determinant of inverse of something?

It is generally useful to start with a simple model (e.g. only fitting a multivariete normal) and fit it to simulated data. Then add complexity one step at a time. This way we would know where the problem could be. In this full form, it is very hard to debug the model.

Best of luck!

3

I have fixed it by replacing lpdf with log. It could go reasonable estimations now. However, I don’t know what is the difference between pdf and log.

Thanks