Stan improper termination due to prior

Modeling
1

When I ran the stan program to fit a model, all chains terminated improperly, with some chains executing several iterations and other chains not iterating at all, and the fit object it returned had all draws of all parameters being zero.

For instance, in the following output, chain 1 had no iteration and chain 2 performed 3 iterations:

SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 1).
Chain 1: 
Chain 1: Gradient evaluation took 0.276904 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 2769.04 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 

SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 2).
Chain 2: 
Chain 2: Gradient evaluation took 0.237885 seconds
Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 2378.85 seconds.
Chain 2: Adjust your expectations accordingly!
Chain 2: 
Chain 2: 
Chain 2: Iteration:    1 / 2000 [  0%]  (Warmup)
Chain 2: Iteration:    2 / 2000 [  0%]  (Warmup)
Chain 2: Iteration:    3 / 2000 [  0%]  (Warmup)

In my model, I had 2 parameters b,\nu\in (0,+\infty). For stability I coded them as logarithm logb and lognu. In the original code which produced the above error, I set logb, lognu both with unit normal prior distribution. After the code returned error, I set b,\nu (i.e. exp(logb), exp(lognu)) to be distributed uniformly on (0.7, 1.5), then the program executed successfully and returned meaningful fit object.

Would there be any solutions to allow the program to run successfully under the prior logb, lognu~ normal(0,1) ?

Would it be possible to see the gradient and every leapfrog of RStan?
2

To diagnose this issue, it would be very helpful to see the complete Stan specification of the model that isn’t working as expected. Better still would be a reproducible example.

3

I have added a print statement at the last of the model block to print the log probability, i.e.

model {
  \\ something
  print("lp = ", target());
}

and the output is as follows. I hope this could help.

SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 1).
Chain 1: lp = -1335.55

Chain 1: lp = -1324.53

Chain 1: 
Chain 1: Gradient evaluation took 0.231407 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 2314.07 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: lp = -1324.53


SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 2).
Chain 2: lp = -442.203

Chain 2: lp = -431.175

Chain 2: 
Chain 2: Gradient evaluation took 0.246433 seconds
Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 2464.33 seconds.
Chain 2: Adjust your expectations accordingly!
Chain 2: 
Chain 2: 
Chain 2: lp = -431.175


SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 3).
Chain 3: lp = -379.296

Chain 3: lp = -368.269

Chain 3: 
Chain 3: Gradient evaluation took 0.217139 seconds
Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 2171.39 seconds.
Chain 3: Adjust your expectations accordingly!
Chain 3: 
Chain 3: 
Chain 3: lp = -368.269


SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 4).
Chain 4: lp = -222.755

Chain 4: lp = -211.728

Chain 4: 
Chain 4: Gradient evaluation took 0.236 seconds
Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 2360 seconds.
Chain 4: Adjust your expectations accordingly!
Chain 4: 
Chain 4: 
Chain 4: lp = -211.728

Chain 4: lp = -211.728

Chain 4: lp = -inf

Chain 4: lp = -211.728

Chain 4: lp = nan

Chain 4: lp = -211.728

Chain 4: lp = -211.728

Chain 4: lp = -146.643

Chain 4: lp = -211.728

Chain 4: lp = -156.312

Chain 4: Iteration:    1 / 2000 [  0%]  (Warmup)
Chain 4: lp = -211.728

Chain 4: lp = -149.854

Chain 4: lp = -154.063

Chain 4: lp = -124.302

Chain 4: Iteration:    2 / 2000 [  0%]  (Warmup)
Chain 4: lp = -154.063

Chain 4: Iteration:    3 / 2000 [  0%]  (Warmup)
Chain 4: lp = -154.063

Besides, if my sample has 100 observations, and I only select 15 of them as data input into the stan function (which means the log-likelihood is less likely to attain -M for dramatically large M>0), then each chain is more likely to perform some iterations (like chain 4 above) instead of performing no iterations at all (like other three chains above).

