7
$\begingroup$

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.

$\endgroup$
  • 2
    $\begingroup$ Could you expand a little on what you mean by "not quite accurate"? What went wrong, exactly? BTW, if the data are truly truncated, then all the samples must lie between the given bounds. Are you sure you actually have truncation and not, say, censored or missing values? $\endgroup$ – whuber Oct 3 '14 at 19:32
  • $\begingroup$ Afterwards I am using this distribution to make a prediction and the results are underestimated, that's what I meant. And yes, the data is truly truncated. I just collected samples larger than 0.05 and eliminated those larger than 3 because I have samples from other distributions in that interval. So yes, all the samples lie in that interval I didn't express myself correctly. Thanks! $\endgroup$ – ms1437 Oct 3 '14 at 20:44
  • 1
    $\begingroup$ I'm afraid I am even more confused now, because you seem to be using the words "sample," "distribution," "prediction," and maybe even "estimated" in unconventional ways. To make your situation clear, consider including an explanation of exactly how you are sampling and what your overall aim is, so we can understand the context of this question. $\endgroup$ – whuber Oct 3 '14 at 21:35
  • $\begingroup$ I am trying to predict the number of samples that are going to be over a threshold in an unknown case. In order to do so, I have to fit a Gamma distribution to a case where we know the number of samples that are above 0.05 and below 3. This samples follow a gamma distribution. I just used MLE, without considering that my samples were only taken in that interval, and the prediction was always lower than the actual number of cases above the threshold. That's why I need to take into account that I am only having samples from a portion of the space. Thanks you! $\endgroup$ – ms1437 Oct 4 '14 at 0:19
  • 1
    $\begingroup$ Yes, MLE on truncated data performed as if the data were not truncated will be biased. You need your ML estimation to take account of the truncation by doing ML on the truncated distribution itself. $\endgroup$ – Glen_b Oct 4 '14 at 20:59
3
$\begingroup$

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.

Figure

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")
$\endgroup$
0
$\begingroup$

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

$\endgroup$
  • 1
    $\begingroup$ Please learn to use $\LaTeX$ notation and update your post accordingly! $\endgroup$ – kjetil b halvorsen Feb 14 '17 at 21:26
  • 1
    $\begingroup$ @kjetil it may be helpful to point to the MathJax tutorial and quick reference at meta.math.SE since it covers the exact subset of LaTeX we're using $\endgroup$ – Glen_b May 31 '17 at 11:26

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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