4
$\begingroup$

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?

$\endgroup$
5
  • $\begingroup$ Is there any reason why you are not using the survreg function from the survival package to do this? $\endgroup$
    – Emma Jean
    Nov 25, 2019 at 19:40
  • $\begingroup$ 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. $\endgroup$
    – AlexR
    Nov 25, 2019 at 19:41
  • 2
    $\begingroup$ Do you know the $\tau_i$? Are they (possibly) different for each $i$? $\endgroup$
    – jbowman
    Nov 25, 2019 at 19:48
  • $\begingroup$ @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 $\endgroup$
    – AlexR
    Nov 25, 2019 at 19:49
  • $\begingroup$ @jbowman Yes, I do (I'll try to clarify - note the logL contains $\tau_i$-terms) And yes, they are. $\endgroup$
    – AlexR
    Nov 25, 2019 at 19:50

1 Answer 1

5
$\begingroup$

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.

Also, you are conflating censoring with truncation. From Wikipedia,

Truncation is similar to but distinct from the concept of statistical censoring. A truncated sample can be thought of as being equivalent to an underlying sample with all values outside the bounds entirely omitted, with not even a count of those omitted being kept. With statistical censoring, a note would be recorded documenting which bound (upper or lower) had been exceeded and the value of that bound. With truncated sampling, no note is recorded.

Finally, I wouldn't call this a penalized likelihood but just a likelihood.

$\endgroup$
0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.