4

Sorry that my model is very complicated and I did not find a simpler example to reproduce the problem. I have tried to let the program print the log probability and the output is in the new reply. Would this help reveal some possible reason to the problem?

5

\mbox{log}(0.7) \approx -0.357 and \mbox{log}(1.5) \approx 0.406. A standard unit normal expands that range by almost a factor of 3 for the 95% limits. Without seeing your model, no one can say for sure, but a reasonable guess would be that the standard normal prior is allowing values that are too extreme for your model and the data.

6

The output that you pasted into your first two messages on this thread doesn’t indicate any obvious problems–it just looks like the sampling hasn’t finished yet. However, you also indicated that the code did return a fit object full of zeros. Then later in the same message, you indicated that the code errored.

Can you explain very specifically what happens when some chains fail to finish, and what code you run to produce a fit object containing all zeros?

7

I am sorry as the model is complicated, I could not post it entirely. I have let the program to give more detailed output as follows. I wish this would help.

Basically, my model has two parameters logb and lognu. They are 3\times 2 matrices and each entries are independently of standard normal distribution. And I have a likelihood function p(x| \log b, \log nu), coded as Modelpdf below. Here x is an ordered pair (n,t) where n\in\{1, 2, 3\} and t>0.

The data, parameters and model block of my code is as follows, where I have set intensive print statement:

data {
  int<lower=0> N; // N is number of samples.
  matrix[N,2] y;
}
parameters {
  matrix[3,2] logb;
  matrix[3,2] lognu;
  }
model {
  print("lpstart=", target());
  to_vector(lognu) ~ normal(0,1);
  to_vector(logb) ~ normal(0,1);
  
  print("lp = ", target());
  print("b = ", exp(logb) );
  print("nu = ", exp(lognu) );
  
  for (n in 1:N) {
      target += log(Modelpdf(to_int(floor(y[n,1]+0.1)), y[n,2], logb, lognu)) ;
      }
  print("lp = ", target());
  print("b = ", exp(logb) );
  print("nu = ", exp(lognu) );
}

I am sure that the Modelpdf is correct as I wrote it in Matlab also and compared the output of the function in Matlab versus the one in stan.

The code I used in rstan to produce the fit object is

test<-list(N=100,y=y)
fit1 <- stan(
  file = "model2.stan", 
  data = test, 
  algorithm="HMC",
  chains = 4,            
  warmup = 1000,          
  iter = 2000,           
  #test_grad = TRUE,
  cores = 1,              
  refresh = 1             
  )

and the output, given the above print statements in the model block, was

SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 1).
Chain 1: lpstart=0
lp = -20.5882
b = [[1.28101,1.11858],[0.342194,0.219518],[4.84288,1.28024]]
nu = [[0.138681,4.57043],[0.187903,6.15098],[0.747808,2.24988]]
lp = -744.816
b = [[1.28101,1.11858],[0.342194,0.219518],[4.84288,1.28024]]
nu = [[0.138681,4.57043],[0.187903,6.15098],[0.747808,2.24988]]

Chain 1: lpstart=0
lp = -9.56093
b = [[1.28101,1.11858],[0.342194,0.219518],[4.84288,1.28024]]
nu = [[0.138681,4.57043],[0.187903,6.15098],[0.747808,2.24988]]
lp = -733.788
b = [[1.28101,1.11858],[0.342194,0.219518],[4.84288,1.28024]]
nu = [[0.138681,4.57043],[0.187903,6.15098],[0.747808,2.24988]]

Chain 1: 
Chain 1: Gradient evaluation took 0.29386 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 2938.6 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: lpstart=0
lp = -9.56093
b = [[1.28101,1.11858],[0.342194,0.219518],[4.84288,1.28024]]
nu = [[0.138681,4.57043],[0.187903,6.15098],[0.747808,2.24988]]
lp = -733.788
b = [[1.28101,1.11858],[0.342194,0.219518],[4.84288,1.28024]]
nu = [[0.138681,4.57043],[0.187903,6.15098],[0.747808,2.24988]]

