I’m trying to use `integrated_1d`

to numerically integrate the following

for t, \gamma > 0 and \alpha, \beta \in \mathbb{R}. Note that this domain means I cannot use the `gamma_p`

function in Stan. For some parameter configurations, `integrate_1d`

does not like this integrand. The only reason I can think of is that in some reasonable configurations u^{\gamma - 1} might overflow (for \gamma < 1 and small u), so I rewrote I as

but this is also produces `integrate_1d`

errors, and occur with such frequency that I can’t get Stan to initialise.

```
suppressMessages(library(cmdstanr))
cmdstan_model <- cmdstan_model(
"some-path-to-make-reprex-happy/submodel-two-phi-step.stan"
)
#> Compiling Stan program...
stan_data <- list(
n_patients = 16L,
baseline_measurement = c(
-0.935695521120025, -2.12761603089136, -0.876325387003469, 0.931350429946891,
0.422731897638725, -1.30375210258178, 1.85158267797611, -0.289448169455647,
0.325903831278964, -1.2757407765375, 1.58794213975781, 1.59307842403321,
-0.0584097720952887, 1.80357260097634, -2.30428952138529, 0.655115279462301
),
log_crude_event_rate = 3.17644992456411
)
res <- cmdstan_model$sample(
data = stan_data,
chains = 1,
iter_warmup = 1,
iter_sampling = 0
)
#> Running MCMC with 1 chain...
#>
#> Chain 1 event_time: 3.15863 hazard_gamma: 3.00293 alpha: -1.14303 beta_one: -1.82362
#> Chain 1 surv_prob_term_two: 2578.11
#> Chain 1 event_time: 0.181123 hazard_gamma: 3.00293 alpha: -1.14303 beta_one: 1.38322
#> Chain 1 surv_prob_term_two: 0.00159067
#> Chain 1 event_time: 0.529865 hazard_gamma: 3.00293 alpha: -1.14303 beta_one: 1.01519
#> Chain 1 surv_prob_term_two: 0.0314054
#> Chain 1 event_time: 4.02702 hazard_gamma: 3.00293 alpha: -1.14303 beta_one: -0.887469
#> Chain 1 surv_prob_term_two: 599.159
#> Chain 1 event_time: 0.581211 hazard_gamma: 3.00293 alpha: -1.14303 beta_one: -0.830041
#> Chain 1 surv_prob_term_two: 0.0992707
#> Chain 1 event_time: 1.74709 hazard_gamma: 3.00293 alpha: -1.14303 beta_one: 0.952865
#> Chain 1 surv_prob_term_two: 0.459929
#> Chain 1 event_time: 0.137359 hazard_gamma: 3.00293 alpha: -1.14303 beta_one: -1.10788
#> Chain 1 surv_prob_term_two: 0.000978182
#> Chain 1 event_time: 0.216454 hazard_gamma: 3.00293 alpha: -1.14303 beta_one: 0.00135613
#> Chain 1 surv_prob_term_two: 0.00336121
#> Chain 1 event_time: 4.41431 hazard_gamma: 3.00293 alpha: -1.14303 beta_one: -1.39034
#> Chain 1 surv_prob_term_two: 10356.2
#> Chain 1 event_time: 0.847361 hazard_gamma: 3.00293 alpha: -1.14303 beta_one: 0.840636
#> Chain 1 surv_prob_term_two: 0.111383
#> Chain 1 event_time: 1.12224 hazard_gamma: 3.00293 alpha: -1.14303 beta_one: 0.780646
#> Chain 1 surv_prob_term_two: 0.226627
#> Chain 1 event_time: 0.628329 hazard_gamma: 3.00293 alpha: -1.14303 beta_one: -1.17817
#> Chain 1 surv_prob_term_two: 0.157637
#> Chain 1 event_time: 0.516754 hazard_gamma: 3.00293 alpha: -1.14303 beta_one: -1.66055
#> Chain 1 surv_prob_term_two: 0.0973701
#> Chain 1 event_time: 4.18989 hazard_gamma: 3.00293 alpha: -1.14303 beta_one: 0.899539
#> Chain 1 surv_prob_term_two: 1.48168
#> Chain 1 event_time: 0.798516 hazard_gamma: 3.00293 alpha: -1.14303 beta_one: -1.63808
#> Chain 1 surv_prob_term_two: 0.540543
#> Chain 1 event_time: 0.204784 hazard_gamma: 3.00293 alpha: -1.14303 beta_one: -1.27114
#> Chain 1 surv_prob_term_two: 0.00356428
#> Chain 1 event_time: 3.15863 hazard_gamma: 3.00293 alpha: -1.14303 beta_one: -1.82362
#> Chain 1 surv_prob_term_two: 2578.11
#> Chain 1 event_time: 0.181123 hazard_gamma: 3.00293 alpha: -1.14303 beta_one: 1.38322
#> Chain 1 Exception: integrate: error estimate of integral 0.000158824 exceeds the given relative tolerance times norm of integral (in '/var/folders/9q/xrp79_nd2d7dhnp5j6nncymd6zq2rx/T/Rtmp12DnEu/model-116b06733f3da.stan', line 16, column 4 to column 16)
#> Chain 1 Exception: integrate: error estimate of integral 0.000158824 exceeds the given relative tolerance times norm of integral (in '/var/folders/9q/xrp79_nd2d7dhnp5j6nncymd6zq2rx/T/Rtmp12DnEu/model-116b06733f3da.stan', line 16, column 4 to column 16)
#> Chain 1 Exception: integrate: error estimate of integral 0.000158824 exceeds the given relative tolerance times norm of integral (in '/var/folders/9q/xrp79_nd2d7dhnp5j6nncymd6zq2rx/T/Rtmp12DnEu/model-116b06733f3da.stan', line 16, column 4 to column 16)
#> Chain 1 Exception: integrate: error estimate of integral 0.000158824 exceeds the given relative tolerance times norm of integral (in '/var/folders/9q/xrp79_nd2d7dhnp5j6nncymd6zq2rx/T/Rtmp12DnEu/model-116b06733f3da.stan', line 16, column 4 to column 16)
#> Warning: Chain 1 finished unexpectedly!
#> Warning: No chains finished successfully. Unable to retrieve the fit.
```

