Truncated Gamma distribution parameter estimation

Question

I need to estimate the parameter of a Gamma distribution from the data, but I only have samples from 0.05 to 3 (most of the samples are concentrated here).

I tried MLE, but due to the truncation it is not quite accurate.

I've been doing some research but I haven't found anything that really works well.

Answer 1

Maximum Likelihood Estimation (MLE) works remarkably well, even for fairly small datasets. For truncation at an upper limit of $U$ and lower limit of $L$, simply divide the Gamma likelihood of any data value $x$ by the total probability of the interval $[L,U]$ to obtain the likelihood for the truncated distribution.

Here are examples of datasets (shown as histograms), the underlying distributions that generated them (colored lines), and their MLEs (black lines) for a range of shape parameters "alpha" and scale parameters "sigma" that might be encountered when truncating to the interval $[0.05, 3]$. The agreement between the fits and distributions is excellent.

Details appear in the following R code.

#
# Log likelihood.
# `theta` is shape, log scale.
#
log.Lambda <- function(theta, x, limits) {
  #
  # Extract the parameters from the arguments.
  #
  alpha <- theta[1]
  log.sigma <- theta[2]
  sigma <- exp(log.sigma)
  n <- length(x)
  #
  # Compute Gamma probabilities for the truncation limits.  Keep them as logs
  # for numerical accuracy and to avoid overflow.
  #
  p.lower <- pgamma(limits[1], alpha, scale=sigma, log.p=TRUE)
  p.upper <- pgamma(limits[2], alpha, scale=sigma, log.p=TRUE)
  #
  # Compute the variable and constant portions of the log likelihood.
  #
  l <- (alpha-1) * log(x) - x/sigma
  const <- alpha * log.sigma + lgamma(alpha) + p.lower + log(exp(p.upper - p.lower) - 1)
  #
  # Return the negative log likelihood.
  #
  return(-(sum(l) - n*const))
}
#
# Truncated Gamma distribution (for plotting).
#
dgammatrunc <- function(x, shape, scale, lower, upper) {
  dgamma(x, shape, scale=scale) / diff(pgamma(c(lower,upper), shape, scale=scale))
}
#
# Test.
#
library(ggplot2)
library(data.table)
upper <- 3
lower <- 0.05
n <- 32      # Sample size
set.seed(17)
#
# Test for a range of shapes and scales.
#
parameters <- expand.grid(alpha=c(1/2, 2, 5), sigma=c(1/2, 6))
x0 <- seq(lower, upper, length.out=101) # Prediction points, for plotting
X <- apply(parameters, 1, function(theta) {
  alpha <- theta[1]
  sigma <- theta[2]
  #
  # Generate data.
  #
  q <- runif(n, pgamma(lower, alpha, scale=sigma), pgamma(upper, alpha, scale=sigma))
  x <- qgamma(q, alpha, scale=sigma)
  # hist(x, freq=FALSE)
  #
  # ML fitting.
  #
  theta.0 <- c(mean(x), 0)
  fit <- nlm(log.Lambda, p=theta.0, x=x, limits=c(lower, upper))
  beta.hat <- fit$estimate
  alpha.hat <- beta.hat[1]
  sigma.hat <- exp(beta.hat[2])
  #
  # Return the data and fits in a form convenient for plotting.
  #
  list(data.table(x=x, alpha=alpha, sigma=sigma),
       data.table(x=x0, alpha=alpha, sigma=sigma,
                  y=dgammatrunc(x0, alpha, sigma, lower, upper),
                  y.hat=dgammatrunc(x0, alpha.hat, sigma.hat, lower, upper))
  )
})
Y <- rbindlist(lapply(X, function(x) x[[2]])) # Data for the graphs
X <- rbindlist(lapply(X, function(x) x[[1]])) # The samples themselves
#
# Plot the results.
#
binwidth <- (upper - lower)/ceiling(n^(0.6))
ggplot(X, aes(x)) +
  geom_histogram(binwidth=binwidth, aes(fill=ordered(alpha)), alpha=1/2, 
                 color="Black", show.legend=FALSE) + 
  geom_path(aes(x, n*y*binwidth, color=ordered(alpha)), size=2, data=Y,
            show.legend=FALSE) +
  geom_path(aes(x, n*y.hat*binwidth), size=1.5, data=Y, show.legend=FALSE) +
  facet_grid(sigma ~ alpha, scales="free_y", labeller = label_both) + 
  ggtitle(paste("MLEs for Samples of Size", n),
          "Fitted Distributions Shown in Black")

Answer 2

First it is more convenient to truncate gamma distribution, utilising the direct method of truncation

G(x) = (F(x,a,b) - F( zmin, a,b))/ (F(zmax,a,b) - F( zmin, a,b)). F are the gamma cumulative function and not (x-zmin).((b-zmin) to truncate in zmin it is te left bound ary zmax the right boundary So, the transformation to number etc may be in a linear way

the truncated gamma distribution keeps its form, when we perform g(x) x^k. We could obtain another gamma function; and this is easier to manage g(x)is the pdf gamma distribution (truncated) Besides, with the truncated gamma we avoid the long tails that spoil the transformations

example the k Moment is b^k * GAMMA(a+k)/Gamma(a)/( F(a,zmax/b)-F(z, zmin/b) it is very easy to manage