Chain 1: lpstart=0
lp = -101664
b = [[0.258684,1.40969e+35],[3.32381e+11,7.5972e+20],[1.58222,2.09183e-16]]
nu = [[0.663955,9.22383e-66],[0.802902,8.30835e-180],[0.660117,1.57487]]

Chain 1: lpstart=0
lp = -9.56093
b = [[1.28101,1.11858],[0.342194,0.219518],[4.84288,1.28024]]
nu = [[0.138681,4.57043],[0.187903,6.15098],[0.747808,2.24988]]
lp = -733.788
b = [[1.28101,1.11858],[0.342194,0.219518],[4.84288,1.28024]]
nu = [[0.138681,4.57043],[0.187903,6.15098],[0.747808,2.24988]]

Chain 1: lpstart=0
lp = -102123
b = [[0.0649057,1.23137e+35],[1.05606e+11,2.18032e+19],[2.40351,9.31028e-17]]
nu = [[0.0928038,3.89123e-66],[0.136592,2.49945e-180],[2.4847,4.22241]]

Chain 1: lpstart=0
lp = -9.56093
b = [[1.28101,1.11858],[0.342194,0.219518],[4.84288,1.28024]]
nu = [[0.138681,4.57043],[0.187903,6.15098],[0.747808,2.24988]]
lp = -733.788
b = [[1.28101,1.11858],[0.342194,0.219518],[4.84288,1.28024]]
nu = [[0.138681,4.57043],[0.187903,6.15098],[0.747808,2.24988]]

Chain 1: lpstart=0
lp = -6195.87
b = [[0.610048,1.07364e+09],[293.233,33501.6],[3.74232,0.000100333]]
nu = [[0.411694,1.79249e-16],[0.203059,4.90851e-45],[0.84552,8.5348]]

Chain 1: lpstart=0
lp = -9.56093
b = [[1.28101,1.11858],[0.342194,0.219518],[4.84288,1.28024]]
nu = [[0.138681,4.57043],[0.187903,6.15098],[0.747808,2.24988]]
lp = -733.788
b = [[1.28101,1.11858],[0.342194,0.219518],[4.84288,1.28024]]
nu = [[0.138681,4.57043],[0.187903,6.15098],[0.747808,2.24988]]

Chain 1: lpstart=0
lp = -347.995
b = [[0.88663,152.109],[2.62581,5.21633],[2.05429,0.0914965]]
nu = [[0.183155,0.000223793],[0.19592,3.19137e-11],[0.86744,2.0867]]

Chain 1: lpstart=0
lp = -9.56093
b = [[1.28101,1.11858],[0.342194,0.219518],[4.84288,1.28024]]
nu = [[0.138681,4.57043],[0.187903,6.15098],[0.747808,2.24988]]
lp = -733.788
b = [[1.28101,1.11858],[0.342194,0.219518],[4.84288,1.28024]]
nu = [[0.138681,4.57043],[0.187903,6.15098],[0.747808,2.24988]]

Chain 1: lpstart=0
lp = -17.0888
b = [[1.13499,3.59844],[0.56805,0.465291],[4.95263,0.656462]]
nu = [[0.172349,0.454657],[0.186865,0.00977451],[0.742156,2.50926]]
lp = -303
b = [[1.13499,3.59844],[0.56805,0.465291],[4.95263,0.656462]]
nu = [[0.172349,0.454657],[0.186865,0.00977451],[0.742156,2.50926]]

Chain 1: lpstart=0
lp = -9.56093
b = [[1.28101,1.11858],[0.342194,0.219518],[4.84288,1.28024]]
nu = [[0.138681,4.57043],[0.187903,6.15098],[0.747808,2.24988]]
lp = -733.788
b = [[1.28101,1.11858],[0.342194,0.219518],[4.84288,1.28024]]
nu = [[0.138681,4.57043],[0.187903,6.15098],[0.747808,2.24988]]

