I am quite excited to start my journey into bayesian PKPD!
I’m trying to reproduce the Friberg–Karlsson semi-mechanistic neutropenia model in Stan using Torsten’s pmx_solve_group_rk45 as documented in this wonderful tutorial. The model compiles fine, but I am encountering fatal issues while sampling as the MCMC consistently stalls due to ODE solver errors. Specifically, I see repeated warnings like:
parameters[X] is NaN/Inf, but must be finite!
Max number of iterations exceeded (1000)
This happens even when running 1 chain with just 1–3 subjects and reduced data. Eventually, the sampler gets stuck in the same place, rejecting proposals and never progressing past early warmup.
Things I’ve tried:
Custom init() function with conservative, plausible values for all parameters (e.g., CLHat = 10, kaHat = 2.0, omega = rep(0.2, ...)), avoiding near-zero or extreme values.
Small subset of data (1–3 subjects, 6–10 timepoints) to isolate issues.
Loosening solver tolerances to 1e-4 and increasing max_num_steps to 2000.
Switching from pmx_solve_group_rk45 to the parallel version pmx_solve_group_rk45_ip.
Confirmed that weight, DV1, DV2, and other input values are clean and positive.
Observations:
Even with a valid init, proposals often produce NaN or Inf in the ODE parameter vector passed to the solver.
Sometimes, parameters[1] is NaN and parameters[3] is Inf, despite being initialized safely.
Reducing warmup iterations too much causes adaptation warnings (not enough warmup iterations to fit three adaptation stages), but increasing warmup doesn’t resolve the core solver failure.
At this point, I suspect something deeper might be going wrong with either how the parameters are transformed or passed into the solver, or that the ODE system is sensitive to small perturbations early in warmup.
I’d appreciate any help reviewing model structure, parameterization strategies, or debugging tips specific to ODE-based models with stiff or constrained dynamics in Torsten. I am attaching the scripts I use to run this toy analysis that includes all the 15 subjects that have been originally used in the analysis.
PS. For visibility to the authors of the tutorial just wanted to tag them specifically — @charlesm93 , @yizhang , and @billg
Could the problem here be that you have circ clamped to machine_precision()? If x[8] + circ0 is negative, circ will be machine_precision(), and (circ0/circ)^gamma is virtually guaranteed to throw an Inf for almost any circ0 > 1 ∧ gamma > 1 due to overflow. fmax() might also be introducing a discontinuity in the derivative if x[8] + circ0 is negative as well? Perhaps this sum is consistently hovering close to 0 in the ODE solver?
Hi, I was able to run your exact code with no problems, so I can’t reproduce your error, but I’d say that reducing the amount of data won’t help - less information will make it harder, not easier. Things you can try:
Simulate your own data in Stan, and play with parameter values and solvers there. It’ll be faster than fitting, and you can test to see what kinds of parameter values cause problems to the solvers (and you can see if Corey’smachine_precision()hypothesis holds weight. When I write this model in Stan I don’t use the machine_precision()part, and it’s fine)
Tighten your priors to something absurd - in your linked model they’ve got something like parameter ~ lognormal(log(x), sd). Just makexthe true value the data was simulated from and sd = 0.001. See if you get the same errors.
Also, if you’re just starting out, the Friberg-Karlsson model probably isn’t the easiest one to play around with. Try one- and two-compartment PK models first. They’ve got single-subject and population models in that tutorial, and I’ve got a whole repo of PK and PK/PD models written in Stan, from simple one-compartment models up to the Friberg-Karlsson model. You can also see in there how to simulate using Stan and how you might be able to play with different solvers.
ODE inference can be quite tricky, so if you are getting started with this formulation or a particular kind of model it may be worth simplifying as much as possible and stepping up complexity to find where things break – that will also be useful when trying to include more interactions in large models, which tends to break at some point, even when the foundations are sound.
Those are mostly good approaches, but normally they shouldn’t all be necessary at once for a “simple enough” model. One thing I’d suggest, especially when getting started with a fresh implementation, is to hard-code some of the parameters and leave only one or two free to see how that behaves. Sometimes it’s just one problematic parameter that tends to go to a negative values when it should be strictly positive, or is too close to zero, and may be straightforward to enforce.
initparameter may be valid, but the neighboring region may not be well-behaved, and as HMC leapfrogs around the parameter space it may encounter an invalid values within an iteration (that is one example of a case where it would be useful to identify precisely which part of the implementation could be sensitive to errors).