Non-numeric argument to pgamma(...) in MakeADFun

19 views
Skip to first unread message

pthomp...@gmail.com

unread,
Nov 12, 2025, 7:15:13 PMNov 12
to TMB Users
Hi there,

I'm trying to build a fairly simple model using RTMB and am running into trouble. It's probably just because I'm forgetting something simple and haven't written a model from scratch in a while, but I figured I'd ask here either way.

The model is designed for categorical data that represent "bins" of numeric values. You could imagine that these values represent coarse observations or measurements for a continuous process (i.e., a simple hierarchical model where the latent values are continuous and the observed values are categorical). My first approach in designing this model has been to assume the data come from a Gamma distribution, and to estimate those parameters based on the probability of a given data point falling within a certain bin from the Gamma distribution (i.e., using the CDF or pgamma() function).

When I wrote up the model, I got a negative log-likelihood function that seems to work - I can call it and plug it into optim() to get reasonable parameter estimates. However I'd like to leverage RTMB for this problem if I can, and when I try to make the model object with MakeADFun(), I get the following error: "Error in pgamma(bin_upper, par1, par2) : Non-numeric argument to mathematical function". I traced this down a bit and for what it's worth, it seems like the parameters (par1 and par2) are the problem, not the data (bin_upper).

Here is some code that generates some fake data and fits the model:

library(tmbstan)

# Latent, continuous values (in real life we do not have these). Note the true parameter values
actual_values = rgamma(100, 5, 0.05)
# Bins for which the real data are classified.
value_bins = c(0, 50, 100, 200, 300, 451, Inf)
# Categorize each data point into one of the bins. This is more like the real data I'd be fitting the model to.
binned_values = cut(actual_values, value_bins)
# Rearrange the data into a 2-column matrix where each row is a datapoint containing the lower and upper limit of its respective bin. This allows us to efficiently calculate pgamma() for the relevant start and end point.
bin_integers = as.integer(binned_depths)
bin_high_low = cbind(dive_bins[bin_integers], dive_bins[bin_integers + 1])

# Objective function
neg_log_like = function(pars) {
 
  par1 = exp(pars[1])
  par2 = exp(pars[2])
 
  bin_upper = as.numeric(bin_high_low[, 2])
  bin_lower = as.numeric(bin_high_low[, 1])
 
  bin_probs = pgamma(bin_upper, par1, par2) - pgamma(bin_lower, par1, par2)
 
  -sum(log(bin_probs))
 
}

# This fails
nll_seals_tmb = MakeADFun(neg_log_like, parameters = log(c(5, 0.05)))

# This works (with a good initial guess - with a bad one the function is equal to Inf)
model_fit = optim(log(c(5, 0.05)), neg_log_like)

Let me know if you're able to diagnose this. Thank you in advance!
Peter R. Thompson

Ben Bolker

unread,
Nov 12, 2025, 9:01:14 PMNov 12
to tmb-...@googlegroups.com
There may be a few typos in your code/non-reproducible ...
binned_depths and dive_bins are both undefined. Are these supposed to be
binned_values and value_bins ? (At a guess, could passing Inf to pgamma
be causing problems ?)
--
Dr. Benjamin Bolker
Professor, Mathematics & Statistics and Biology, McMaster University
Associate chair (graduate), Mathematics & Statistics
Director, School of Computational Science and Engineering
* E-mail is sent at my convenience; I don't expect replies outside of
working hours.

pthomp...@gmail.com

unread,
Nov 13, 2025, 3:00:35 PMNov 13
to TMB Users
Eek! So sorry about that, Ben. You're absolutely right about the typos - the code below should run.

