Creating a custom distribution with flexsurvreg I'm interested in fitting a parametric survival model but would like to explore the use of a Beta distribution for this purpose rather than a Weibull, Exponential, etc. model.
The Beta distribution is useful in that it can be parameterized with a min and a max parameter which means that, for my purposes, I can constrain my resulting survival curve to only take values within those limits.  My hope is that the Beta is flexible enough to fit a survival curve that is comparable to other, more common distributions.
According to the flexsurvreg reference manual on the use of custom distributions, I need only create an input list that has a few features: a name, names parameters, etc.  My understanding is that I also need to define a probability density function and a cumulative density function for my distribution.
A Beta density parameterized with $\theta_{1}$/$\theta_{2}$/min/max parameters can be written as:
$$
Pr(X=x)=\frac{1}{Beta(\theta_{1},\theta_{2})}\cdot \frac{(x-Min)^{\theta_{1}-1} \cdot(Max-x)^{\theta_{2}-1}}{(Max-Min)^{\theta_{1}+\theta_{2}-1}}
$$
The next thing flexsurveg( ) needs is a list with some metadata about my distribution (initializing values, etc).  I'm trying to get a simple example working before I start adding predictors, etc., but I can't get a basic model to fit without hearing about "non-finite" function values.  I currently suspect that, when $\theta_1$ and/or $\theta_{2}$ are possibly large values, the exponents are "blowing up" and causing this issue.
Here's a workable example:
#Simulate some right-censored data
set.seed(8765309)
myData <- data.frame(i = 1:100,
                 time = ceiling(runif(100, min = 1, max = 60)),
                 cens = sample(0:1, 100, replace = TRUE),
                 age = ceiling(rnorm(100, mean = 45, sd = 15)))

#Define the beta density
dbetamax <- Vectorize(function(x, theta1, theta2, log = FALSE){
  d_y <- ifelse(x <= 0, 0,
            ifelse(x >= 120, 0,
                   (1 / beta(theta1, theta2)) * ( (x - 0)^(theta1- 1) ) * ( (120 - x)^(theta2-1) ) * (1 / ((120 - 0)^(theta1+theta2-1)))))
  ifelse(log, log(d_y), d_y)
})

#Define the custom distribution
custom.betamax <- list(name = "betamax",
                   pars = c("theta1", "theta2"),
                   location = "theta1",
                   transforms = c(log, log),
                   inv.transforms = c(exp, exp),
                   inits = function(t){
                     c(1, 1)
                   })

#Pass to flexsurv
flexsurv::flexsurvreg(survival::Surv(myData$time, myData$cens) ~ age, data = myData, dist = custom.betamax)

Which results in the following:

Forming cumulative distribution function...
Forming integrated rmst function...
Forming integrated mean function...
Error in integrate(f = function (u, ...)  : non-finite function value
Error: $ operator is invalid for atomic vectors

Does anybody have much experience with custom distributions in flexsurvreg( ).  If so?  If so, do you know what the nature of this error is? And, if so, how does one avoid integration issues like these?
 A: I know this is not flexsurvreg, or even R, but if you are interested in Python, you can do this in Python's lifelines package¹ also using custom model creation. 
Instead of defining a PDF, we define a cumulative hazard (often easier to express and reason about vs a PDF). It's all wrapped into a class, subclassed from ParametricUnivariateFitter. 
from lifelines.fitters import ParametricUnivariateFitter

from autograd_gamma import betainc
from autograd import numpy as np

# create some fake data
n = 100
T_actual = 10 * np.random.random(n)**2
T_censor = 10 * np.random.random(n)**2
E = T_actual < T_censor
T_obs = np.minimum(T_actual, T_censor)

MAX = 120
MIN = 0


class BetaFitter(ParametricUnivariateFitter):
    _fitted_parameter_names = ['alpha_', 'beta_']
    _bounds = [(0, None), (0, None)]

    def _cumulative_density(self, params, times):
        alpha_, beta_ = params
        return betainc(alpha_, beta_, (times - MIN) / (MAX-MIN))

    def _cumulative_hazard(self, params, times):
        return -np.log(1-self._cumulative_density(params, times))

Inference is done:
bf = BetaFitter()
bf.fit(T_obs, E)
bf.print_summary()

Plotting:
bf.plot_survival_function()

More documentation on custom models here.

¹I'm the author of lifelines.
