Problems with using Dirichlet Prior in Stan

Modeling
,
1

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);
}"
2

Hello

Could you give a bit of insight into which values of alpha you are using? How big is p in your experiments?