Chain 1: lpstart=0
lp = -6.96571
b = [[1.1569,1.51236],[0.356334,0.291846],[4.6475,1.20018]]
nu = [[0.131149,2.60687],[0.186482,1.19334],[0.67535,2.32727]]
lp = -181.823
b = [[1.1569,1.51236],[0.356334,0.291846],[4.6475,1.20018]]
nu = [[0.131149,2.60687],[0.186482,1.19334],[0.67535,2.32727]]

Chain 1: Iteration:    1 / 2000 [  0%]  (Warmup)
Chain 1: lpstart=0
lp = -9.56093
b = [[1.28101,1.11858],[0.342194,0.219518],[4.84288,1.28024]]
nu = [[0.138681,4.57043],[0.187903,6.15098],[0.747808,2.24988]]
lp = -733.788
b = [[1.28101,1.11858],[0.342194,0.219518],[4.84288,1.28024]]
nu = [[0.138681,4.57043],[0.187903,6.15098],[0.747808,2.24988]]

Chain 1: lpstart=0
lp = -6.62342
b = [[1.35014,1.58544],[0.355369,0.315736],[4.82317,1.16233]]
nu = [[0.143983,2.4156],[0.191915,1.22097],[0.778794,2.28856]]
lp = -176.541
b = [[1.35014,1.58544],[0.355369,0.315736],[4.82317,1.16233]]
nu = [[0.143983,2.4156],[0.191915,1.22097],[0.778794,2.28856]]

Chain 1: lpstart=0
lp = -6.98545
b = [[1.47061,2.37171],[0.486639,0.504675],[4.77374,0.78277]]
nu = [[0.150132,1.31862],[0.19451,0.207015],[0.811868,2.3221]]
lp = -172.839
b = [[1.47061,2.37171],[0.486639,0.504675],[4.77374,0.78277]]
nu = [[0.150132,1.31862],[0.19451,0.207015],[0.811868,2.3221]]

Chain 1: lpstart=0
lp = -10.7785
b = [[1.67228,2.40545],[0.843087,0.738026],[4.69247,0.550689]]
nu = [[0.157202,0.861084],[0.198997,0.0354028],[0.847158,1.90158]]
lp = -187.928
b = [[1.67228,2.40545],[0.843087,0.738026],[4.69247,0.550689]]
nu = [[0.157202,0.861084],[0.198997,0.0354028],[0.847158,1.90158]]

Chain 1: lpstart=0
lp = -17.9047
b = [[1.94346,1.53232],[1.49993,0.802787],[4.56981,0.410436]]
nu = [[0.165466,0.645804],[0.208855,0.00615597],[0.885124,1.28946]]
lp = -169.085
b = [[1.94346,1.53232],[1.49993,0.802787],[4.56981,0.410436]]
nu = [[0.165466,0.645804],[0.208855,0.00615597],[0.885124,1.28946]]

Chain 1: lpstart=0
lp = -28.7183
b = [[2.2922,0.630646],[2.18052,0.621359],[4.40839,0.314732]]
nu = [[0.175195,0.555789],[0.224462,0.00109276],[0.925936,0.803149]]
lp = -146.783
b = [[2.2922,0.630646],[2.18052,0.621359],[4.40839,0.314732]]
nu = [[0.175195,0.555789],[0.224462,0.00109276],[0.925936,0.803149]]

Chain 1: lpstart=0
lp = -43.3923
b = [[2.71318,0.291384],[3.08589,0.446546],[4.21095,0.257959]]
nu = [[0.186695,0.498598],[0.242963,0.00019923],[0.96985,0.481822]]
lp = -200.806
b = [[2.71318,0.291384],[3.08589,0.446546],[4.21095,0.257959]]
nu = [[0.186695,0.498598],[0.242963,0.00019923],[0.96985,0.481822]]