Re: Inf passed to pgamma(), I'm not sure that's entirely it, although may be an issue - if you replace "par1" and "par2" in the pgamma() call with numeric values (e.g., 5 for par1 and 0.05 for par2), the function works (even with some of the Inf's still in there). But I may be misunderstanding some of the subtle complexity there!

library(RTMB)
set.seed(123)

# Latent, continuous values (in real life we do not have these). Note the true parameter values
actual_values = rgamma(100, 5, 0.05)
# Bins for which the real data are classified.
value_bins = c(0, 50, 100, 200, 300, 451, Inf)
# Categorize each data point into one of the bins. This is more like the real data I'd be fitting the model to.
binned_values = cut(actual_values, value_bins)
# Rearrange the data into a 2-column matrix where each row is a datapoint containing the lower and upper limit of its respective bin. This allows us to efficiently calculate pgamma() for the relevant start and end point.
bin_integers = as.integer(binned_values)
bin_high_low = cbind(value_bins[bin_integers], value_bins[bin_integers + 1])


# Objective function
neg_log_like = function(pars) {
 
  par1 = exp(pars[1])
  par2 = exp(pars[2])
 
  bin_upper = as.numeric(bin_high_low[, 2])
  bin_lower = as.numeric(bin_high_low[, 1])
 
  bin_probs = pgamma(bin_upper, par1, par2) - pgamma(bin_lower, par1, par2)
 
  -sum(log(bin_probs))
 
}

# This fails
nll_seals_tmb = MakeADFun(neg_log_like, parameters = log(c(5, 0.05)))

# This works (with a good initial guess - with a bad one the function is equal to Inf)
model_fit = optim(log(c(5, 0.05)), neg_log_like)

Ben Bolker

unread,
Nov 13, 2025, 9:09:59 PMNov 13
to tmb-...@googlegroups.com
I don't have an answer, but I think I've narrowed it down.

Here's a super-minimal example, based on debugging the objective
function and/or options(error=recover)

library(RTMB)
p1 <- advector(5.0)
p2 <- advector(0.5)

try(pgamma(1, p1, p2))
## Error in pgamma(1, p1, p2) :
## Non-numeric argument to mathematical function

pgamma(1, 5.0, 0.5)
## [1] 0.0001721156

dgamma(1, p1, p2)
## class='advector'
## [1] 0.0007897535

So: pgamma() on plain R numeric values works, pgamma() on advectors
fails, dgamma() on advectors works.

Different sets of methods are defined for the two functions, but I'm not
quite sure what to fix.

On the other hand, things work a little better (although not
perfectly) if we specify scale instead of rate:


library(RTMB)
set.seed(123)

# Latent, continuous values (in real life we do not have these). Note
the true parameter values
actual_values = rgamma(100, 5, 0.05)
# Bins for which the real data are classified.
value_bins = c(0, 50, 100, 200, 300, 451, Inf)
# Categorize each data point into one of the bins. This is more like the
real data I'd be fitting the model to.
binned_values = cut(actual_values, value_bins)
# Rearrange the data into a 2-column matrix where each row is a
datapoint containing the lower and upper limit of its respective bin.
This allows us to efficiently calculate pgamma() for the relevant start
and end point.
bin_integers = as.integer(binned_values)
bin_high_low = cbind(value_bins[bin_integers], value_bins[bin_integers + 1])

# Objective function
neg_log_like = function(pars) {

par1 = exp(pars[1])
par2 = exp(pars[2])

bin_upper = as.numeric(bin_high_low[, 2])
bin_lower = as.numeric(bin_high_low[, 1])

bin_probs = pgamma(bin_upper, par1, scale = 1/par2) -
pgamma(bin_lower, par1, scale = 1/par2)

-sum(log(bin_probs))

}

# This fails
nll_seals_tmb = MakeADFun(neg_log_like, parameters = log(c(5, 0.05)))

nll_seals_tmb$fn()
nll_seals_tmb$gr()
neg_log_like(log(c(5, 0.05)))

# This works (with a good initial guess - with a bad one the function is
equal to Inf)
model_fit = optim(log(c(5, 0.05)), neg_log_like)
model_fit2 = optim(log(c(5, 0.05)), fn = nll_seals_tmb$fn,
gr = nll_seals_tmb$gr, method = "BFGS")
model_fit$par
model_fit2$par

But:

