Hi Stan folks,
This might be less of a Stan question (although maybe there’s a Stan solution) and more of a model-fitting question, but here goes.
In part of my model, I have
\pi_j = 1 - \exp{\left(-\dfrac{r_j^2}{\tau^2}\right)}
Where I want \tau \sim MVN(\mu, \Sigma).
Of course, there is the problem where \tau cannot be 0; in fact it is known that \tau is strictly positive. However, there does not exist a truncated multinormal distribution in Stan. When I go to run the model, Stan gives me the "Initialize between (-2,2) failed after 100 attempts.
I’m wondering how I might go about fixing this. In the univariate version, I set \tau \sim lognormal(1,1), which resolves the issue of \tau being sample as 0. Is there a way to force Stan to try initializing between, say, (1,3)?
The following is my model, for reference:
data {
int<lower = 1> n_samples; // total number of sampling events i
int<lower = 1> n_species; // total number of specise
int<lower = 2> max_intervals; // maximum number of intervals being considered
int species[n_samples]; // species being considered for each sample
int abund_per_band[n_samples, max_intervals];// abundance in distance band k for sample i
vector[n_samples] abund_per_sample; // total abundnace for sample i
int bands_per_sample[n_samples]; // number of distance bands for sample i
matrix[n_samples, max_intervals] max_dist; // max distance for distance band k
corr_matrix[n_species] phylo_corr; // correlation matrix of phylogeny
}
parameters {
row_vector<lower = 0>[n_species] mu;
vector<lower = 0>[n_species] sigma;
vector[n_species] log_tau;
}
model {
matrix[n_samples, max_intervals] Pi; // probabilities
sigma ~ cauchy(0, 2.5);
mu ~ lognormal(1, 1);
log_tau ~ multi_normal(mu, quad_form_diag(phylo_corr, sigma));
Pi = rep_matrix(0, n_samples, max_intervals);
for (i in 1:n_samples)
{
for (k in 2:bands_per_sample[i])
{
Pi[i,k] = ((1 - exp(-(max_dist[i,k]^2 / exp(log_tau[species[i]])^2))) -
(1 - exp(-(max_dist[i,k - 1]^2 / exp(log_tau[species[i]])^2)))) /
(1 - exp(-(max_dist[i,bands_per_sample[i]]^2 / exp(log_tau[species[i]])^2)));
}
Pi[i,1] = 1 - sum(Pi[i,]);
abund_per_band[i,] ~ multinomial(to_vector(Pi[i,]));
}
}
Thanks for any assistance!