# Problems with using Dirichlet Prior in Stan

So I’m trying to fit a model through Stan which was initially working fine, however after adding the phi parameter and giving it a Dirichlet prior I’m getting a large amount of divergent samples and also some transitions exceeding maximum treedepth. It seems that the Dirichlet prior is causing a lot of issues. Is there a possible workaround for this, or do I just have to use a different prior distribution? The model is included below.

``````functions{
real pi1_lpdf(real theta_1,real lambda, data vector lambdas, data int n, data int p) {

real ht = 1/(1 + exp(-theta_1)); // h(theta_1)
real dht = exp(-theta_1)/(1 + exp(-theta_1))^2;

if (p < n) {
real d = sqrt(sum(ht*lambdas + 1 - ht) - p - sum(log(ht * lambdas + 1 - ht)) +
(n-p)*(-ht + theta_1 + log(1 + exp(-theta_1) )) ); // d(R2)

real ddR = 1/(2*d)*(sum(lambdas) - n - sum((lambdas - 1)./(ht*lambdas + 1 - ht)) + (n-p)/(1-ht));
real ret = lambda * exp(-lambda*d)*ddR*dht;
return log(ret);
}
else {
real d = sqrt(sum(ht*lambdas + 1 - ht) - n - sum(log(ht*lambdas + 1 -ht)));
real ddR = 1/(2*d)*(sum(lambdas -1 ) - sum((lambdas -1)./(ht*lambdas + 1 - ht)));
real ret = lambda * exp(-lambda*d)*ddR*dht;
return log(ret);
}

}

real pi2_lpdf(real theta_2) {

real gt = sqrt(exp(theta_2));

real dsig = 0.5*gt;

real ret = log(2)*exp(-log(2)*gt)*dsig;

return log(ret);

}
}

data {

int<lower = 1> n; // number of data points
int<lower = 1> p; // number of predictors
matrix[n, p] X;   // predictor matrix
vector[n] y;      // outcome vector
real<lower=0> lambda;  // rate parameter for pi(R^2)
vector[min(n,p)] lambdas; // eigenvalues of X %*% t(X)
vector<lower = 0> [p] alpha; // parameter vector for the dirichlet prior

}

parameters {

vector[p] beta;
real theta_1;
real<lower=0> theta_2;
simplex[p] phi;
}

transformed parameters {
real<lower = 0,upper = 1> R2 = inv_logit(theta_1);
real<lower=0> sigma2 = exp(theta_2);
vector <lower=0>[p] sigmab = R2 * sigma2 * phi;
real<lower=0> sigma_n = (1-R2)*sigma2;
}

model {

target += pi1_lpdf(theta_1 | lambda, lambdas,n,p);
target += pi2_lpdf(theta_2);
target += dirichlet_lpdf(phi | alpha); // this is what I believe is causing the issues
beta ~ normal(0,sigmab);
y ~ normal(X * beta, sigma_n);
}"

``````