Log probability evaluates to negative infinity

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