Chain 1: Iteration:    2 / 2000 [  0%]  (Warmup)
Chain 1: lpstart=0
lp = -43.3923
b = [[2.71318,0.291384],[3.08589,0.446546],[4.21095,0.257959]]
nu = [[0.186695,0.498598],[0.242963,0.00019923],[0.96985,0.481822]]
lp = -200.806
b = [[2.71318,0.291384],[3.08589,0.446546],[4.21095,0.257959]]
nu = [[0.186695,0.498598],[0.242963,0.00019923],[0.96985,0.481822]]

Chain 1: lpstart=0
lp = -2.66669e+07
b = [[4.38545e-24,inf],[2.38016e-44,inf],[2.4329e-152,inf]]
nu = [[2.72374e+66,4.58482e-67],[6.01771e+71,inf],[3.18612e+20,9.90321e-195]]


SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 2).
Chain 2: lpstart=0
lp = -15.3714
b = [[0.256505,0.268295],[1.0411,0.493499],[0.743356,0.619902]]
nu = [[0.309038,1.76606],[4.21477,0.772639],[1.94766,0.924812]]
lp = -452.584
b = [[0.256505,0.268295],[1.0411,0.493499],[0.743356,0.619902]]
nu = [[0.309038,1.76606],[4.21477,0.772639],[1.94766,0.924812]]

Chain 2: lpstart=0
lp = -4.34415
b = [[0.256505,0.268295],[1.0411,0.493499],[0.743356,0.619902]]
nu = [[0.309038,1.76606],[4.21477,0.772639],[1.94766,0.924812]]
lp = -441.557
b = [[0.256505,0.268295],[1.0411,0.493499],[0.743356,0.619902]]
nu = [[0.309038,1.76606],[4.21477,0.772639],[1.94766,0.924812]]

Chain 2: 
Chain 2: Gradient evaluation took 0.316785 seconds
Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 3167.85 seconds.
Chain 2: Adjust your expectations accordingly!
Chain 2: 
Chain 2: 
Chain 2: lpstart=0
lp = -4.34415
b = [[0.256505,0.268295],[1.0411,0.493499],[0.743356,0.619902]]
nu = [[0.309038,1.76606],[4.21477,0.772639],[1.94766,0.924812]]
lp = -441.557
b = [[0.256505,0.268295],[1.0411,0.493499],[0.743356,0.619902]]
nu = [[0.309038,1.76606],[4.21477,0.772639],[1.94766,0.924812]]

Chain 2: lpstart=0
lp = -20968.8
b = [[1.05277e+17,1.71171e+11],[2.31555e+44,2436.64],[1.10887e+29,2.96038e+12]]
nu = [[0.00210357,4.01181e-17],[2.57692e-59,0.00234713],[3.53367e-28,1.66718e-07]]


SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 3).
Chain 3: lpstart=0
lp = -21.123
b = [[0.22941,0.226122],[0.899819,6.68798],[2.61169,1.753]]
nu = [[1.32054,1.42598],[6.61398,4.62156],[6.84219,2.91001]]
lp = -1009.9
b = [[0.22941,0.226122],[0.899819,6.68798],[2.61169,1.753]]
nu = [[1.32054,1.42598],[6.61398,4.62156],[6.84219,2.91001]]

Chain 3: lpstart=0
lp = -10.0957
b = [[0.22941,0.226122],[0.899819,6.68798],[2.61169,1.753]]
nu = [[1.32054,1.42598],[6.61398,4.62156],[6.84219,2.91001]]
lp = -998.877
b = [[0.22941,0.226122],[0.899819,6.68798],[2.61169,1.753]]
nu = [[1.32054,1.42598],[6.61398,4.62156],[6.84219,2.91001]]

