Hi,
I’m trying to fit a model but I can’t figure out how to implement a constraint on a variable which is some function of some parameters. Right at the top of the function df_lpdf
I have this E
variable, and I want it to be positive. I have tried defining it inside the transformed parameters block instead but the constraints there are only for checks, and not actual constraints. Do I need to reparameterize this, and if so, how can I do so? Right now the model sometimes runs and sometimes doesn’t because of bad starting values. If it runs then I’m getting n_eff
< 100 from 1000 draws.
Also, is the way I currently have it set up the proper way to pass multiple inputs to a custom distribution function (the real[] y
part)? It seems a little clunky but I couldn’t find any examples in the documentation.
Sorry if this is an obvious problem, this is my first time using Stan so I am unfamiliar with it.
Thanks! :)
functions {
real df_lpdf(real[] y, real phi0, real g, real b, real a) {
// E should be > 0
real E = -square(y[1]) / 2. + phi0 / pow(y[3], g);
real L = y[3] * y[2];
real num_t1 = -2 * b * log(L);
real num_t2 = (b * (g - 2.) / g + a / g - 3. / 2.) * log(E);
real num_t3 = lgamma(a / g - 2. * b / g + 1.);
real denom_t1 = 0.5 * log(8. * pow(pi(), 3.) * pow(2., -2. * b));
real denom_t2 = (-2. * b / g + a / g) * log(phi0);
real denom_t3 = lgamma(b * (g - 2.) / g + a / g - 1. / 2.);
real numerator = num_t1 + num_t2 + num_t3;
real denominator = denom_t1 + denom_t2 + denom_t3;
// this ~should~ never happen but it does.
if (E < 0.) {
return negative_infinity();
}
return numerator - denominator;
}
vector get_angles(real x, real y, real z) {
real projR = sqrt(square(x) + square(y));
real R = sqrt(square(x) + square(y) + square(z));
real sintheta = y / projR;
real costheta = x / projR;
real sinphi = projR / R;
real cosphi = z / R;
return to_vector([sintheta, costheta, sinphi, cosphi]);
}
vector transform_vels(real ra, real dec, real pmra, real pmdec, real vlos, real plx, real sintheta, real costheta, real sinphi, real cosphi) {
matrix[3,3] T = [
[-0.06699, -0.87276, -0.48354],
[ 0.49273, -0.45035, 0.74458],
[-0.86760, -0.18837, 0.46020]
];
real deg2rad = pi() / 180.;
real sinra = sin(ra * deg2rad);
real cosra = cos(ra * deg2rad);
real sindec = sin(dec * deg2rad);
real cosdec = cos(dec * deg2rad);
matrix[3,3] A = [
[cosra * cosdec, -sinra, -cosra * sindec],
[sinra * cosdec, cosra , -sinra * sindec],
[sindec, 0, cosdec]
];
matrix[3, 3] B = T * A;
real k = 4.74057;
vector[3] solarmotion = to_vector([11.1, 232.24, 7.25]);
vector[3] dat = to_vector([vlos, k * pmra / plx, k * pmdec / plx]);
vector[3] uvw = B * dat + solarmotion;
matrix[3, 3] ptz_mat = [
[costheta, sintheta, 0],
[-sintheta, costheta, 0],
[0, 0, 1]
];
vector[3] ptz = ptz_mat * uvw;
matrix[3, 3] rtp_mat = [
[sinphi, 0, cosphi],
[0, 1, 0],
[cosphi, 0, -sinphi]
];
vector[3] rtp = rtp_mat * ptz;
return rtp / 100;
}
}
data {
int<lower=1> N;
vector[N] ra;
vector[N] dec;
vector[N] plx;
vector[N] Xgc;
vector[N] Ygc;
vector[N] Zgc;
// measurements
vector[N] pmra_true;
vector[N] pmdec_true;
vector[N] vlos_true;
vector[N] r_true;
}
transformed data {
vector[N] sintheta;
vector[N] costheta;
vector[N] sinphi;
vector[N] cosphi;
real<lower=0> y[N, 3];
for (i in 1:N) {
vector[4] angles = get_angles(Xgc[i], Ygc[i], Zgc[i]);
sintheta[i] = angles[1];
costheta[i] = angles[2];
sinphi[i] = angles[3];
cosphi[i] = angles[4];
}
{
for (i in 1:N) {
vector[3] vels_sph = transform_vels(ra[i], dec[i], pmra_true[i], pmdec_true[i], vlos_true[i], plx[i], sintheta[i], costheta[i], sinphi[i], cosphi[i]);
y[i, 2] = sqrt(square(vels_sph[2]) + square(vels_sph[3]));
y[i, 1] = sqrt(square(vels_sph[1]) + square(y[i, 2]));
y[i, 3] = r_true[i];
}
}
}
parameters {
real<lower=0> p_phi0;
real<lower=0> p_gamma;
real<upper=1> p_beta;
real<lower=(p_beta * (2. - p_gamma) + p_gamma / 2.)> p_alpha;
}
model {
// priors ----------------
p_phi0 ~ uniform(1., 200);
p_gamma ~ normal(0.5, 0.06);
p_beta ~ uniform(-0.5, 1.);
p_alpha ~ gamma(2.99321604, 2.82409927);
// likelihood ---------------------
for (i in 1:N) {
y[i] ~ df(p_phi0, p_gamma, p_beta, p_alpha);
}
}