Here are some examples I’ve encountered. Initially I thought the asymptote at u = 0 was going to be the issue, but this also happens when the integrand is equal to zero at u = 0.

```
# some prefab examples:
# Example - decaying behaviour
event_time <- 0.333964
hazard_gamma <- 0.280632
alpha <- 1.71436
beta_one <- 1.641
# Chain 1 Rejecting initial value:
# Chain 1 Error evaluating the log probability at the initial value.
# Chain 1 Exception: integrate: error estimate of integral 0.0568744 exceeds the
# given relative tolerance times norm of integral
f_integrand <- function(x) {
x^(hazard_gamma - 1) * exp(alpha * beta_one * x)
}
f_alt <- function(x) {
exp((hazard_gamma - 1) * log(x) + alpha * beta_one * x)
}
integrate(f_integrand, lower = 0, upper = event_time)
#> 3.338441 with absolute error < 0.00011
integrate(f_alt, lower = 0, upper = event_time)
#> 3.338441 with absolute error < 0.00011
curve(f_integrand, from = 0, to = event_time)
curve(f_alt, from = 0, to = event_time, add = T, col = "blue")
```

```
# And another - increasing behaviour
event_time <- 0.141916
hazard_gamma <- 4.43969
alpha <- -1.20662
beta_one <- -0.81199
#Chain 1 Exception: integrate: error estimate of integral 7.81929e-08 exceeds
f_integrand <- function(x) {
x^(hazard_gamma - 1) * exp(alpha * beta_one * x)
}
f_alt <- function(x) {
exp((hazard_gamma - 1) * log(x) + alpha * beta_one * x)
}
integrate(f_integrand, lower = 0, upper = event_time)
#> 4.33822e-05 with absolute error < 1e-14
integrate(f_alt, lower = 0, upper = event_time)
#> 4.33822e-05 with absolute error < 1e-14
curve(f_integrand, from = 0, to = event_time)
curve(f_alt, from = 0, to = event_time, add = T, col = "blue")
```

^{Created on 2021-01-26 by the reprex package (v0.3.0)}

The Stan code is attached if you want to play around:

submodel-two-phi-step.stan (2.4 KB)