# Integrate_1d in likelihood computation produces NaNs in the gradient

Dear Stan Community,

I am trying to use the integrate_1d function to integrate parameters out of my likelihood, but I get the error Exception: gradient_of_f: The gradient of f is nan for parameter 0.

I vastly simplified the model to only one parameter, but the problem persists. The code is as follows:

functions {
real pcm_likelihood(
real x,                    // Function argument
real xc,                   // Complement of function argument
real[] pars,               // parameters
data real[] data_real,     // data (real)
data int[] data_int        // data (integer)
) {
int J = data_int;
int K = data_int;

int mu_index = 1;

vector[J + 1] item_prob;
vector[K] log_lik;

matrix[J, K] xmbetas = x - to_matrix(data_real, J, K);

for (k in 1:K) {
item_prob = log_softmax(
append_row(0, cumulative_sum(xmbetas[, k]))
);
log_lik[k] = item_prob[data_int[2 + k]];
}

print("pcm_likelihood x = ", x);

return exp(sum(log_lik) + normal_lpdf(x | pars[mu_index], 1));
}
}

data {
// Counts
int<lower=1> I; // No. persons
int<lower=1> K; // No. items
int<lower=1> J; // No. parameters per item

// Data
int<lower=1, upper=J+1> y[I, K];  // Data as Persons by Items Matrix

// Integration parameters
real tol;

// Item Parameters
real beta[K, J];
}

transformed data {
int<lower=1> integration_data[I, K + 2];

for (i in 1:I) {
integration_data[i] = append_array({ J, K }, to_array_1d(y[i,]));
}
}

parameters {
// Person Parameters
real mu_theta;
}

model {
// Priors
target += std_normal_lpdf(mu_theta);

print("mu_theta = ", mu_theta, " ----------------------------------");

// Likelihood
for (i in 1:I) {
real marg_lik = integrate_1d(
pcm_likelihood,
negative_infinity(),
positive_infinity(),
{ mu_theta },
to_array_1d(beta),
integration_data[i],
tol
);

print("marg_lik = ", marg_lik, " ----------------------------------");
target += log(marg_lik);
}
}


When I run this model with I=1, I get the following (truncated output):

Chain 1: mu_theta = -0.534418 ----------------------------------
pcm_likelihood x = 1.79769e+308
pcm_likelihood x = -1.79769e+308
pcm_likelihood x = 0
pcm_likelihood x = 3.08829
pcm_likelihood x = -3.08829
[...]
pcm_likelihood x = 87.715
pcm_likelihood x = -87.715
pcm_likelihood x = 124.21
pcm_likelihood x = -124.21
marg_lik = 0.0333858 ----------------------------------

Chain 1: mu_theta    = -0.534418 ----------------------------------
pcm_likelihood x = 1.79769e+308
pcm_likelihood x = -1.79769e+308
pcm_likelihood x = 0
pcm_likelihood x = 3.08829
pcm_likelihood x = -3.08829
[...]
pcm_likelihood x = 87.715
pcm_likelihood x = -87.715
pcm_likelihood x = 124.21
pcm_likelihood x = -124.21
pcm_likelihood x = 1.79769e+308

Chain 1: Exception: gradient_of_f: The gradient of f is nan for parameter 0  (in 'modele66b9d35cc3_pcm_bs_marg_nobeta_discourse' at line 69)


So what happens, essentially, is that the integrator runs through clean (the computed integral is correct btw, so the integrator works and the function is defined correctly), then starts again with exactly the same parameter with the same nodes and then tries to start again (probably for autodiff?) with \approx+\infty and fails. I would be happy if someone could explain what is going on there?

From @nhuurre’s suggestion in this thread I figured it could be log_softmax and friends. So I would like to use

for (k in 1:K) {
log_lik[k] = categorical_logit_lpmf(
data_int[2 + k] |
append_row(0, cumulative_sum(xmbetas[, k]))
);
}


But since I want to integrate over the whole real number line, I get an error saying that in categorical_logit_lpmf, probabilities can’t be infinite (don’t recall the exact wording).

Following @bbbales2 suggestion in the aforementioned thread, I don’t see what could go wrong with
\displaystyle\frac{\partial}{\partial \mu_{\theta}}\int f(x, \mu_{\theta}) dx,
since in my case,
f(x, \mu_{\theta})=g(x)N(x;\mu_{\theta},1)
with g a fraction of sums of exponentials, so my gradient is
\displaystyle\frac{\partial}{\partial \mu_{\theta}}\int f(x, \mu_{\theta}) dx=\int (x-\mu_{\theta}) f(x, \mu_{\theta}) dx,
which should be well behaved (?).

For now, I am using adaptive Gauss-Hermite, which seems to work but is incredibly slow, plus I think the boost integrator has been shown to be optimal, so I’m hoping there is a way to make this work.

Other questions:

• Is there another way to program my integrand pcm_likelihood without the use of categorical_logit_lpmf, or, alternatively, can I make categorical_logit_lpmf work with \infty?
• How can I “debug” this better? How can I see what happens with the autodiff to find out where exactly it goes wrong?
• Is there a way to provide an analytical derivative? Can I program my own integrated_categorical_lpmf and add the proper derivatives to avoid numerical differentiation?

Thank you!

Hmm this is probably some at infinity sorta problem yeah.

It happens that out at the limits of double precision, you’ll sometimes get NaNs in the gradients but at infinities if the integral is gonna exist it better go to zero.

We have a little logic to allow for this:

if (is_nan(gradient)) {
if (fx.val() == 0) {
} else {
"is nan for parameter ", "");
}
}


I’m guessing what is happening is that the function is evaluating somewhere way out and is returning something very near zero but the gradient is going to NaN.

Can you do this as a print:

real retval = exp(sum(log_lik) + normal_lpdf(x | pars[mu_index], 1));
print("pcm_likelihood x = ", x, ", retval = ", retval);
return retval;


This way we know what the input and output are as everything explodes. Hopefully we get all the prints and Stan isn’t eating them somewhere internally.

Code is from here

It’s a popular idea but hasn’t been done yet. Stan’s autodiff implemented in Stan is the dream!

1 Like

My guess is that at extremely large values normal_lpdf(x| ...) becomes negative infinity. exp(...) takes that to zero but cannot propagate the gradient correctly. You could try to cut off the bad derivative like this

real retval = exp(sum(log_lik) + normal_lpdf(x | pars[mu_index], 1));
if (retval == 0.0) {
return 0.0;
} else {
return retval;
}

1 Like

Thank you for the quick response!

Here is the output:

Chain 1: mu_theta = 1.58279 ----------------------------------
pcm_likelihood x = 1.79769e+308, retval = -nan
pcm_likelihood x = -1.79769e+308, retval = 0
pcm_likelihood x = 0, retval = 0.00381785
pcm_likelihood x = 3.08829, retval = 3.97431e-09
pcm_likelihood x = -3.08829, retval = 2.10173e-08
pcm_likelihood x = 148.993, retval = 0
pcm_likelihood x = -148.993, retval = 0
pcm_likelihood x = 3.41229e+06, retval = 0
pcm_likelihood x = -3.41229e+06, retval = 0
pcm_likelihood x = 2.06933e+18, retval = 0
pcm_likelihood x = -2.06933e+18, retval = 0
pcm_likelihood x = 2.087e+50, retval = 0
pcm_likelihood x = -2.087e+50, retval = 0
pcm_likelihood x = 2.01977e+137, retval = 0
pcm_likelihood x = -2.01977e+137, retval = 0
pcm_likelihood x = 0.913049, retval = 0.000781258
pcm_likelihood x = -0.913049, retval = 0.000971339
pcm_likelihood x = 14.1579, retval = 1.00068e-71
pcm_likelihood x = -14.1579, retval = 5.22716e-67
pcm_likelihood x = 6704.22, retval = 0
pcm_likelihood x = -6704.22, retval = 0
pcm_likelihood x = 9.64173e+10, retval = 0
pcm_likelihood x = -9.64173e+10, retval = 0
pcm_likelihood x = 2.50895e+30, retval = 0
pcm_likelihood x = -2.50895e+30, retval = 0
pcm_likelihood x = 1.44726e+83, retval = 0
pcm_likelihood x = -1.44726e+83, retval = 0
pcm_likelihood x = 0.407298, retval = 0.00271993
pcm_likelihood x = -0.407298, retval = 0.00290594
pcm_likelihood x = 1.68207, retval = 2.86832e-05
pcm_likelihood x = -1.68207, retval = 5.37362e-05
pcm_likelihood x = 6.1509, retval = 4.64787e-21
pcm_likelihood x = -6.1509, retval = 3.05067e-19
pcm_likelihood x = 40.0396, retval = 0
pcm_likelihood x = -40.0396, retval = 0
pcm_likelihood x = 792.92, retval = 0
pcm_likelihood x = -792.92, retval = 0
pcm_likelihood x = 0.198135, retval = 0.00349845
pcm_likelihood x = -0.198135, retval = 0.00360028
pcm_likelihood x = 0.640156, retval = 0.00170476
pcm_likelihood x = -0.640156, retval = 0.00192914
pcm_likelihood x = 1.24893, retval = 0.000220718
pcm_likelihood x = -1.24893, retval = 0.000320677
pcm_likelihood x = 2.26608, retval = 1.02819e-06
pcm_likelihood x = -2.26608, retval = 2.88847e-06
pcm_likelihood x = 4.29646, retval = 2.56333e-13
pcm_likelihood x = -4.29646, retval = 3.61444e-12
pcm_likelihood x = 9.13029, retval = 1.17296e-36
pcm_likelihood x = -9.13029, retval = 9.2339e-34
pcm_likelihood x = 23.1111, retval = 2.33338e-161
pcm_likelihood x = -23.1111, retval = 2.14002e-153
pcm_likelihood x = 74.2771, retval = 0
pcm_likelihood x = -74.2771, retval = 0
pcm_likelihood x = 326.721, retval = 0
pcm_likelihood x = -326.721, retval = 0
pcm_likelihood x = 0.0983968, retval = 0.00372374
pcm_likelihood x = -0.0983968, retval = 0.00377558
pcm_likelihood x = 0.300606, retval = 0.00315605
pcm_likelihood x = -0.300606, retval = 0.00330366
pcm_likelihood x = 0.519858, retval = 0.00222271
pcm_likelihood x = -0.519858, retval = 0.00243573
pcm_likelihood x = 0.770362, retval = 0.0012104
pcm_likelihood x = -0.770362, retval = 0.00142726
pcm_likelihood x = 1.07131, retval = 0.000447644
pcm_likelihood x = -1.07131, retval = 0.000595836
pcm_likelihood x = 1.45057, retval = 8.99144e-05
pcm_likelihood x = -1.45057, retval = 0.000146111
pcm_likelihood x = 1.95078, retval = 6.6724e-06
pcm_likelihood x = -1.95078, retval = 1.49604e-05
pcm_likelihood x = 2.64003, retval = 9.19968e-08
pcm_likelihood x = -2.64003, retval = 3.42542e-07
pcm_likelihood x = 3.63137, retval = 6.36456e-11
pcm_likelihood x = -3.63137, retval = 5.20831e-10
pcm_likelihood x = 5.11992, retval = 1.44761e-16
pcm_likelihood x = -5.11992, retval = 4.03225e-15
pcm_likelihood x = 7.45666, retval = 2.02808e-27
pcm_likelihood x = -7.45666, retval = 3.9524e-25
pcm_likelihood x = 11.3023, retval = 1.84534e-50
pcm_likelihood x = -11.3023, retval = 8.89663e-47
pcm_likelihood x = 17.9641, retval = 1.41078e-105
pcm_likelihood x = -17.9641, retval = 1.76493e-99
pcm_likelihood x = 30.1781, retval = 1.08137e-256
pcm_likelihood x = -30.1781, retval = 3.6088e-246
pcm_likelihood x = 54.0388, retval = 0
pcm_likelihood x = -54.0388, retval = 0
pcm_likelihood x = 104.108, retval = 0
pcm_likelihood x = -104.108, retval = 0
marg_lik = 0.00509291 ----------------------------------

