# Why my estimations not converged?

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){
}
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;

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

}


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

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!

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