My priors are not known distribution, how should I represent them in the code?
The likelihood is y \sim N(x\beta, \{tr(v_{\lambda}^{-1})W_c\}^{-1}), and \rho=\sigma^2/\sigma_e^2.
the priors are \sigma_e^2 \sim \frac{1}{\sigma_e^{1.5}}, \rho \sim \rho^{-0.5}g(\rho), and g(\rho)=\sum_{1}^{m}(1+n_{i}\rho)^{-1}
Any insight is appreciated.
data{
int<lower=0> m; //m=27
int<lower=0> n; //n=50
vector[m] ni; //vector n=(n1,n2,..,n_m),size for each region
matrix[n,2] x; //log(x)
vector[n] y; //log(y)
matrix[m,n] N;
matrix[n,n] Q;
matrix[n,n] I; //identity matrix with rank n
}
parameters{
vector[2] beta;
real<lower=0> sigma;
real<lower=0> sigmae;
real<lower=0,upper=1> alpha;
}
transformed parameters{
real<lower=0> rho = sigma / sigmae;
vector[m] B2;
vector[m] nirho;
matrix[n,n] vlambda = sigmae * I + sigma * Q;
matrix[n,n] wj;
matrix[n,n] wm = inverse(vlambda) / sum(diagonal(inverse(vlambda)));
matrix[n,n] wc;
for(i in 1:m){
B2[i] = (sigmae/(sigmae+ni[i]*sigma))*(sigmae/(sigmae+ni[i]*sigma));
nirho[i] = 1/(1+ni[i]*rho);
}
wj = (N'*diag_matrix(B2)*N) / sum(diagonal(N'*diag_matrix(B2)*N));
wc = alpha*wj + (1-alpha)*wm;
}
model{
target += -1.5*log(sigmae)-0.5*log(rho)+log(sum(nirho)); //sigmae ~ sigmae^-1.5; rho ~ rho^-0.5*sum(1+ni*rho)^-1
y ~ multi_normal(x*beta,inverse_spd(trace(vlambdainv)*wc));
}