Chain 3: 
Chain 3: Gradient evaluation took 0.305883 seconds
Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 3058.83 seconds.
Chain 3: Adjust your expectations accordingly!
Chain 3: 
Chain 3: 
Chain 3: lpstart=0
lp = -10.0957
b = [[0.22941,0.226122],[0.899819,6.68798],[2.61169,1.753]]
nu = [[1.32054,1.42598],[6.61398,4.62156],[6.84219,2.91001]]
lp = -998.877
b = [[0.22941,0.226122],[0.899819,6.68798],[2.61169,1.753]]
nu = [[1.32054,1.42598],[6.61398,4.62156],[6.84219,2.91001]]

Chain 3: lpstart=0
lp = -414028
b = [[3.60783e+17,0.333854],[5.96048e+104,0.736455],[3.84746e+105,1.62398e+25]]
nu = [[0.000421326,0.0001302],[0,0.0761439],[1.2502e-160,7.40588e-25]]


SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 4).
Chain 4: lpstart=0
lp = -19.9926
b = [[1.25489,5.86404],[6.19436,0.831945],[0.284353,2.42538]]
nu = [[5.13431,1.48601],[3.02993,0.404332],[0.15197,2.15936]]
lp = -1453.3
b = [[1.25489,5.86404],[6.19436,0.831945],[0.284353,2.42538]]
nu = [[5.13431,1.48601],[3.02993,0.404332],[0.15197,2.15936]]

Chain 4: lpstart=0
lp = -8.96531
b = [[1.25489,5.86404],[6.19436,0.831945],[0.284353,2.42538]]
nu = [[5.13431,1.48601],[3.02993,0.404332],[0.15197,2.15936]]
lp = -1442.27
b = [[1.25489,5.86404],[6.19436,0.831945],[0.284353,2.42538]]
nu = [[5.13431,1.48601],[3.02993,0.404332],[0.15197,2.15936]]

Chain 4: 
Chain 4: Gradient evaluation took 0.240945 seconds
Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 2409.45 seconds.
Chain 4: Adjust your expectations accordingly!
Chain 4: 
Chain 4: 
Chain 4: lpstart=0
lp = -8.96531
b = [[1.25489,5.86404],[6.19436,0.831945],[0.284353,2.42538]]
nu = [[5.13431,1.48601],[3.02993,0.404332],[0.15197,2.15936]]
lp = -1442.27
b = [[1.25489,5.86404],[6.19436,0.831945],[0.284353,2.42538]]
nu = [[5.13431,1.48601],[3.02993,0.404332],[0.15197,2.15936]]

Chain 4: lpstart=0
lp = -1.08691e+06
b = [[1.19098e+63,3.07881],[0,8.71634e+07],[9.24704e+19,668397]]
nu = [[1.62845e-148,9.38263],[9.98345e+140,0.0108514],[0.00685156,0.000196538]]

Here the data y was generated by simulation with ground truth \log b= \log \begin{pmatrix} 1& 1.1 \\ 1.2 & 1.32\\ 1.3 & 1.43 \end{pmatrix} and \log nu (i,j)=0 (i.e. nu(i,j)= 1 for all i,j). Hence in the last output of each chain, the b and nu are really ridiculously large or small. This might be the reason why it terminated unexpectedly. But how to solve this?

8

The output you’ve shown does not look complete. It looks like either the model fitting was still in progress, or alternatively the model fitting halted unexpectedly, in which case I would have expected some kind of warning or error message. What part of the output makes you say that it terminated unexpectedly?

9

Yes this is exactly where my problem is strange. This is the screenshot of the Rstudio console. The execution just stopped unexpectedly with a > sign in the new command line meaning the current execution is over, and a fit object with all draws of all parameters being zero is returned.

image
image1310×190 14.7 KB

10

Have you tried running the same model via CmdStanR/CmdStanPy/CmdStan? This thread has a similar story with Rstudio and rstan, but it’s very old and doesn’t seem to have any resolution.

I assume that the function Modelpdf is pretty complicated. Is it a continuous function of the parameters? Are the values of logb or lognu ever involved in control flow?

11

Thanks. I have tried to run it with cmdStanR and the error message is as follows:

