How to fit a generalized logistic function? I'm trying to fit models that predict probabilities: $h(X,B) \rightarrow (0,1)$
It struck me that for for a great many cases logistic regression seems like a bad link function, I expect an S-curve, but I also expect bounds on the predicted probability to be much larger than 0, and much less than 1.
For example, let's say I want to predict the likelihood of a first bankruptcy within 10 years, and I have some feature $x$ (e.g. number of years a loan in is default). Here I might expect, that for all x,  $p>>0$ because many people go bankrupt for reasons unrelated to loan defaults, and that $p<<1$ because beyond some point more years in loan default doesn't increase the risk of new bankruptcy, because the number of years a loan in is default doesn't completely determine probability, you'd expect $p<<1$.
Fitting such a probability function with logistic regression leads to a 
very poor fit:

The target function above is a (special case) of "generalized logistic function". In this case:
$$
{prob} = p_{min} + (p_{max}-p_{min})*logistic(X,B)
$$
with logistic inputs $X$ and coefficients $B$.
Are there good ways to optimize $p_{min}$, $p_{max}$ and $B$  for this type of function such that a measure of regression error is minimized (log error, squared error, ...), or a similar 'S'-shaped model that's likely to provide a better fit than a logistic curve?
 A: Given the binary response $y_i$ and the covariate $x_i$, $i=1,2,\dots,n$, the likelihood for your model is
$$
L(\beta_0,\beta_1,p_\text{min},p_\text{max})=\prod_{i=1}^n p_i^{y_i}(1-p_i)^{1-y_i}
$$
where each
$$
p_i=p_\text{min} + (p_\text{max} - p_\text{min})\frac1{1+\exp(-(\beta_0 + \beta_1 x_i)}.
$$
Just write a function computing the log of this an apply some general purpose optimization algorithm to maximise this numerically with respect to the four parameters.  For example, in R do:
# the log likelihood
loglik <- function(par,y,x) {
  beta0 <- par[1]
  beta1 <- par[2]
  pmin <- par[3]
  pmax <- par[4]
  p <- pmin + (pmax - pmin)*plogis(beta0 + beta1*x)
  sum(dbinom(y, size=1, prob=p, log=TRUE))
}
# simulated data
x <- seq(-10,10,len=1000)
y <- rbinom(n=length(x),size=1,prob=.2 + .6*plogis(.5*x))
# fit the model
optim(c(0, 0.5 ,.1, .9), loglik, control=list(fnscale=-1), y=y,x=x, lower=c(-Inf,-Inf,0,0),upper=c(Inf,Inf,1,1))

Note that to test for evidence of a lower plateau at $p_\text{min}$ in your data, your $H_0:p_\text{min}=0$ is at the boundary of the parameter space and the approximate/asymptotic distribution of $2(\log L(\hat\theta_1)-\log L(\hat\theta_0))$ is going to be a mixture of chi-square distributions with 1 and 0 degrees of freedom, see Self, S. G. & Liang, K. Asymptotic properties of maximum likelihood estimators and likelihood ratio tests under nonstandard conditions J. Amer. Statist. Assoc., 1987, 82, 605-610.
In the simpler case that there is only one plateau (so $p_\text{max}=1$ or $p_\text{min}=0$) the model is equivalent to a zero-inflated binary regression model that can be fitted with e.g. the glmmTMB R-package.
