# Debugging Log probability evaluates to log(0), i.e. negative infinity

I need help in interpreting the diagnostics I have printed out from a troublesome run of my model.

As my model includes a large set of user-defined functions, I want to try to explain my problem without including all of the functions block, but if any of you very helpful people need to see it just ask.
On some, but not all, runs of my model I get the (by me) dreaded ‘Log probability evaluates to log(0)’ error.
In my efforts to debug I have inserted some print statements into the model block as shown. You will see that I am using user-defined probability distributions in three of which the parameters are supplied as data, and a fourth where the parameters are found by transforming sampled parameters.

One of the quantities I print is denoted `tester`. The log probability is undefined if `tester < 0`. When I run the model I supply initial values for the three parameters `qtilde1, qtilde2, qtilde3` that are guaranteed to give a positive value for `tester`.

``````data{
// parameters of distribution 1
int<lower = 0> N1;
vector[N1] knots1;
vector[N1] weights1;
// parameters of distribution 1
int<lower = 0> N2;
vector[N2] knots2;
vector[N2] weights2;
// parameters of distribution 1
int<lower = 0> N3;
vector[N3] knots3;
vector[N3] weights3;

int<lower=0> N; // number of observations
vector[N] y; // observations

real u; // high threshold

}
parameters {
real<lower = 0> qtilde1;
real<lower = 0>  qtilde2;
real<lower = 0>  qtilde3;
}
transformed parameters{
real mu;
real<lower = 0> nu = qtilde3 / qtilde2;
real xi = log10(nu);
real<lower = 0> sigma;
sigma  = qtilde2 * xi / (nu * (nu - 1));
mu = qtilde1 - qtilde2 / nu;
}
model {
// priors
qtilde1 ~ weighted_hat_lpdf(knots1, weights1, N1);
print("qt1=",qtilde1);
qtilde2 ~ weighted_hat_lpdf(knots2, weights2, N2);
print("qt2=",qtilde2);
qtilde3 ~ weighted_hat_lpdf(knots3, weights3, N3);
print("qt3=",qtilde3);
print("tester= ", min(1 + xi * (y - mu)/ sigma));
print("lp= ", gpoispoint_lpdf(y | mu, xi, sigma,  u))
// log probability
target += gpoispoint_lpdf(y | mu, xi, sigma,  u) ;
print("target = ", target());
}
``````

Here is the start of the printed output (the message is repeated 100 times):

``````Chain 1: qt1=0.878467
qt2=0.703037
qt3=0.497315
tester= 0.622135
lp= -56.6921
target = nan

Chain 1: Rejecting initial value:
Chain 1:   Log probability evaluates to log(0), i.e. negative infinity.
Chain 1:   Stan can't start sampling from this initial value.
Chain 1: qt1=0.878467
qt2=0.703037
qt3=0.497315
tester= 0.622135
lp= -56.6921
target = nan
``````

The values of `qt1, qt2, qt3` are those I supplied. The positive value of `tester` shows that they are valid and the printout shows that they do indeed lead to a permissible value for the log probability (`lp`).

So why is `target = nan`?

I am running this in rstan 2.19.3 on Rstudio 1.2.5033 on a Macbook Pro running Catalina 10.15.5

Can you share the code for this?

Here is the whole function block:

``````functions {

real gpoispoint_lpdf(vector y, real mu, real k, real sigma, real u) {
// generalised Pareto log pdf
int N = rows(y);
real inv_k = inv(k);
if (k<0 && max(y-mu)/sigma > -inv_k)
reject("k<0 and max(y-mu)/sigma > -1/k; found k, sigma =", k, sigma)
if (sigma<=0)
reject("sigma<=0; found sigma =", sigma)
if (fabs(k) > 1e-15)
return -(1+inv_k)*sum(log1p((y-mu) * (k/sigma))) -N*log(sigma) -
54 * (1 + k * (u - mu) / sigma) ^ -inv_k;
else
return -sum(y-mu)/sigma -N*log(sigma) -
(54 * exp(u - mu) / sigma); // limit k->0
}

real gpoispoint_cdf(vector y, real mu, real k, real sigma) {
//   cdf
real inv_k = inv(k);
if (k<0 && max(y-mu)/sigma > -inv_k)
reject("k<0 and max(y-mu)/sigma > -1/k; found k, sigma =", k, sigma)
if (sigma<=0)
reject("sigma<=0; found sigma =", sigma)
if (fabs(k) > 1e-15)
return exp(sum(log1m_exp((-inv_k)*(log1p((y-mu) * (k/sigma))))));
else
return exp(sum(log1m_exp(-(y-mu)/sigma))); // limit k->0
}
real gpoispoint_lcdf(vector y, real mu, real k, real sigma) {
//  log cdf
real inv_k = inv(k);
if (k<0 && max(y-mu)/sigma > -inv_k)
reject("k<0 and max(y-mu)/sigma > -1/k; found k, sigma =", k, sigma)
if (sigma<=0)
reject("sigma<=0; found sigma =", sigma)
if (fabs(k) > 1e-15)
return sum(log1m_exp((-inv_k)*(log1p((y-mu) * (k/sigma)))));
else
return sum(log1m_exp(-(y-mu)/sigma)); // limit k->0
}

real gpoispoint_lccdf(vector y, real mu, real k, real sigma) {
//   log ccdf
real inv_k = inv(k);
if (k<0 && max(y-mu)/sigma > -inv_k)
reject("k<0 and max(y-mu)/sigma > -1/k; found k, sigma =", k, sigma)
if (sigma<=0)
reject("sigma<=0; found sigma =", sigma)
if (fabs(k) > 1e-15)
return (-inv_k)*sum(log1p((y-mu) * (k/sigma)));
else
return -sum(y-mu)/sigma; // limit k->0
}
real gpoispoint_rng(real mu, real k, real sigma) {
//  rng
if (sigma<=0)
reject("sigma<=0; found sigma =", sigma)
if (fabs(k) > 1e-15)
return mu + (uniform_rng(0,1)^-k -1) * sigma / k;
else
return mu - sigma*log(uniform_rng(0,1)); // limit k->0
}

real felt_p1(real x){
return (1 - cos(pi() * (x + 1))) / 2;
}

real felt_p2(real x) {
return (1 + cos(pi() * x)) / 2;
}

real felt_c1(real x) {
return (x + sin(pi() * x) / pi()) / 2;
}

real felt_c2(real x) {
return (x + sin(pi() * x) / pi()) / 2;
}

real cf(real x, real lo, real mid, real hi){
if (x < lo)
return 0;
else if (x <= mid)
return (mid - lo) * (felt_c1((x - mid)/(mid - lo)) - felt_c1(-1)) ;
else if (x  < hi)
return (hi - mid) * (felt_c2((x - mid)/(hi - mid))) + (mid - lo) / 2;
else
return (hi - lo) / 2;
}

real pf(real x, real lo, real mid, real hi){
real x1 = (x - mid)/(mid - lo);
real x2 = (x - mid) / (hi - mid);
if ((x < lo) || (x > hi))
return 0;
else if (x <= mid)
return felt_p1(x1);
else
return felt_p2(x2);
}

real wted_cdf(real x, int j, vector knots, vector weights, int N){
real lo;
real mid = knots[j];
real hi;
if(j == 1)
lo = 0;
else
lo = knots[j - 1];
if(j == N)
hi = 1;
else
hi = knots[j + 1];
return weights[j] * cf(x , lo, mid, hi);
}

real wted_pdf(real x, int j, vector knots, vector weights, int N){
real lo;
real mid = knots[j];
real hi;
if (j == 1)
lo = 0;
else
lo = knots[j - 1];
if (j == N)
hi = 1;
else
hi = knots[j + 1];
return weights[j] * pf(x , lo, mid, hi);
}

real weighted_hat_lpdf(real y, vector knots, vector weights, int N){
vector[N] pfn;
real z = y / (y + 1);
for (n in 1:N){
pfn[n] = wted_pdf(z, n, knots, weights, N);
}
return log(sum(pfn));
}
}

``````

Did you print `target` before adding `gpoispoint`? Are you sure all values in `knots` are distinct? If two consecutive `knots` are equal that `pf()` is going to divide by zero.

Thank you for the question. The knots are distinct by construction.

I don’t quite understand your question about printing. Can you put it another way?

You can print the `target` at any point. Since the `gpoispoint_lpdf` appears finite I wonder if `target` was NaN way before the last line was added.

``````model {
print("target = ", target());
// priors
qtilde1 ~ weighted_hat_lpdf(knots1, weights1, N1);
print("qt1=",qtilde1);
print("target = ", target());
qtilde2 ~ weighted_hat_lpdf(knots2, weights2, N2);
print("qt2=",qtilde2);
print("target = ", target());
qtilde3 ~ weighted_hat_lpdf(knots3, weights3, N3);
print("qt3=",qtilde3);
print("target = ", target());
print("tester= ", min(1 + xi * (y - mu)/ sigma));
print("lp= ", gpoispoint_lpdf(y | mu, xi, sigma,  u))
// log probability
target += gpoispoint_lpdf(y | mu, xi, sigma,  u) ;
print("target = ", target());
}
``````

I have just tried out your idea. Unfortunately if I insert the `print(target)` command before `target += ...` the parser rejects it :

``````SYNTAX ERROR, MESSAGE(S) FROM PARSER:
Variable "target" does not exist.
``````

That’s a weird error message…
Anyway, it’s `print(target())`, with extra `()` after `target`.

I typed the command correctly in the code itself.

I can now reproduce the `target = nan` message even if I replace `gpoispoint_lpdf(y | mu, xi, sigma, u)` with `normal_lpdf(y| mu, sigma)` so the problem does not lie in that user-defined distribution function.

In that case the most likely culprit is one of `qtildeN ~ weighted_hat`. Replace those one-by-one and see if `target` is still nan.