Running MCMC with 1 chain...

Chain 1 lpstart=0 
Chain 1 lp = -18.9524 
Chain 1 b = [[6.2492,0.646089],[4.46616,0.754293],[4.25773,0.44894]] 
Chain 1 nu = [[0.306784,0.291134],[0.186534,2.88347],[1.77894,1.2521]] 
Chain 1 lp = -185.353 
Chain 1 b = [[6.2492,0.646089],[4.46616,0.754293],[4.25773,0.44894]] 
Chain 1 nu = [[0.306784,0.291134],[0.186534,2.88347],[1.77894,1.2521]] 
Chain 1 lpstart=0 
Chain 1 lp = -7.92512 
Chain 1 b = [[6.2492,0.646089],[4.46616,0.754293],[4.25773,0.44894]] 
Chain 1 nu = [[0.306784,0.291134],[0.186534,2.88347],[1.77894,1.2521]] 
Chain 1 lp = -174.325 
Chain 1 b = [[6.2492,0.646089],[4.46616,0.754293],[4.25773,0.44894]] 
Chain 1 nu = [[0.306784,0.291134],[0.186534,2.88347],[1.77894,1.2521]] 
Chain 1 lpstart=0 
Chain 1 lp = -7.92512 
Chain 1 b = [[6.2492,0.646089],[4.46616,0.754293],[4.25773,0.44894]] 
Chain 1 nu = [[0.306784,0.291134],[0.186534,2.88347],[1.77894,1.2521]] 
Chain 1 lp = -174.325 
Chain 1 b = [[6.2492,0.646089],[4.46616,0.754293],[4.25773,0.44894]] 
Chain 1 nu = [[0.306784,0.291134],[0.186534,2.88347],[1.77894,1.2521]] 
Chain 1 lpstart=0 
Chain 1 lp = -3547.67 
Chain 1 b = [[9.47789,8033.76],[5.9068,6.0524e+19],[0.0413913,8778.36]] 
Chain 1 nu = [[0.760337,87.5313],[1.32905,1.30386e-30],[19.7712,0.000175897]] 
Chain 1 Exception initializing step size.
Chain 1 Exception: Exception: Exception: Error in function tanh_sinh<double>::integrate: The tanh_sinh quadrature found your function to be non-finite everywhere! Please check your function for singularities.
# (with information indicating the specific location like line 154, column 6 to column 79)

Warning: Chain 1 finished unexpectedly!

I believe my model is sufficiently regular wrt logb and lognu. Besides, the last output of logb and lognu was so extreme that program terminated with error. For the second last output, the loglikelihood of the sample evaluated for this logb and lognu was -166.4002, which is not too extreme.

12

This is great progress, since there is now an error message. It seems like you are doing some numerical integration and the call is failing because function is non-finite (everywhere?!). I can’t help much more without seeing more of the code, but at least you have a place to start debugging.

13

Yes, my model likelihood involves some integration like \int_0^y f(\log b,\log nu, u)du (upper limit being the sample and logb, lognu are parameters of the intergrand). Given such extreme logb and lognu, the intergral would be too close to 0+ and becomes negative infinity after taking logarithm. So my question now reduces to how did the sampler go to such extreme value of logb and lognu.

14

If you are expecting the integral to be close to 0, then the error message still doesn’t quite make sense to me, since the failure occurred during integration, not while taking the log of the result.

I can’t speculate on why the sampler got to the parameter values that triggered this crash. During warmup, my potentially flawed mental model is that the sampler explores regions of the latent space that are potentially very extreme, so if you want consistent runs, your functions need to be able to handle most inputs that are valid given the model constraints. You can include reject statements for particularly gnarly and rare cases.

It seems like you are computing probability density on the natural scale. It’s probably possible and potentially more numerically stable to do all computations on the log scale. This section of the user guide has been very helpful to me (and not just in the context of Stan), especially the use of the built-in log_sum_exp() function and its cousins.