nll_seals_tmb$gr()
outer mgc: NaN
[,1] [,2]
[1,] NaN -15.49354
Warning message:
In EvalADFunObject(ADFun, theta, order = order, hessiancols = cols, :
incpl_gamma (indef) integrate unreliable: x=0.000000 shape=5.000000
n=1.000000 ier=1

?? hope that helps a bit ...
> --
> To post to this group, send email to us...@tmb-project.org. Before
> posting, please check the wiki and issuetracker at https://github.com/
> kaskr/adcomp/ <https://github.com/kaskr/adcomp/>. Please try to create a
> simple repeatable example to go with your question (e.g issues 154, 134,
> 51). Use the issuetracker to report bugs.
> ---
> You received this message because you are subscribed to the Google
> Groups "TMB Users" group.
> To unsubscribe from this group and stop receiving emails from it, send
> an email to tmb-users+...@googlegroups.com <mailto:tmb-
> users+un...@googlegroups.com>.
> To view this discussion visit https://groups.google.com/d/msgid/tmb-
> users/f238ee71-e44e-4e0f-8b50-0d8444ffc1d8n%40googlegroups.com <https://
> groups.google.com/d/msgid/tmb-users/f238ee71-
> e44e-4e0f-8b50-0d8444ffc1d8n%40googlegroups.com?
> utm_medium=email&utm_source=footer>.

Mark Clements

unread,
Nov 14, 2025, 3:42:36 AMNov 14
to Ben Bolker, tmb-...@googlegroups.com
Based on https://github.com/kaskr/RTMB/blob/master/RTMB/R/distributions.R#L319, one should use the shape and scale arguments. The following works:

pgamma(1, shape=advector(5), scale=1/advector(0.5))

In my limited experience with RTMB, it is helpful to keep the implementation's code close at hand.

Kindly, Mark.

From: tmb-...@googlegroups.com <tmb-...@googlegroups.com> on behalf of Ben Bolker <bbo...@gmail.com>
Sent: 14 November 2025 03:09
To: tmb-...@googlegroups.com <tmb-...@googlegroups.com>
Subject: Re: [TMB users] Non-numeric argument to pgamma(...) in MakeADFun
 
> posting, please check the wiki and issuetracker at https://eur01.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgithub.com%2F&data=05%7C02%7Cmark.clements%40ki.se%7C7cdfa57713d5401c8bd908de2322ec19%7Cbff7eef1cf4b4f32be3da1dda043c05d%7C0%7C0%7C638986830081851544%7CUnknown%7CTWFpbGZsb3d8eyJFbXB0eU1hcGkiOnRydWUsIlYiOiIwLjAuMDAwMCIsIlAiOiJXaW4zMiIsIkFOIjoiTWFpbCIsIldUIjoyfQ%3D%3D%7C0%7C%7C%7C&sdata=cG0Lja9TMu0nwO%2Bzuecbjn6aQgMyuoS%2B7SamJsVP904%3D&reserved=0
> kaskr/adcomp/ <https://eur01.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgithub.com%2Fkaskr%2Fadcomp%2F&data=05%7C02%7Cmark.clements%40ki.se%7C7cdfa57713d5401c8bd908de2322ec19%7Cbff7eef1cf4b4f32be3da1dda043c05d%7C0%7C0%7C638986830081874402%7CUnknown%7CTWFpbGZsb3d8eyJFbXB0eU1hcGkiOnRydWUsIlYiOiIwLjAuMDAwMCIsIlAiOiJXaW4zMiIsIkFOIjoiTWFpbCIsIldUIjoyfQ%3D%3D%7C0%7C%7C%7C&sdata=fPCv38QwbmYRbmOJJiTa8a9aPNdRC%2BCpSX4Cm%2BL1MiA%3D&reserved=0>. Please try to create a
> simple repeatable example to go with your question (e.g issues 154, 134,
> 51). Use the issuetracker to report bugs.
> ---
> You received this message because you are subscribed to the Google
> Groups "TMB Users" group.
> To unsubscribe from this group and stop receiving emails from it, send
> an email to tmb-users+...@googlegroups.com <mailto:tmb-
> users+un...@googlegroups.com>.
> users/f238ee71-e44e-4e0f-8b50-0d8444ffc1d8n%40googlegroups.com <https://
> groups.google.com/d/msgid/tmb-users/f238ee71-
> e44e-4e0f-8b50-0d8444ffc1d8n%40googlegroups.com?
> utm_medium=email&utm_source=footer>.

--
Dr. Benjamin Bolker
Professor, Mathematics & Statistics and Biology, McMaster University
Associate chair (graduate), Mathematics & Statistics
Director, School of Computational Science and Engineering
* E-mail is sent at my convenience; I don't expect replies outside of
working hours.


---
You received this message because you are subscribed to the Google Groups "TMB Users" group.



När du skickar e-post till Karolinska Institutet (KI) innebär detta att KI kommer att behandla dina personuppgifter. Här finns information om hur KI behandlar personuppgifter. 


Sending email to Karolinska Institutet (KI) will result in KI processing your personal data. You can read more about KI’s processing of personal data here. 

Kasper Kristensen

unread,
Nov 14, 2025, 4:17:51 AMNov 14
to TMB Users
Thanks for isolating the problem. I've posted a github issue here

pthomp...@gmail.com

unread,
Nov 18, 2025, 6:31:50 PM (14 days ago) Nov 18
to TMB Users
Thanks all for your comments on this! It works perfectly with the pgamma() call changed to specify scale instead of rate.

Peter
Reply all
Reply to author
Forward
0 new messages