Chain 1: mu_theta = 1.58279 ----------------------------------
pcm_likelihood x = 1.79769e+308, retval = -nan
pcm_likelihood x = -1.79769e+308, retval = 0
pcm_likelihood x = 0, retval = 0.00381785
pcm_likelihood x = 3.08829, retval = 3.97431e-09
pcm_likelihood x = -3.08829, retval = 2.10173e-08
pcm_likelihood x = 148.993, retval = 0
pcm_likelihood x = -148.993, retval = 0
pcm_likelihood x = 3.41229e+06, retval = 0
pcm_likelihood x = -3.41229e+06, retval = 0
pcm_likelihood x = 2.06933e+18, retval = 0
pcm_likelihood x = -2.06933e+18, retval = 0
pcm_likelihood x = 2.087e+50, retval = 0
pcm_likelihood x = -2.087e+50, retval = 0
pcm_likelihood x = 2.01977e+137, retval = 0
pcm_likelihood x = -2.01977e+137, retval = 0
pcm_likelihood x = 0.913049, retval = 0.000781258
pcm_likelihood x = -0.913049, retval = 0.000971339
pcm_likelihood x = 14.1579, retval = 1.00068e-71
pcm_likelihood x = -14.1579, retval = 5.22716e-67
pcm_likelihood x = 6704.22, retval = 0
pcm_likelihood x = -6704.22, retval = 0
pcm_likelihood x = 9.64173e+10, retval = 0
pcm_likelihood x = -9.64173e+10, retval = 0
pcm_likelihood x = 2.50895e+30, retval = 0
pcm_likelihood x = -2.50895e+30, retval = 0
pcm_likelihood x = 1.44726e+83, retval = 0
pcm_likelihood x = -1.44726e+83, retval = 0
pcm_likelihood x = 0.407298, retval = 0.00271993
pcm_likelihood x = -0.407298, retval = 0.00290594
pcm_likelihood x = 1.68207, retval = 2.86832e-05
pcm_likelihood x = -1.68207, retval = 5.37362e-05
pcm_likelihood x = 6.1509, retval = 4.64787e-21
pcm_likelihood x = -6.1509, retval = 3.05067e-19
pcm_likelihood x = 40.0396, retval = 0
pcm_likelihood x = -40.0396, retval = 0
pcm_likelihood x = 792.92, retval = 0
pcm_likelihood x = -792.92, retval = 0
pcm_likelihood x = 0.198135, retval = 0.00349845
pcm_likelihood x = -0.198135, retval = 0.00360028
pcm_likelihood x = 0.640156, retval = 0.00170476
pcm_likelihood x = -0.640156, retval = 0.00192914
pcm_likelihood x = 1.24893, retval = 0.000220718
pcm_likelihood x = -1.24893, retval = 0.000320677
pcm_likelihood x = 2.26608, retval = 1.02819e-06
pcm_likelihood x = -2.26608, retval = 2.88847e-06
pcm_likelihood x = 4.29646, retval = 2.56333e-13
pcm_likelihood x = -4.29646, retval = 3.61444e-12
pcm_likelihood x = 9.13029, retval = 1.17296e-36
pcm_likelihood x = -9.13029, retval = 9.2339e-34
pcm_likelihood x = 23.1111, retval = 2.33338e-161
pcm_likelihood x = -23.1111, retval = 2.14002e-153
pcm_likelihood x = 74.2771, retval = 0
pcm_likelihood x = -74.2771, retval = 0
pcm_likelihood x = 326.721, retval = 0
pcm_likelihood x = -326.721, retval = 0
pcm_likelihood x = 0.0983968, retval = 0.00372374
pcm_likelihood x = -0.0983968, retval = 0.00377558
pcm_likelihood x = 0.300606, retval = 0.00315605
pcm_likelihood x = -0.300606, retval = 0.00330366
pcm_likelihood x = 0.519858, retval = 0.00222271
pcm_likelihood x = -0.519858, retval = 0.00243573
pcm_likelihood x = 0.770362, retval = 0.0012104
pcm_likelihood x = -0.770362, retval = 0.00142726
pcm_likelihood x = 1.07131, retval = 0.000447644
pcm_likelihood x = -1.07131, retval = 0.000595836
pcm_likelihood x = 1.45057, retval = 8.99144e-05
pcm_likelihood x = -1.45057, retval = 0.000146111
pcm_likelihood x = 1.95078, retval = 6.6724e-06
pcm_likelihood x = -1.95078, retval = 1.49604e-05
pcm_likelihood x = 2.64003, retval = 9.19968e-08
pcm_likelihood x = -2.64003, retval = 3.42542e-07
pcm_likelihood x = 3.63137, retval = 6.36456e-11
pcm_likelihood x = -3.63137, retval = 5.20831e-10
pcm_likelihood x = 5.11992, retval = 1.44761e-16
pcm_likelihood x = -5.11992, retval = 4.03225e-15
pcm_likelihood x = 7.45666, retval = 2.02808e-27
pcm_likelihood x = -7.45666, retval = 3.9524e-25
pcm_likelihood x = 11.3023, retval = 1.84534e-50
pcm_likelihood x = -11.3023, retval = 8.89663e-47
pcm_likelihood x = 17.9641, retval = 1.41078e-105
pcm_likelihood x = -17.9641, retval = 1.76493e-99
pcm_likelihood x = 30.1781, retval = 1.08137e-256
pcm_likelihood x = -30.1781, retval = 3.6088e-246
pcm_likelihood x = 54.0388, retval = 0
pcm_likelihood x = -54.0388, retval = 0
pcm_likelihood x = 104.108, retval = 0
pcm_likelihood x = -104.108, retval = 0
pcm_likelihood x = 1.79769e+308, retval = -nan

Chain 1: Exception: gradient_of_f: The gradient of f is nan for parameter 0  (in 'modele66b29d3a6f9_pcm_bs_marg_nobeta_discourse' at line 71)


Can you elaborate a little bit more on what exactly happens in these first 2 evaluations? Where exactly does the gradient_of_f function come in?

pcm_likelihood x = 1.79769e+308, retval = -nan


Oh wow so it’s evaluating something for the gradient it isn’t evaluating for the main function. Can you turn your tolerances down and kick it off again to see if the behavior changes?

Also what version of rstan/cmdstanr?

1 Like

First we integrate the function itself here

Then for each input parameter, we do a separate integral. Code here

1 Like

