I was hoping to be able to print target
, to figure out where it’s getting messed up… but target
is not a variable.
Any thoughts?
functions {
int real_to_int(real x) {
int i = 0;
while(i < x) {
i += 1;
}
return i;
}
real interpolate(vector F, real x, int T) {
int i;
real res;
real flo = floor(x*T); // possible values: {0 ,..., T-1}
real lambda = x*T - flo;
i = real_to_int(flo) + 1;
res = (1 - lambda) * F[i] + lambda * F[i+1];
return(res);
}
}
data {
int N;
int T; // number of line segments (granularity of Gaussian Process)
int<lower=0> Y[N+1]; // Number of Y samples in each bucket
}
transformed data{
int N_Y = sum(Y);
}
parameters {
simplex[N+1] intervals;
positive_ordered[T+1] F; // Cumulative density ratio.
real<lower=0> sigma; //how quickly F can vary
}
model {
real l = 0; // left end of interval
real r; // right end of interval
real rate; // Poisson rate
real F_l = 0; // F(0) = 0
real F_r;
intervals ~ dirichlet(rep_vector(1.0, N+1));
for (k in 1:(N+1)) { // for each bucket
r = l + intervals[k];
F_r = interpolate(F, r, T); //interpolate
rate = F_r - F_l; // the Poisson rate is the integral
Y[k] ~ poisson(N_Y * rate);
l = r;
F_l = F_r;
}
}
1 Like
It can be accessed as target()
.
2 Likes
I’ve commented out all ~
statements, and all references to target
… but it’s still returning negative infinity.
If you can post your model code and some example data I can debug locally
3 Likes
Thanks so much for the offer! I made some changes, and it now looks like the main issue is divergent transitions. I’ll work on it again tomorrow, and update this thread.
My issue now is that all my transitions are divergent.
I’ve tried increasing adapt_delta
to 0.99 and 0.999 but it doesn’t help.
functions {
int real_to_int(real x) {
int i = 0;
while(i < x) {
i += 1;
}
return i;
}
// Should we vectorize this and have it done all at once, i.e. real[] x?
real interpolate(vector F, real x, int T) {
// i and i+1
int i;
real res;
real flo = floor(x*T); // possible values: {0 ,..., T-1}
real lambda = x*T - flo;
i = real_to_int(flo) + 1;
res = (1 - lambda) * F[i] + lambda * F[i+1];
return(res);
}
}
data {
int<lower=0> N; // number of X obs, yielding N+1 buckets
int<lower=0> T; // number of line segments (granularity of Gaussian Process)
int<lower=0> Y[N+1]; // Number of Y samples in bucket k
}
transformed data{
int N_Y = sum(Y);
real empty_r[0];
int empty_i[0];
}
parameters {
simplex[N+1] intervals; //dirichlet
positive_ordered[T+1] F; // Cumulative density ratio.
real<lower=0> sigma; //how quickly the log-density ratio can vary
}
transformed parameters {
real<lower=0> dF[T]; // density ratio
real ddF[T-1]; // 1st deriv of density ratio
for (t in 1:T)
dF[t] = F[t+1] - F[t];
for (t in 1:(T-1))
ddF[t] = dF[t+1] - dF[t];
}
// ToDo: change this to a non-centered implementation?
model {
print("BEFORE: target = ", target());
{
real l = 0; // left end of interval
real r; // right end of interval
real rate; // Poisson rate
real F_l = 0; // F(0) = 0
real F_r;
intervals ~ dirichlet(rep_vector(1.0, N+1));
for (k in 1:(N+1)) { // for each bucket
print("----------");
print("k = ", k);
r = l + intervals[k];
print("l = ", l);
print("r = ", r);
F_r = interpolate(F, r, T); //interpolate
print("F_l = ", F_l);
print("F_r = ", F_r);
rate = F_r - F_l; // the Poisson rate is the integral
print("rate = ", rate);
print("N_Y = ", N_Y);
Y[k] ~ poisson(N_Y * rate);
l = r;
F_l = F_r;
}
/////////////////////
// Smoothness prior
// ToDo: Define the adjacent differences of logs, in 'transformed_parameters'
for (t in 1:(T-1)) { // for each section of the function approximation
dF[t+1] ~ lognormal(dF[t], sigma);
}
sigma ~ exponential(1);
// Set F(0) = 0.1. Doing this here avoid having to use 'min' above.
F[1] ~ normal(0.1, 1e-9);
}
}
1 Like
Sorry, can’t look into this very deeply, but plese check out Divergent transitions - a primer for some tips. If you have trouble applying the advice there or can’t make progress for any other reason, you are welcome to ask for clarifications here (this stuff is hard).
Best of luck with your model!
1 Like