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

Modeling
1

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

2

Can you share the code for this?

3

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));
    }
}
4

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.

5

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

6

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

7

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());
}
8

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.
9

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

10

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.

11

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.