Hello,

When I run this model, all my transitions are divergent. I’m trying to figure out what a “non-centered” model would be here.

This model is a Poisson Process, in which we have multiple Poisson counts with rate \int_l^r f(x) = F(r) – F(l), for many different intervals [l,r]. These intervals are also unknown: we only know that they partition the unit interval and that their lengths can be assumed to follow a Dirichlet(1, …, 1) prior.

I’m using a `positive_ordered`

to model the cumulative rate function F. Would it be better to instead use `real<lower=0> f[]`

to represent the rate directly? Is one of these parametrizations more likely to cause problems with divergent transitions?

```
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;
print("rate = ", rate);
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;
}
/////////////////////
// Smoothing prior
for (t in 1:(T-1)) { // for each section of the function approximation
dF[t+1] ~ lognormal(dF[t], sigma);
}
// Set F(0) = 0.1. Doing this here avoid having to use 'min' above.
F[1] ~ normal(0.1, 1e-9);
}
}
```