# Fitting a gamma distribution to truncated data

I am faced with the following truncation problem:

$$X_i \sim \Gamma(\alpha, \beta) \\ \delta_i = \chi(X_i \le \tau_i)$$

I can observe only $$(X_i, \tau_i)$$ where $$\delta_i = 1$$ and I have no a-priori information about the total number of observations.

Trying to fit $$X_i$$ leads to a log-likelihood with an additional term

$$\log\mathcal L(\alpha, \beta ; x_i) = \sum_{i=1}^{N}\log f_X(x_i;\alpha, \beta) - \sum_{i=1}^N \log 1 - P_X(\tau_i; \alpha, \beta)$$

Where $$f_X$$ and $$P_X$$ are PDF resp. CDF of the gamma distribution. Note that the first term is "just" the regular log-likelihood for $$Y\sim\Gamma(\alpha,\beta)$$ and the second term relates to the truncated samples.

Now, for the exponential case ($$\alpha\equiv 1$$) I managed to efficiently estimate this likelihood in R using

set.seed(42)
r <- 2
Ntrue <- 1000
xtrue <- rexp(Ntrue, r)
ytrue <- runif(Ntrue)
tau  <- 2
obs <- xtrue + ytrue < tau
x <- xtrue[obs]
t <- tau - ytrue[obs]
rhat <- optimize(
function(r) {log(r) - r * mean(x) - mean(log(1 - exp(-r * t)))},
c(1e-4, 1/mean(x)),
maximum = TRUE
)$maximum # = 1.927  A straight-forward generalization of this approach using general gamma distributions bears no success, however: set.seed(42) r <- 2 Ntrue <- 1000 xtrue <- rexp(Ntrue, r) ytrue <- runif(Ntrue) tau <- 2 obs <- xtrue + ytrue < tau x <- xtrue[obs] t <- tau - ytrue[obs] llik <- function(shape, rate) { mean(dgamma(x, shape = shape, rate = rate, log = TRUE)) - mean(pgamma(t, shape = shape, rate = rate, lower.tail = FALSE, log.p = TRUE)) } rhat <- optim(c(shape = 1, rate = 1), fn = function(x){-llik(shape = x["shape"], rate = x["rate"])})$par
# Way off (2.03e50, 2.39e52) and warnings related to NaNs in {d,p}gamma.


Do you have any guidance on how to optimize the penalized likelihood?

• Is there any reason why you are not using the survreg function from the survival package to do this? Nov 25, 2019 at 19:40
• Pardon my ignorance, but does survreg deal with unknown sample sizes? I don't have censoring in the individual observations, only non-random completely missing observations. Nov 25, 2019 at 19:41
• Do you know the $\tau_i$? Are they (possibly) different for each $i$? Nov 25, 2019 at 19:48
• @EmmaJean: Looks like plain unpenalized ML to me: > 1/unique(predict(survival::survreg(survival::Surv(x)~1, dist = "exponential"))) [1] 2.355137 > 1/mean(x) [1] 2.355137 Nov 25, 2019 at 19:49
• @jbowman Yes, I do (I'll try to clarify - note the logL contains $\tau_i$-terms) And yes, they are. Nov 25, 2019 at 19:50

The contribution to your log-likelihood function due to the truncation should be $$\log P_X(\tau_i;\alpha,\beta)$$ not $$\log 1 - P_X(\tau_i;\alpha,\beta)$$. Thus, I think you just need to change lower.tail = FALSE to lower.tail = TRUE in your llik function.