Thank you! I just tried your suggestion, at while it didn’t work with if (retval == 0.0) {, I changed it to if (is_nan(retval)) { (which is what you meant, I guess?) and it finishes clean! That’s great! I’m just wondering, how can we distinguish between -nan and +nan, i.e. how can we make sure there aren’t fringe cases where this gives us bad values?

1 Like

Oh wait, that line:

pcm_likelihood x = 1.79769e+308, retval = -nan


I wonder why the first thing doesn’t blow up? It’s just the second has two nans (one at the top and one at the bottom).

So it just looks like evaluating this function at 1.8e308 gets us a NaN. Probably something like @nhuurre was saying.

1 Like

The first check is ignoring 0.0 values that lead to NaN gradients, which is something safe enough we do it in the code.

This is ignoring NaNs in values, which under the assumption that this is only happening out at crazy large values of the lpdf, seems fair.

You could make sure you aren’t dropping NaNs you care about (ones that actually indicate problems) with an extra condition like if(is_nan(retval) && abs(pars[mu_index]) > 1e20)

Or maybe it makes more sense to just bound these integrals to some large range like -1e4, 1e4. I’m not sure if that’s a bad idea with these quadratures or not. You could compute the full integral too. (Edit: Not sure last comment makes sense – you could like integrate [-1e5, 1e5] in generated quantities and compute against the [-1e4, 1e4] integral to make sure the cutoff is behaving fine)

1 Like

In the original example, I used sqrt(.Machine$double.eps)=1.490116e-08, but with sqrt(sqrt(.Machine$double.eps))=0.0001220703 the error persists.

I also wonder why

pcm_likelihood x = 1.79769e+308, retval = -nan
pcm_likelihood x = -1.79769e+308, retval = 0


it only happens for +1.8e308 and not negative that.

I changed the debug to

    sum_log_lik = sum(log_lik);
normal_lpdf_val = normal_lpdf(x | pars[mu_index], 1);
retval = exp(sum_log_lik + normal_lpdf_val);

print("pcm_likelihood x = ", x, ", sum_log_lik = ", sum_log_lik, ", normal_lpdf_val = ", normal_lpdf_val, ", retval = ", retval);

// if (is_nan(retval)) {
//   return 0.0;
// } else {
return retval;
// }


and get

Chain 1: mu_theta = -0.69882 ----------------------------------
pcm_likelihood x = 1.79769e+308, sum_log_lik = -nan, normal_lpdf_val = -inf, retval = -nan
pcm_likelihood x = -1.79769e+308, sum_log_lik = -inf, normal_lpdf_val = -inf, retval = 0
pcm_likelihood x = 0, sum_log_lik = -3.39651, normal_lpdf_val = -1.16311, retval = 0.010466
pcm_likelihood x = 3.08829, sum_log_lik = -17.2912, normal_lpdf_val = -8.09003, retval = 9.48558e-12
pcm_likelihood x = -3.08829, sum_log_lik = -5.84948, normal_lpdf_val = -3.77372, retval = 6.61757e-05
pcm_likelihood x = 148.993, sum_log_lik = -892.509, normal_lpdf_val = -11204.8, retval = 0
pcm_likelihood x = -148.993, sum_log_lik = -297.486, normal_lpdf_val = -10996.5, retval = 0
pcm_likelihood x = 3.41229e+06, sum_log_lik = -2.04737e+07, normal_lpdf_val = -5.82186e+12, retval = 0
pcm_likelihood x = -3.41229e+06, sum_log_lik = -6.82458e+06, normal_lpdf_val = -5.82186e+12, retval = 0
pcm_likelihood x = 2.06933e+18, sum_log_lik = -1.2416e+19, normal_lpdf_val = -2.14105e+36, retval = 0
pcm_likelihood x = -2.06933e+18, sum_log_lik = -4.13865e+18, normal_lpdf_val = -2.14105e+36, retval = 0
pcm_likelihood x = 2.087e+50, sum_log_lik = -1.2522e+51, normal_lpdf_val = -2.17779e+100, retval = 0
pcm_likelihood x = -2.087e+50, sum_log_lik = -4.174e+50, normal_lpdf_val = -2.17779e+100, retval = 0
pcm_likelihood x = 2.01977e+137, sum_log_lik = -1.21186e+138, normal_lpdf_val = -2.03973e+274, retval = 0
pcm_likelihood x = -2.01977e+137, sum_log_lik = -4.03953e+137, normal_lpdf_val = -2.03973e+274, retval = 0
pcm_likelihood x = 0.913049, sum_log_lik = -6.01139, normal_lpdf_val = -2.218, retval = 0.0002667
pcm_likelihood x = -0.913049, sum_log_lik = -2.90328, normal_lpdf_val = -0.941886, retval = 0.0213828
pcm_likelihood x = 14.1579, sum_log_lik = -83.4974, normal_lpdf_val = -111.28, retval = 2.56633e-85
pcm_likelihood x = -14.1579, sum_log_lik = -27.8158, normal_lpdf_val = -91.4923, retval = 1.53169e-52
pcm_likelihood x = 6704.22, sum_log_lik = -40223.8, normal_lpdf_val = -2.24779e+07, retval = 0
pcm_likelihood x = -6704.22, sum_log_lik = -13407.9, normal_lpdf_val = -2.24686e+07, retval = 0
pcm_likelihood x = 9.64173e+10, sum_log_lik = -5.78504e+11, normal_lpdf_val = -4.64814e+21, retval = 0
pcm_likelihood x = -9.64173e+10, sum_log_lik = -1.92835e+11, normal_lpdf_val = -4.64814e+21, retval = 0
pcm_likelihood x = 2.50895e+30, sum_log_lik = -1.50537e+31, normal_lpdf_val = -3.14742e+60, retval = 0
pcm_likelihood x = -2.50895e+30, sum_log_lik = -5.0179e+30, normal_lpdf_val = -3.14742e+60, retval = 0
pcm_likelihood x = 1.44726e+83, sum_log_lik = -8.68358e+83, normal_lpdf_val = -1.04729e+166, retval = 0
pcm_likelihood x = -1.44726e+83, sum_log_lik = -2.89453e+83, normal_lpdf_val = -1.04729e+166, retval = 0
pcm_likelihood x = 0.407298, sum_log_lik = -4.29732, normal_lpdf_val = -1.53069, retval = 0.00294395
pcm_likelihood x = -0.407298, sum_log_lik = -2.94183, normal_lpdf_val = -0.961431, retval = 0.0201761
pcm_likelihood x = 1.68207, sum_log_lik = -9.53533, normal_lpdf_val = -3.75325, retval = 1.69372e-06
pcm_likelihood x = -1.68207, sum_log_lik = -3.58283, normal_lpdf_val = -1.40233, retval = 0.00683872
pcm_likelihood x = 6.1509, sum_log_lik = -35.4651, normal_lpdf_val = -24.3783, retval = 1.02409e-26
pcm_likelihood x = -6.1509, sum_log_lik = -11.8098, normal_lpdf_val = -15.7815, retval = 1.04048e-12
pcm_likelihood x = 40.0396, sum_log_lik = -238.788, normal_lpdf_val = -830.729, retval = 0
pcm_likelihood x = -40.0396, sum_log_lik = -79.5792, normal_lpdf_val = -774.768, retval = 0
pcm_likelihood x = 792.92, sum_log_lik = -4756.07, normal_lpdf_val = -314916, retval = 0
pcm_likelihood x = -792.92, sum_log_lik = -1585.34, normal_lpdf_val = -313808, retval = 0
pcm_likelihood x = 0.198135, sum_log_lik = -3.77786, normal_lpdf_val = -1.3212, retval = 0.00610248
pcm_likelihood x = -0.198135, sum_log_lik = -3.12195, normal_lpdf_val = -1.04428, retval = 0.0155106
pcm_likelihood x = 0.640156, sum_log_lik = -5.01111, normal_lpdf_val = -1.81537, retval = 0.00108467
pcm_likelihood x = -0.640156, sum_log_lik = -2.86099, normal_lpdf_val = -0.920659, retval = 0.0227851
pcm_likelihood x = 1.24893, sum_log_lik = -7.44395, normal_lpdf_val = -2.8158, retval = 3.50143e-05
pcm_likelihood x = -1.24893, sum_log_lik = -3.11681, normal_lpdf_val = -1.07025, retval = 0.0151909
pcm_likelihood x = 2.26608, sum_log_lik = -12.6353, normal_lpdf_val = -5.31426, retval = 1.60174e-08
pcm_likelihood x = -2.26608, sum_log_lik = -4.42893, normal_lpdf_val = -2.14709, retval = 0.00139338
pcm_likelihood x = 4.29646, sum_log_lik = -24.3914, normal_lpdf_val = -13.3954, retval = 3.8854e-17
pcm_likelihood x = -4.29646, sum_log_lik = -8.14432, normal_lpdf_val = -7.39046, retval = 1.79198e-07
pcm_likelihood x = 9.13029, sum_log_lik = -53.3322, normal_lpdf_val = -49.2246, retval = 2.88478e-45
pcm_likelihood x = -9.13029, sum_log_lik = -17.761, normal_lpdf_val = -36.4638, retval = 2.82149e-24
pcm_likelihood x = 23.1111, sum_log_lik = -137.216, normal_lpdf_val = -284.375, retval = 8.04202e-184
pcm_likelihood x = -23.1111, sum_log_lik = -45.7222, normal_lpdf_val = -252.074, retval = 4.66619e-130
pcm_likelihood x = 74.2771, sum_log_lik = -444.212, normal_lpdf_val = -2811.61, retval = 0
pcm_likelihood x = -74.2771, sum_log_lik = -148.054, normal_lpdf_val = -2707.8, retval = 0
pcm_likelihood x = 326.721, sum_log_lik = -1958.88, normal_lpdf_val = -53602.8, retval = 0
pcm_likelihood x = -326.721, sum_log_lik = -652.942, normal_lpdf_val = -53146.1, retval = 0
pcm_likelihood x = 0.0983968, sum_log_lik = -3.57237, normal_lpdf_val = -1.23672, retval = 0.00815531
pcm_likelihood x = -0.0983968, sum_log_lik = -3.24706, normal_lpdf_val = -1.09919, retval = 0.0129553
pcm_likelihood x = 0.300606, sum_log_lik = -4.01749, normal_lpdf_val = -1.41836, retval = 0.0043575
pcm_likelihood x = -0.300606, sum_log_lik = -3.02019, normal_lpdf_val = -0.998226, retval = 0.0179815
pcm_likelihood x = 0.519858, sum_log_lik = -4.62517, normal_lpdf_val = -1.66153, retval = 0.00186089
pcm_likelihood x = -0.519858, sum_log_lik = -2.888, normal_lpdf_val = -0.934952, retval = 0.0218632
pcm_likelihood x = 0.770362, sum_log_lik = -5.46785, normal_lpdf_val = -1.99819, retval = 0.000572194
pcm_likelihood x = -0.770362, sum_log_lik = -2.86439, normal_lpdf_val = -0.921498, retval = 0.0226887
pcm_likelihood x = 1.07131, sum_log_lik = -6.66177, normal_lpdf_val = -2.48562, retval = 0.000106498
pcm_likelihood x = -1.07131, sum_log_lik = -2.98447, normal_lpdf_val = -0.988314, retval = 0.018821
pcm_likelihood x = 1.45057, sum_log_lik = -8.38897, normal_lpdf_val = -3.22888, retval = 9.00394e-06
pcm_likelihood x = -1.45057, sum_log_lik = -3.31156, normal_lpdf_val = -1.2015, retval = 0.0109648
pcm_likelihood x = 1.95078, sum_log_lik = -10.9309, normal_lpdf_val = -4.42912, retval = 2.13419e-07
pcm_likelihood x = -1.95078, sum_log_lik = -3.9481, normal_lpdf_val = -1.70264, retval = 0.00351491
pcm_likelihood x = 2.64003, sum_log_lik = -14.7237, normal_lpdf_val = -6.4929, retval = 6.10587e-10
pcm_likelihood x = -2.64003, sum_log_lik = -5.05181, normal_lpdf_val = -2.80309, retval = 0.000387848
pcm_likelihood x = 3.63137, sum_log_lik = -20.4604, normal_lpdf_val = -10.2942, retval = 4.39979e-14
pcm_likelihood x = -3.63137, sum_log_lik = -6.86289, normal_lpdf_val = -5.21887, retval = 5.66182e-06
pcm_likelihood x = 5.11992, sum_log_lik = -29.2969, normal_lpdf_val = -17.8478, retval = 3.35219e-21
pcm_likelihood x = -5.11992, sum_log_lik = -9.76236, normal_lpdf_val = -10.692, retval = 1.30856e-09
pcm_likelihood x = 7.45666, sum_log_lik = -43.2926, normal_lpdf_val = -34.1749, retval = 2.27144e-34
pcm_likelihood x = -7.45666, sum_log_lik = -14.4155, normal_lpdf_val = -23.7531, retval = 2.65197e-17
pcm_likelihood x = 11.3023, sum_log_lik = -66.3636, normal_lpdf_val = -72.9319, retval = 3.19682e-61
pcm_likelihood x = -11.3023, sum_log_lik = -22.1046, normal_lpdf_val = -57.1354, retval = 3.85929e-35
pcm_likelihood x = 17.9641, sum_log_lik = -106.335, normal_lpdf_val = -175.071, retval = 6.1225e-123
pcm_likelihood x = -17.9641, sum_log_lik = -35.4282, normal_lpdf_val = -149.964, retval = 3.05619e-81
pcm_likelihood x = 30.1781, sum_log_lik = -179.619, normal_lpdf_val = -477.611, retval = 3.70415e-286
pcm_likelihood x = -30.1781, sum_log_lik = -59.8562, normal_lpdf_val = -435.433, retval = 7.91719e-216
pcm_likelihood x = 54.0388, sum_log_lik = -322.783, normal_lpdf_val = -1499.02, retval = 0
pcm_likelihood x = -54.0388, sum_log_lik = -107.578, normal_lpdf_val = -1423.49, retval = 0
pcm_likelihood x = 104.108, sum_log_lik = -623.196, normal_lpdf_val = -5493.13, retval = 0
pcm_likelihood x = -104.108, sum_log_lik = -207.715, normal_lpdf_val = -5347.62, retval = 0
pcm_likelihood x = 0.0491151, sum_log_lik = -3.48097, normal_lpdf_val = -1.19864, retval = 0.00928264
pcm_likelihood x = -0.0491151, sum_log_lik = -3.31864, normal_lpdf_val = -1.13, retval = 0.0116945
pcm_likelihood x = 0.148013, sum_log_lik = -3.67116, normal_lpdf_val = -1.2775, retval = 0.00709289
pcm_likelihood x = -0.148013, sum_log_lik = -3.18155, normal_lpdf_val = -1.07063, retval = 0.0142331
pcm_likelihood x = 0.248939, sum_log_lik = -3.89307, normal_lpdf_val = -1.36806, retval = 0.00518944
pcm_likelihood x = -0.248939, sum_log_lik = -3.06817, normal_lpdf_val = -1.02013, retval = 0.0167676
pcm_likelihood x = 0.353325, sum_log_lik = -4.15194, normal_lpdf_val = -1.47244, retval = 0.0036088
pcm_likelihood x = -0.353325, sum_log_lik = -2.97804, normal_lpdf_val = -0.978622, retval = 0.0191269
pcm_likelihood x = 0.462734, sum_log_lik = -4.45467, normal_lpdf_val = -1.59354, retval = 0.00236208
pcm_likelihood x = -0.462734, sum_log_lik = -2.91173, normal_lpdf_val = -0.946807, retval = 0.0210988
pcm_likelihood x = 0.578912, sum_log_lik = -4.81016, normal_lpdf_val = -1.73524, retval = 0.00143672
pcm_likelihood x = -0.578912, sum_log_lik = -2.87095, normal_lpdf_val = -0.926127, retval = 0.0224363
pcm_likelihood x = 0.70387, sum_log_lik = -5.22971, normal_lpdf_val = -1.90271, retval = 0.000798783
pcm_likelihood x = -0.70387, sum_log_lik = -2.85861, normal_lpdf_val = -0.918951, retval = 0.0228784
pcm_likelihood x = 0.839966, sum_log_lik = -5.72762, normal_lpdf_val = -2.10287, retval = 0.000397432
pcm_likelihood x = -0.839966, sum_log_lik = -2.87902, normal_lpdf_val = -0.9289, retval = 0.0221943
pcm_likelihood x = 0.990015, sum_log_lik = -6.32179, normal_lpdf_val = -2.34502, retval = 0.000172208
pcm_likelihood x = -0.990015, sum_log_lik = -2.93809, normal_lpdf_val = -0.961336, retval = 0.0202536
pcm_likelihood x = 1.15743, sum_log_lik = -7.0346, normal_lpdf_val = -2.64177, retval = 6.27484e-05
pcm_likelihood x = -1.15743, sum_log_lik = -3.0436, normal_lpdf_val = -1.0241, retval = 0.0171166
pcm_likelihood x = 1.34641, sum_log_lik = -7.8939, normal_lpdf_val = -3.01043, retval = 1.83786e-05
pcm_likelihood x = -1.34641, sum_log_lik = -3.20558, normal_lpdf_val = -1.12863, retval = 0.0131123
pcm_likelihood x = 1.56217, sum_log_lik = -8.93424, normal_lpdf_val = -3.47497, retval = 4.08084e-06
pcm_likelihood x = -1.56217, sum_log_lik = -3.43662, normal_lpdf_val = -1.29162, retval = 0.00884196
pcm_likelihood x = 1.81124, sum_log_lik = -10.1985, normal_lpdf_val = -4.06914, retval = 6.36248e-07
pcm_likelihood x = -1.81124, sum_log_lik = -3.75247, normal_lpdf_val = -1.53768, retval = 0.00504101
pcm_likelihood x = 2.10192, sum_log_lik = -11.7402, normal_lpdf_val = -4.84102, retval = 6.29302e-08
pcm_likelihood x = -2.10192, sum_log_lik = -4.17254, normal_lpdf_val = -1.90329, retval = 0.00229773
pcm_likelihood x = 2.44484, sum_log_lik = -13.6261, normal_lpdf_val = -5.86025, retval = 3.44492e-09
pcm_likelihood x = -2.44484, sum_log_lik = -4.72073, normal_lpdf_val = -2.44324, retval = 0.000773978
pcm_likelihood x = 2.85372, sum_log_lik = -15.9406, normal_lpdf_val = -7.22921, retval = 8.65883e-11
pcm_likelihood x = -2.85372, sum_log_lik = -5.42646, normal_lpdf_val = -3.24074, retval = 0.000172141
pcm_likelihood x = 3.34646, sum_log_lik = -18.7916, normal_lpdf_val = -9.10108, retval = 7.6975e-13
pcm_likelihood x = -3.34646, sum_log_lik = -6.32627, normal_lpdf_val = -4.42394, retval = 2.14411e-05
pcm_likelihood x = 3.94665, sum_log_lik = -22.3188, normal_lpdf_val = -11.7091, retval = 1.66669e-15
pcm_likelihood x = -3.94665, sum_log_lik = -7.46628, normal_lpdf_val = -6.19313, retval = 1.16895e-06
pcm_likelihood x = 4.68567, sum_log_lik = -26.7064, normal_lpdf_val = -15.4153, retval = 5.09067e-19
pcm_likelihood x = -4.68567, sum_log_lik = -8.90614, normal_lpdf_val = -8.86644, retval = 1.9119e-08
pcm_likelihood x = 5.60576, sum_log_lik = -32.2014, normal_lpdf_val = -20.7928, retval = 9.6582e-24
pcm_likelihood x = -5.60576, sum_log_lik = -10.7254, normal_lpdf_val = -12.958, retval = 5.18141e-11
pcm_likelihood x = 6.76433, sum_log_lik = -39.1413, normal_lpdf_val = -28.7683, retval = 3.21561e-30
pcm_likelihood x = -6.76433, sum_log_lik = -13.033, normal_lpdf_val = -19.3142, retval = 8.94956e-15
pcm_likelihood x = 8.24038, sum_log_lik = -47.9935, normal_lpdf_val = -40.8736, retval = 2.5439e-39
pcm_likelihood x = -8.24038, sum_log_lik = -15.9818, normal_lpdf_val = -29.3565, retval = 2.04095e-20
pcm_likelihood x = 10.1439, sum_log_lik = -59.4138, normal_lpdf_val = -59.7017, retval = 1.85686e-52
pcm_likelihood x = -10.1439, sum_log_lik = -19.788, normal_lpdf_val = -45.5241, retval = 4.31808e-29
pcm_likelihood x = 12.6302, sum_log_lik = -74.3315, normal_lpdf_val = -89.751, retval = 5.49412e-72
pcm_likelihood x = -12.6302, sum_log_lik = -24.7605, normal_lpdf_val = -72.0984, retval = 8.60405e-43
pcm_likelihood x = 15.9213, sum_log_lik = -94.0778, normal_lpdf_val = -139.033, retval = 5.76994e-102
pcm_likelihood x = -15.9213, sum_log_lik = -31.3426, normal_lpdf_val = -116.781, retval = 4.68549e-65
pcm_likelihood x = 20.3392, sum_log_lik = -120.585, normal_lpdf_val = -222.218, retval = 1.32498e-149
pcm_likelihood x = -20.3392, sum_log_lik = -40.1784, normal_lpdf_val = -193.792, retval = 2.44408e-102
pcm_likelihood x = 26.3585, sum_log_lik = -156.701, normal_lpdf_val = -366.967, retval = 3.74853e-228
pcm_likelihood x = -26.3585, sum_log_lik = -52.2169, normal_lpdf_val = -330.128, retval = 8.90986e-167
pcm_likelihood x = 34.6893, sum_log_lik = -206.686, normal_lpdf_val = -627.077, retval = 0
pcm_likelihood x = -34.6893, sum_log_lik = -68.8785, normal_lpdf_val = -578.594, retval = 6.40067e-282
pcm_likelihood x = 46.4129, sum_log_lik = -277.027, normal_lpdf_val = -1110.68, retval = 0
pcm_likelihood x = -46.4129, sum_log_lik = -92.3258, normal_lpdf_val = -1045.81, retval = 0
pcm_likelihood x = 63.2055, sum_log_lik = -377.783, normal_lpdf_val = -2042.8, retval = 0
pcm_likelihood x = -63.2055, sum_log_lik = -125.911, normal_lpdf_val = -1954.46, retval = 0
pcm_likelihood x = 87.715, sum_log_lik = -524.84, normal_lpdf_val = -3909.42, retval = 0
pcm_likelihood x = -87.715, sum_log_lik = -174.93, normal_lpdf_val = -3786.82, retval = 0
pcm_likelihood x = 124.21, sum_log_lik = -743.808, normal_lpdf_val = -7801.99, retval = 0
pcm_likelihood x = -124.21, sum_log_lik = -247.919, normal_lpdf_val = -7628.39, retval = 0
marg_lik = 0.0340928 ----------------------------------

Chain 1: mu_theta = -0.69882 ----------------------------------
pcm_likelihood x = 1.79769e+308, sum_log_lik = -nan, normal_lpdf_val = -inf, retval = -nan
pcm_likelihood x = -1.79769e+308, sum_log_lik = -inf, normal_lpdf_val = -inf, retval = 0
pcm_likelihood x = 0, sum_log_lik = -3.39651, normal_lpdf_val = -1.16311, retval = 0.010466
pcm_likelihood x = 3.08829, sum_log_lik = -17.2912, normal_lpdf_val = -8.09003, retval = 9.48558e-12
pcm_likelihood x = -3.08829, sum_log_lik = -5.84948, normal_lpdf_val = -3.77372, retval = 6.61757e-05
pcm_likelihood x = 148.993, sum_log_lik = -892.509, normal_lpdf_val = -11204.8, retval = 0
pcm_likelihood x = -148.993, sum_log_lik = -297.486, normal_lpdf_val = -10996.5, retval = 0
pcm_likelihood x = 3.41229e+06, sum_log_lik = -2.04737e+07, normal_lpdf_val = -5.82186e+12, retval = 0
pcm_likelihood x = -3.41229e+06, sum_log_lik = -6.82458e+06, normal_lpdf_val = -5.82186e+12, retval = 0
pcm_likelihood x = 2.06933e+18, sum_log_lik = -1.2416e+19, normal_lpdf_val = -2.14105e+36, retval = 0
pcm_likelihood x = -2.06933e+18, sum_log_lik = -4.13865e+18, normal_lpdf_val = -2.14105e+36, retval = 0
pcm_likelihood x = 2.087e+50, sum_log_lik = -1.2522e+51, normal_lpdf_val = -2.17779e+100, retval = 0
pcm_likelihood x = -2.087e+50, sum_log_lik = -4.174e+50, normal_lpdf_val = -2.17779e+100, retval = 0
pcm_likelihood x = 2.01977e+137, sum_log_lik = -1.21186e+138, normal_lpdf_val = -2.03973e+274, retval = 0
pcm_likelihood x = -2.01977e+137, sum_log_lik = -4.03953e+137, normal_lpdf_val = -2.03973e+274, retval = 0
pcm_likelihood x = 0.913049, sum_log_lik = -6.01139, normal_lpdf_val = -2.218, retval = 0.0002667
pcm_likelihood x = -0.913049, sum_log_lik = -2.90328, normal_lpdf_val = -0.941886, retval = 0.0213828
pcm_likelihood x = 14.1579, sum_log_lik = -83.4974, normal_lpdf_val = -111.28, retval = 2.56633e-85
pcm_likelihood x = -14.1579, sum_log_lik = -27.8158, normal_lpdf_val = -91.4923, retval = 1.53169e-52
pcm_likelihood x = 6704.22, sum_log_lik = -40223.8, normal_lpdf_val = -2.24779e+07, retval = 0
pcm_likelihood x = -6704.22, sum_log_lik = -13407.9, normal_lpdf_val = -2.24686e+07, retval = 0
pcm_likelihood x = 9.64173e+10, sum_log_lik = -5.78504e+11, normal_lpdf_val = -4.64814e+21, retval = 0
pcm_likelihood x = -9.64173e+10, sum_log_lik = -1.92835e+11, normal_lpdf_val = -4.64814e+21, retval = 0
pcm_likelihood x = 2.50895e+30, sum_log_lik = -1.50537e+31, normal_lpdf_val = -3.14742e+60, retval = 0
pcm_likelihood x = -2.50895e+30, sum_log_lik = -5.0179e+30, normal_lpdf_val = -3.14742e+60, retval = 0
pcm_likelihood x = 1.44726e+83, sum_log_lik = -8.68358e+83, normal_lpdf_val = -1.04729e+166, retval = 0
pcm_likelihood x = -1.44726e+83, sum_log_lik = -2.89453e+83, normal_lpdf_val = -1.04729e+166, retval = 0
pcm_likelihood x = 0.407298, sum_log_lik = -4.29732, normal_lpdf_val = -1.53069, retval = 0.00294395
pcm_likelihood x = -0.407298, sum_log_lik = -2.94183, normal_lpdf_val = -0.961431, retval = 0.0201761
pcm_likelihood x = 1.68207, sum_log_lik = -9.53533, normal_lpdf_val = -3.75325, retval = 1.69372e-06
pcm_likelihood x = -1.68207, sum_log_lik = -3.58283, normal_lpdf_val = -1.40233, retval = 0.00683872
pcm_likelihood x = 6.1509, sum_log_lik = -35.4651, normal_lpdf_val = -24.3783, retval = 1.02409e-26
pcm_likelihood x = -6.1509, sum_log_lik = -11.8098, normal_lpdf_val = -15.7815, retval = 1.04048e-12
pcm_likelihood x = 40.0396, sum_log_lik = -238.788, normal_lpdf_val = -830.729, retval = 0
pcm_likelihood x = -40.0396, sum_log_lik = -79.5792, normal_lpdf_val = -774.768, retval = 0
pcm_likelihood x = 792.92, sum_log_lik = -4756.07, normal_lpdf_val = -314916, retval = 0
pcm_likelihood x = -792.92, sum_log_lik = -1585.34, normal_lpdf_val = -313808, retval = 0
pcm_likelihood x = 0.198135, sum_log_lik = -3.77786, normal_lpdf_val = -1.3212, retval = 0.00610248
pcm_likelihood x = -0.198135, sum_log_lik = -3.12195, normal_lpdf_val = -1.04428, retval = 0.0155106
pcm_likelihood x = 0.640156, sum_log_lik = -5.01111, normal_lpdf_val = -1.81537, retval = 0.00108467
pcm_likelihood x = -0.640156, sum_log_lik = -2.86099, normal_lpdf_val = -0.920659, retval = 0.0227851
pcm_likelihood x = 1.24893, sum_log_lik = -7.44395, normal_lpdf_val = -2.8158, retval = 3.50143e-05
pcm_likelihood x = -1.24893, sum_log_lik = -3.11681, normal_lpdf_val = -1.07025, retval = 0.0151909
pcm_likelihood x = 2.26608, sum_log_lik = -12.6353, normal_lpdf_val = -5.31426, retval = 1.60174e-08
pcm_likelihood x = -2.26608, sum_log_lik = -4.42893, normal_lpdf_val = -2.14709, retval = 0.00139338
pcm_likelihood x = 4.29646, sum_log_lik = -24.3914, normal_lpdf_val = -13.3954, retval = 3.8854e-17
pcm_likelihood x = -4.29646, sum_log_lik = -8.14432, normal_lpdf_val = -7.39046, retval = 1.79198e-07
pcm_likelihood x = 9.13029, sum_log_lik = -53.3322, normal_lpdf_val = -49.2246, retval = 2.88478e-45
pcm_likelihood x = -9.13029, sum_log_lik = -17.761, normal_lpdf_val = -36.4638, retval = 2.82149e-24
pcm_likelihood x = 23.1111, sum_log_lik = -137.216, normal_lpdf_val = -284.375, retval = 8.04202e-184
pcm_likelihood x = -23.1111, sum_log_lik = -45.7222, normal_lpdf_val = -252.074, retval = 4.66619e-130
pcm_likelihood x = 74.2771, sum_log_lik = -444.212, normal_lpdf_val = -2811.61, retval = 0
pcm_likelihood x = -74.2771, sum_log_lik = -148.054, normal_lpdf_val = -2707.8, retval = 0
pcm_likelihood x = 326.721, sum_log_lik = -1958.88, normal_lpdf_val = -53602.8, retval = 0
pcm_likelihood x = -326.721, sum_log_lik = -652.942, normal_lpdf_val = -53146.1, retval = 0
pcm_likelihood x = 0.0983968, sum_log_lik = -3.57237, normal_lpdf_val = -1.23672, retval = 0.00815531
pcm_likelihood x = -0.0983968, sum_log_lik = -3.24706, normal_lpdf_val = -1.09919, retval = 0.0129553
pcm_likelihood x = 0.300606, sum_log_lik = -4.01749, normal_lpdf_val = -1.41836, retval = 0.0043575
pcm_likelihood x = -0.300606, sum_log_lik = -3.02019, normal_lpdf_val = -0.998226, retval = 0.0179815
pcm_likelihood x = 0.519858, sum_log_lik = -4.62517, normal_lpdf_val = -1.66153, retval = 0.00186089
pcm_likelihood x = -0.519858, sum_log_lik = -2.888, normal_lpdf_val = -0.934952, retval = 0.0218632
pcm_likelihood x = 0.770362, sum_log_lik = -5.46785, normal_lpdf_val = -1.99819, retval = 0.000572194
pcm_likelihood x = -0.770362, sum_log_lik = -2.86439, normal_lpdf_val = -0.921498, retval = 0.0226887
pcm_likelihood x = 1.07131, sum_log_lik = -6.66177, normal_lpdf_val = -2.48562, retval = 0.000106498
pcm_likelihood x = -1.07131, sum_log_lik = -2.98447, normal_lpdf_val = -0.988314, retval = 0.018821
pcm_likelihood x = 1.45057, sum_log_lik = -8.38897, normal_lpdf_val = -3.22888, retval = 9.00394e-06
pcm_likelihood x = -1.45057, sum_log_lik = -3.31156, normal_lpdf_val = -1.2015, retval = 0.0109648
pcm_likelihood x = 1.95078, sum_log_lik = -10.9309, normal_lpdf_val = -4.42912, retval = 2.13419e-07
pcm_likelihood x = -1.95078, sum_log_lik = -3.9481, normal_lpdf_val = -1.70264, retval = 0.00351491
pcm_likelihood x = 2.64003, sum_log_lik = -14.7237, normal_lpdf_val = -6.4929, retval = 6.10587e-10
pcm_likelihood x = -2.64003, sum_log_lik = -5.05181, normal_lpdf_val = -2.80309, retval = 0.000387848
pcm_likelihood x = 3.63137, sum_log_lik = -20.4604, normal_lpdf_val = -10.2942, retval = 4.39979e-14
pcm_likelihood x = -3.63137, sum_log_lik = -6.86289, normal_lpdf_val = -5.21887, retval = 5.66182e-06
pcm_likelihood x = 5.11992, sum_log_lik = -29.2969, normal_lpdf_val = -17.8478, retval = 3.35219e-21
pcm_likelihood x = -5.11992, sum_log_lik = -9.76236, normal_lpdf_val = -10.692, retval = 1.30856e-09
pcm_likelihood x = 7.45666, sum_log_lik = -43.2926, normal_lpdf_val = -34.1749, retval = 2.27144e-34
pcm_likelihood x = -7.45666, sum_log_lik = -14.4155, normal_lpdf_val = -23.7531, retval = 2.65197e-17
pcm_likelihood x = 11.3023, sum_log_lik = -66.3636, normal_lpdf_val = -72.9319, retval = 3.19682e-61
pcm_likelihood x = -11.3023, sum_log_lik = -22.1046, normal_lpdf_val = -57.1354, retval = 3.85929e-35
pcm_likelihood x = 17.9641, sum_log_lik = -106.335, normal_lpdf_val = -175.071, retval = 6.1225e-123
pcm_likelihood x = -17.9641, sum_log_lik = -35.4282, normal_lpdf_val = -149.964, retval = 3.05619e-81
pcm_likelihood x = 30.1781, sum_log_lik = -179.619, normal_lpdf_val = -477.611, retval = 3.70415e-286
pcm_likelihood x = -30.1781, sum_log_lik = -59.8562, normal_lpdf_val = -435.433, retval = 7.91719e-216
pcm_likelihood x = 54.0388, sum_log_lik = -322.783, normal_lpdf_val = -1499.02, retval = 0
pcm_likelihood x = -54.0388, sum_log_lik = -107.578, normal_lpdf_val = -1423.49, retval = 0
pcm_likelihood x = 104.108, sum_log_lik = -623.196, normal_lpdf_val = -5493.13, retval = 0
pcm_likelihood x = -104.108, sum_log_lik = -207.715, normal_lpdf_val = -5347.62, retval = 0
pcm_likelihood x = 0.0491151, sum_log_lik = -3.48097, normal_lpdf_val = -1.19864, retval = 0.00928264
pcm_likelihood x = -0.0491151, sum_log_lik = -3.31864, normal_lpdf_val = -1.13, retval = 0.0116945
pcm_likelihood x = 0.148013, sum_log_lik = -3.67116, normal_lpdf_val = -1.2775, retval = 0.00709289
pcm_likelihood x = -0.148013, sum_log_lik = -3.18155, normal_lpdf_val = -1.07063, retval = 0.0142331
pcm_likelihood x = 0.248939, sum_log_lik = -3.89307, normal_lpdf_val = -1.36806, retval = 0.00518944
pcm_likelihood x = -0.248939, sum_log_lik = -3.06817, normal_lpdf_val = -1.02013, retval = 0.0167676
pcm_likelihood x = 0.353325, sum_log_lik = -4.15194, normal_lpdf_val = -1.47244, retval = 0.0036088
pcm_likelihood x = -0.353325, sum_log_lik = -2.97804, normal_lpdf_val = -0.978622, retval = 0.0191269
pcm_likelihood x = 0.462734, sum_log_lik = -4.45467, normal_lpdf_val = -1.59354, retval = 0.00236208
pcm_likelihood x = -0.462734, sum_log_lik = -2.91173, normal_lpdf_val = -0.946807, retval = 0.0210988
pcm_likelihood x = 0.578912, sum_log_lik = -4.81016, normal_lpdf_val = -1.73524, retval = 0.00143672
pcm_likelihood x = -0.578912, sum_log_lik = -2.87095, normal_lpdf_val = -0.926127, retval = 0.0224363
pcm_likelihood x = 0.70387, sum_log_lik = -5.22971, normal_lpdf_val = -1.90271, retval = 0.000798783
pcm_likelihood x = -0.70387, sum_log_lik = -2.85861, normal_lpdf_val = -0.918951, retval = 0.0228784
pcm_likelihood x = 0.839966, sum_log_lik = -5.72762, normal_lpdf_val = -2.10287, retval = 0.000397432
pcm_likelihood x = -0.839966, sum_log_lik = -2.87902, normal_lpdf_val = -0.9289, retval = 0.0221943
pcm_likelihood x = 0.990015, sum_log_lik = -6.32179, normal_lpdf_val = -2.34502, retval = 0.000172208
pcm_likelihood x = -0.990015, sum_log_lik = -2.93809, normal_lpdf_val = -0.961336, retval = 0.0202536
pcm_likelihood x = 1.15743, sum_log_lik = -7.0346, normal_lpdf_val = -2.64177, retval = 6.27484e-05
pcm_likelihood x = -1.15743, sum_log_lik = -3.0436, normal_lpdf_val = -1.0241, retval = 0.0171166
pcm_likelihood x = 1.34641, sum_log_lik = -7.8939, normal_lpdf_val = -3.01043, retval = 1.83786e-05
pcm_likelihood x = -1.34641, sum_log_lik = -3.20558, normal_lpdf_val = -1.12863, retval = 0.0131123
pcm_likelihood x = 1.56217, sum_log_lik = -8.93424, normal_lpdf_val = -3.47497, retval = 4.08084e-06
pcm_likelihood x = -1.56217, sum_log_lik = -3.43662, normal_lpdf_val = -1.29162, retval = 0.00884196
pcm_likelihood x = 1.81124, sum_log_lik = -10.1985, normal_lpdf_val = -4.06914, retval = 6.36248e-07
pcm_likelihood x = -1.81124, sum_log_lik = -3.75247, normal_lpdf_val = -1.53768, retval = 0.00504101
pcm_likelihood x = 2.10192, sum_log_lik = -11.7402, normal_lpdf_val = -4.84102, retval = 6.29302e-08
pcm_likelihood x = -2.10192, sum_log_lik = -4.17254, normal_lpdf_val = -1.90329, retval = 0.00229773
pcm_likelihood x = 2.44484, sum_log_lik = -13.6261, normal_lpdf_val = -5.86025, retval = 3.44492e-09
pcm_likelihood x = -2.44484, sum_log_lik = -4.72073, normal_lpdf_val = -2.44324, retval = 0.000773978
pcm_likelihood x = 2.85372, sum_log_lik = -15.9406, normal_lpdf_val = -7.22921, retval = 8.65883e-11
pcm_likelihood x = -2.85372, sum_log_lik = -5.42646, normal_lpdf_val = -3.24074, retval = 0.000172141
pcm_likelihood x = 3.34646, sum_log_lik = -18.7916, normal_lpdf_val = -9.10108, retval = 7.6975e-13
pcm_likelihood x = -3.34646, sum_log_lik = -6.32627, normal_lpdf_val = -4.42394, retval = 2.14411e-05
pcm_likelihood x = 3.94665, sum_log_lik = -22.3188, normal_lpdf_val = -11.7091, retval = 1.66669e-15
pcm_likelihood x = -3.94665, sum_log_lik = -7.46628, normal_lpdf_val = -6.19313, retval = 1.16895e-06
pcm_likelihood x = 4.68567, sum_log_lik = -26.7064, normal_lpdf_val = -15.4153, retval = 5.09067e-19
pcm_likelihood x = -4.68567, sum_log_lik = -8.90614, normal_lpdf_val = -8.86644, retval = 1.9119e-08
pcm_likelihood x = 5.60576, sum_log_lik = -32.2014, normal_lpdf_val = -20.7928, retval = 9.6582e-24
pcm_likelihood x = -5.60576, sum_log_lik = -10.7254, normal_lpdf_val = -12.958, retval = 5.18141e-11
pcm_likelihood x = 6.76433, sum_log_lik = -39.1413, normal_lpdf_val = -28.7683, retval = 3.21561e-30
pcm_likelihood x = -6.76433, sum_log_lik = -13.033, normal_lpdf_val = -19.3142, retval = 8.94956e-15
pcm_likelihood x = 8.24038, sum_log_lik = -47.9935, normal_lpdf_val = -40.8736, retval = 2.5439e-39
pcm_likelihood x = -8.24038, sum_log_lik = -15.9818, normal_lpdf_val = -29.3565, retval = 2.04095e-20
pcm_likelihood x = 10.1439, sum_log_lik = -59.4138, normal_lpdf_val = -59.7017, retval = 1.85686e-52
pcm_likelihood x = -10.1439, sum_log_lik = -19.788, normal_lpdf_val = -45.5241, retval = 4.31808e-29
pcm_likelihood x = 12.6302, sum_log_lik = -74.3315, normal_lpdf_val = -89.751, retval = 5.49412e-72
pcm_likelihood x = -12.6302, sum_log_lik = -24.7605, normal_lpdf_val = -72.0984, retval = 8.60405e-43
pcm_likelihood x = 15.9213, sum_log_lik = -94.0778, normal_lpdf_val = -139.033, retval = 5.76994e-102
pcm_likelihood x = -15.9213, sum_log_lik = -31.3426, normal_lpdf_val = -116.781, retval = 4.68549e-65
pcm_likelihood x = 20.3392, sum_log_lik = -120.585, normal_lpdf_val = -222.218, retval = 1.32498e-149
pcm_likelihood x = -20.3392, sum_log_lik = -40.1784, normal_lpdf_val = -193.792, retval = 2.44408e-102
pcm_likelihood x = 26.3585, sum_log_lik = -156.701, normal_lpdf_val = -366.967, retval = 3.74853e-228
pcm_likelihood x = -26.3585, sum_log_lik = -52.2169, normal_lpdf_val = -330.128, retval = 8.90986e-167
pcm_likelihood x = 34.6893, sum_log_lik = -206.686, normal_lpdf_val = -627.077, retval = 0
pcm_likelihood x = -34.6893, sum_log_lik = -68.8785, normal_lpdf_val = -578.594, retval = 6.40067e-282
pcm_likelihood x = 46.4129, sum_log_lik = -277.027, normal_lpdf_val = -1110.68, retval = 0
pcm_likelihood x = -46.4129, sum_log_lik = -92.3258, normal_lpdf_val = -1045.81, retval = 0
pcm_likelihood x = 63.2055, sum_log_lik = -377.783, normal_lpdf_val = -2042.8, retval = 0
pcm_likelihood x = -63.2055, sum_log_lik = -125.911, normal_lpdf_val = -1954.46, retval = 0
pcm_likelihood x = 87.715, sum_log_lik = -524.84, normal_lpdf_val = -3909.42, retval = 0
pcm_likelihood x = -87.715, sum_log_lik = -174.93, normal_lpdf_val = -3786.82, retval = 0
pcm_likelihood x = 124.21, sum_log_lik = -743.808, normal_lpdf_val = -7801.99, retval = 0
pcm_likelihood x = -124.21, sum_log_lik = -247.919, normal_lpdf_val = -7628.39, retval = 0
pcm_likelihood x = 1.79769e+308, sum_log_lik = -nan, normal_lpdf_val = -inf, retval = -nan

Chain 1: Exception: gradient_of_f: The gradient of f is nan for parameter 0  (in 'modele66b68dbfe35_pcm_bs_marg_nobeta_discourse' at line 81)


So I guess the nan comes from the log_softmax etc. rather than the “correct” -inf of the normal_lpdf, right?

I am using rstan_2.21.2 by the way. I know you guys have updated the integrator in the latest version of CmdStan, but I haven’t tried it yet (since on remote I am dependent on my sysadmin for updates and locally, I have Big Sur → Big Trouble with everything regarding compilation)

1 Like

Oh! 1.79769e+308+1.79769e+308 overflows so cumulative_sum() produces infinities and log_softmax() is all nan.

The exp(sum(log_lik)) part that comes from softmax is always \leq 1 so x sufficiently far in the tails of the normal_lpdf will always yield zero total. You can put this at the start of the function

if (abs(x - pars[mu_index]) > 1000) {
return 0.0;
}


Similar idea should work with categorical_logit: check at the start if the value of x is so far from reasonable that no calculation is needed to see that the result is zero and there’s no need to worry about infinities.

3 Likes

Even more debug:

Chain 1: mu_theta = 1.90732 ----------------------------------
pcm_likelihood x = 1.79769e+308, log_lik = [-inf,-nan], normal_lpdf_val = -inf, retval = -nan
pcm_likelihood x = -1.79769e+308, log_lik = [0,-inf], normal_lpdf_val = -inf, retval = 0
[...]
pcm_likelihood x = 104.108, log_lik = [-416.331,-206.865], normal_lpdf_val = -5223.38, retval = 0
pcm_likelihood x = -104.108, log_lik = [0,-207.715], normal_lpdf_val = -5620.51, retval = 0
marg_lik = 0.0029028 ----------------------------------

Chain 1: mu_theta = 1.90732 ----------------------------------
pcm_likelihood x = 1.79769e+308, log_lik = [-inf,-nan], normal_lpdf_val = -inf, retval = -nan
pcm_likelihood x = -1.79769e+308, log_lik = [0,-inf], normal_lpdf_val = -inf, retval = 0
[...]
pcm_likelihood x = 1.79769e+308, log_lik = [-nan,-nan], normal_lpdf_val = -inf, retval = -nan

Chain 1: Exception: gradient_of_f: The gradient of f is nan for parameter 0  (in 'modele66b52768eeb_pcm_bs_marg_nobeta_discourse' at line 81)


It looks like that in the first 1.79769e+308 we get [-inf,-nan] and in the second we get [-nan,-nan]. For x=\infty, the likelihood contribution of the first k should be 0 and the second 1, so I understand whrere the -inf in [-inf, -nan] comes from, but not the -nan (it should be \log(1)=0).

Of course, this is a great suggestion! This works, even in the big model (I=1000, K=10, J=4 + mixture). However, it is insanely slow, even slower than the adaptive Gauss Hermite quadrature:

AGH:

Chain 1: Gradient evaluation took 1.20227 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 12022.7 seconds.


integrate_1d:

Chain 1: Gradient evaluation took 86.1009 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 861009 seconds.


Do you have any suggestion how this could be sped up?

How many parameters passed to the integrator?

Can you write out the Math for the integral you’re trying to evaluate? Not clear to me looking at the code but maybe @nhuurre can see something in the Eqs.

2 Likes

I noticed that the term

log_softmax(
append_row(0, cumulative_sum(xmbetas[, k]))
)


is needlessly computed for every i\in I. So I basically copy-pasted the whole sinh-sinh integrator from the boost GitHub, dropping everything unnecessary for my specific function and further employing a cutoff criterion by @nhuurre’s suggestion.

I managed to bring it down from

1000 transitions using 10 leapfrog steps per transition would take 924607 seconds.


to

1000 transitions using 10 leapfrog steps per transition would take 1982 seconds.


so a couple hundred times faster, even though still very slow (without integrating out the person parameters, this model takes about 2 hours for 2000 iter/chain instead of a couple of days now).

Since I have a lot of cores at my disposal, I want to try to implement reduce_sum now, maybe that brings along another big acceleration.

Thanks again both of you for your help, I have one more question that came up which I will investigate empirically, but I would be curious of your opinion, too:
When testing my accelerated quadrature, I picked out some specific log-lik terms and every time, the difference to the integrate_1d function turned out to be zero. However, there must be some slight differences since the root mean square error of the whole 1000-vector of differences between the likelihood terms of my integrator and the integrate_1d was always around the 5e-14 range. On the other hand, benchmarks say I can double the speed again with only 5 levels of refinement, but then the RMSE of difference between my integrator and the integrate_1d is only about 1.5e-6.

Since the first is slightly above machine epsilon I was wondering: What level of precision still matters for stan? Probably it heavily depends on the model and the data, but do you have an intuition where you say “okay after x significant digits it doesn’t really matter anymore”?

It is good these this had high payoff, but can you talk a bit about how it was difficult/not possible to do it with integrate_1d? Ideally integrate_1d would just work and no need to send people off to implement their own sinh-sinh quadratures :P.

If you have y \sim N(f, \sigma) and f is the thing you’re computing, my intuition is you’re probably okay as long as the scale of the error in f is much less than either the difference between f and y or \sigma you’ll be okay.

The other thing you can do is compute a log density with the high precision calculation, compute it with the low precision approximation, and use the psis stuff in the loo package to check if an importance sampling correction is possible. The idea is that if p_{high}/p_{low} is close to 1.0 then the importance sampling will work, so you can just sampling using the low precision model, compute the log density of the high precision model in generated quantities, check the ratios, and then resample according to the importance sampling weights. Programmatically it gets a little awkward, but it’s a way to check if you’re off/do a correction if possible that should be pretty hands-off.

Ask if you want more details on the second thing @jtimonen

1 Like

After using the suggestions above, I managed to make it work and in the end it was just a matter of speed. With integrate_1d I broke off the estimation with 2 chains and 2000 iter/chain after a week or so, because it never reached 200/2000 iterations. The implementation I have now finishes in 1.5h. I would be curious how fast it will go when I reduce_sum the model with integrate_1d and what will happen in comparison to my implementation with reduce_sum.

Just some empirical results here: I estimated my model with my implementation with varying levels of refinement and I noticed that as soon as I have >=5 levels, my posterior means don’t change any more. This is in contrast to the maximum of 7 levels that the boost integrator uses.

Interestingly, using 5 levels gets me errors (so |I_1 - I_0|/||f_0|| as implemented in the boost function, I know we talked about this somewhere else, I haven’t managed to contact the boost guys with my question yet) as big as \approx 0.006, which corresponds to setting the tolerance to 0.006 and using 7 levels in integrate_1d (edit: of course it does).

Lastly, thank you for the suggestion, this seems to be really interesting and I might come back to it in the future when I have more time to get into the details of this approach.

Well if you can get a suitable solution and it is this much faster than integrate_1d it makes me wonder if we’re doing something wrong in integrate_1d.

I don’t think there is anything wrong with the implementation in Stan, it’s just that I have a likelihood with 1000 factors each consisting of two integrals in a mixture, which means I have to compute 2000 integrals per iteration. At the core, the likelihood is a categorical_lpmf and what I did is just compute the probabilities at every abscissa value of the integrator beforehand and then use the loops to only pick the correct probability based on the data out of the pre-computed probabilities.

Also, from the printing I did above, I know that the integrator works as intended. The only thing I don’t understand is the number of abscissas the boost quadrature takes at every refinement level, because I don’t understand the function get_abscissa_row() in the boost code, but I don’t speak enough C++ to figure this out.

Maybe there is a way to do this way faster with the on-board tools but I didn’t see it. The code is rather long, but go ahead if you want to take a look. This is an old version, I improved it some more since then, but the effect is in the millisecond range now.

data {
// Counts
int<lower=1> I; // No. persons
int<lower=1> K; // No. items
int<lower=1> J; // No. parameters per item
int<lower=1> C; // No. Classes

// Data
int<lower=1, upper=J+1> y[I, K];  // Data as Persons by Items Matrix

// Priors
// Mixture Parameters
int<lower=1> lambda_prior;

// Scale Parameter
real scale_prior_mu;
real<lower=0> scale_prior_sigma;

// Integration parameters
int N_levels;
int max_N_nodes;
vector[max_N_nodes] abscissas[N_levels];
vector[max_N_nodes] weights[N_levels];
int node_indices[N_levels];
real tol;
}

transformed data {
int<lower=J> Jp1 = J + 1;
row_vector[K] zerovec_K = rep_row_vector(0, K);
row_vector[C] zerovec_C = rep_row_vector(0, C);
int max_node_indices = max(node_indices);
real log_half_pi = log(pi()) - log2();
}

parameters {
// Mixture Parameters
simplex[C] lambda;

// Item Parameters
matrix[J, K] beta[C];

// Person Parameters
ordered[C] mu_theta;
vector<lower=0>[C] sigma_theta;
}

transformed parameters {
vector[I] log_lik;
{
matrix[Jp1, K] item_probs_pos[C, N_levels, max_node_indices];
matrix[Jp1, K] item_probs_neg[C, N_levels, max_node_indices];
matrix[Jp1, K] item_probs_zero[C];
vector[max_node_indices] log_dnorm_vals_pos[C, N_levels];
vector[max_node_indices] log_dnorm_vals_neg[C, N_levels];
row_vector[C] log_dnorm_vals_zero;
int nonzero_node_indices_pos[C, N_levels, max_node_indices];
int nonzero_node_indices_neg[C, N_levels, max_node_indices];

for (c in 1:C) {
item_probs_zero[c] = append_row(zerovec_K, -beta[c]);

for (k in 1:K) {
item_probs_zero[c][, k] = log_softmax(cumulative_sum(item_probs_zero[c][, k]));
}

log_dnorm_vals_zero[c] = normal_lpdf(0 | mu_theta[c], sigma_theta[c]);

for (n in 1:N_levels) {
for (eval_index in 1:node_indices[n]) {
// print("eval_index = ", eval_index, ", dnormarg_pos = ", fabs((abscissas[n][eval_index] - mu_theta[c]) / sigma_theta[c]));
if (fabs((abscissas[n][eval_index] - mu_theta[c]) / sigma_theta[c]) < 100) {
item_probs_pos[c, n, eval_index] = append_row(zerovec_K, abscissas[n][eval_index] - beta[c]);

for (k in 1:K) {
item_probs_pos[c, n, eval_index][, k] = log_softmax(cumulative_sum(item_probs_pos[c, n, eval_index][, k]));
}

log_dnorm_vals_pos[c, n, eval_index] = normal_lpdf(abscissas[n][eval_index] | mu_theta[c], sigma_theta[c]);

nonzero_node_indices_pos[c, n, eval_index] = 1;
} else {
nonzero_node_indices_pos[c, n, eval_index] = 0;
}
}

for (eval_index in 1:node_indices[n]) {
// print("eval_index = ", eval_index, ", dnormarg_neg = ", fabs((abscissas[n][eval_index] + mu_theta[c]) / sigma_theta[c]));
if (fabs((abscissas[n][eval_index] + mu_theta[c]) / sigma_theta[c]) < 100) {
item_probs_neg[c, n, eval_index] = append_row(zerovec_K, -abscissas[n][eval_index] - beta[c]);

for (k in 1:K) {
item_probs_neg[c, n, eval_index][, k] = log_softmax(cumulative_sum(item_probs_neg[c, n, eval_index][, k]));
}

log_dnorm_vals_neg[c, n, eval_index] = normal_lpdf(-abscissas[n][eval_index] | mu_theta[c], sigma_theta[c]);

nonzero_node_indices_neg[c, n, eval_index] = 1;
} else {
nonzero_node_indices_neg[c, n, eval_index] = 0;
}
}
}
}

for (i in 1:I) {
row_vector[C] int_by_group = log_dnorm_vals_zero;
// real val;

for (c in 1:C) {
for (k in 1:K) {
int_by_group[c] += item_probs_zero[c][y[i, k], k];
}

// val = exp(int_by_group[c]);
// print(i, ",", c, ",", 0, ",", -mu_theta[c] / sigma_theta[c], ",", log_dnorm_vals_zero[c], ",", val);

int_by_group[c] = exp(int_by_group[c] + log_half_pi);

for (n in 1:2) {
row_vector[node_indices[n]] this_f_eval_pos;
row_vector[node_indices[n]] this_f_eval_neg;

for (eval_index in 1:node_indices[n]) {
if (nonzero_node_indices_pos[c, n, eval_index]) {
real this_log_f_evals_pos = log_dnorm_vals_pos[c, n, eval_index];

for (k in 1:K) {
this_log_f_evals_pos += item_probs_pos[c, n, eval_index][y[i, k], k];
}

// print(i, ",", c, ",",  abscissas[n][eval_index], ",", ( abscissas[n][eval_index] - mu_theta[c]) / sigma_theta[c], ",", log_dnorm_vals_pos[c, n, eval_index], ",", exp(this_log_f_evals_pos));

this_f_eval_pos[eval_index] = exp(this_log_f_evals_pos);
} else {
this_f_eval_pos[eval_index] = 0;
}
}

for (eval_index in 1:node_indices[n]) {
if (nonzero_node_indices_neg[c, n, eval_index]) {
real this_log_f_evals_neg = log_dnorm_vals_neg[c, n, eval_index];

for (k in 1:K) {
this_log_f_evals_neg += item_probs_neg[c, n, eval_index][y[i, k], k];
}

// print(i, ",", c, ",", -abscissas[n][eval_index], ",", (-abscissas[n][eval_index] - mu_theta[c]) / sigma_theta[c], ",", log_dnorm_vals_neg[c, n, eval_index], ",", exp(this_log_f_evals_neg));

this_f_eval_neg[eval_index] = exp(this_log_f_evals_neg);
} else {
this_f_eval_neg[eval_index] = 0;
}
}

int_by_group[c] += (this_f_eval_pos + this_f_eval_neg) * weights[n][1:node_indices[n]];
}

int_by_group[c] /= 2;

for (n in 3:N_levels) {
row_vector[node_indices[n]] this_f_eval_pos;
row_vector[node_indices[n]] this_f_eval_neg;

real I0 = int_by_group[c];
int_by_group[c] /= 2;

for (eval_index in 1:node_indices[n]) {
if (nonzero_node_indices_pos[c, n, eval_index]) {
real this_log_f_evals_pos = log_dnorm_vals_pos[c, n, eval_index];

for (k in 1:K) {
this_log_f_evals_pos += item_probs_pos[c, n, eval_index][y[i, k], k];
}

// print(i, ",", c, ",",  abscissas[n][eval_index], ",", ( abscissas[n][eval_index] - mu_theta[c]) / sigma_theta[c], ",", log_dnorm_vals_pos[c, n, eval_index], ",", exp(this_log_f_evals_pos));

this_f_eval_pos[eval_index] = exp(this_log_f_evals_pos);
} else {
this_f_eval_pos[eval_index] = 0;
}
}

for (eval_index in 1:node_indices[n]) {
if (nonzero_node_indices_neg[c, n, eval_index]) {
real this_log_f_evals_neg = log_dnorm_vals_neg[c, n, eval_index];

for (k in 1:K) {
this_log_f_evals_neg += item_probs_neg[c, n, eval_index][y[i, k], k];
}

// print(i, ",", c, ",", -abscissas[n][eval_index], ",", (-abscissas[n][eval_index] - mu_theta[c]) / sigma_theta[c], ",", log_dnorm_vals_neg[c, n, eval_index], ",", exp(this_log_f_evals_neg));

this_f_eval_neg[eval_index] = exp(this_log_f_evals_neg);
} else {
this_f_eval_neg[eval_index] = 0;
}
}

int_by_group[c] += (this_f_eval_pos + this_f_eval_neg) * weights[n][1:node_indices[n]] / 2^(n - 1);

if (fabs(I0 - int_by_group[c]) < int_by_group[c] * tol) {
break;
}
}
}

log_lik[i] = log(int_by_group * lambda); //sum(int_by_group); //
}
}
}

model {
// Priors
for (c in 1:C) {
real beta_sum = 0;
target += std_normal_lpdf(to_vector(beta[c]));
beta_sum += sum(beta[c]);

target += normal_lpdf(beta_sum | 0, 0.01 * J * K);
}

target += dirichlet_lpdf(lambda | rep_vector(lambda_prior, C));
target += std_normal_lpdf(mu_theta);
target += lognormal_lpdf(sigma_theta | scale_prior_mu, scale_prior_sigma);

target += sum(log_lik);
}

2 Likes