Hi, I have a model which basically generates multi-variate gaussians where the mean vector and the covariance matrix are dynamically computed using some custom functions as shown bellow:

```
functions {
vector computeWeights(int nterm, int nreg, real alpha, int[] nbranch, real[] reg_profile, real[] epochs, vector theta) {
int index = 1;
real y[nterm];
matrix[nterm, nreg] W;
for (term in 1:nterm) {
int nb = nbranch[term];
for (br in 1:nb) {
real t = epochs[index] - epochs[index+br-1];
y[br] = exp(-alpha*t);
}
for (br in 1:(nb-1)) {
y[br] = y[br] - y[br+1];
}
for (reg in 1:nreg) {
real dotp = 0.0;
for (br in 1:nb) {
dotp += y[br]*(reg_profile[index+br-1] == reg ? 1 : 0);
}
W[term,reg] = dotp;
}
index += nb;
}
return (W*theta);
}
matrix computeCovars(int nterm, real alpha, matrix bt, real sigmasq) {
matrix[nterm, nterm] V;
for (termi in 1:nterm) {
for (termj in 1:nterm) {
real ii = bt[termi, termi];
real jj = bt[termj, termj];
real ij = bt[termi, termj];
real part1 = exp(alpha*(-ii - jj + 2*ij));
real part2 = (1-exp(-2*alpha*ij));
V[termi, termj] = part1*part2/(2*alpha);
}
}
return (V*sigmasq);
}
}
// The input data is a vector 'y' of length 'nterm'.
data {
int<lower=0> nterm;
int<lower=0> nreg;
int<lower=0> nbranch[nterm];
int<lower=0> nentries;
matrix[nterm,nterm] branch_times;
real epochs[nentries];
real reg_profiles[nentries];
vector[nterm] y;
}
parameters {
real<lower=0, upper=1e100> alpha;
real<lower=0, upper=1e100> sigmasq;
vector[nreg] theta;
}
transformed parameters {
vector[nterm] mu = computeWeights(nterm, nreg, alpha, nbranch, reg_profiles, epochs, theta);
matrix[nterm, nterm] covar = computeCovars(nterm, alpha, branch_times, sigmasq);
}
model {
target += multi_normal_lpdf(y | mu, covar);
}
```

I just wanted to check if the mean of the samples generated using stan function is similar to the MLE computed using optimizing function. It seems there is a huge difference. Interestingly the values from optimizing somewhat matches with my own implementation using R optim function.

```
fit = stan(file = 'EvoMCMC.stan', data = model_data, init = initf, verbose = TRUE, chains=nchain, iter=10000)
model_obj = stan_model(file='EvoMCMC.stan')
mle = optimizing(model_obj, data=model_data, init=initf)
mymle = myMLE(data, alpha=1, sigmasq=1, theta=rep(prior_theta_mean,1))
cat("\n**** my MLE ****\n")
print(setNames(mymle$par, c("alpha", "sigmasq", "theta")))
cat("\n**** rstan MLE ****\n")
print(mle$par[c('alpha', 'sigmasq', 'theta[1]')])
cat("\n**** rstan MCMC based solution ****\n")
print(setNames(c(mean(extract(fit)$alpha), mean(extract(fit)$sigmasq), mean(extract(fit)$theta)),
c("alpha", "sigmasq", "theta")))
```

Gives the following result:

```
**** my MLE ****
alpha sigmasq theta
0.5520956 214.2985145 10.0782889
**** rstan MLE ****
alpha sigmasq theta[1]
0.554084 214.682116 10.098615
**** rstan MCMC based solution ****
alpha sigmasq theta
5.593688e+97 6.614111e+99 7.592474e+00
```

Could somebody provide me some suggestions on what is going wrong here?

Thank you.