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