3-parameter Gamma distribution stats not matching with Minitab

General
1

hi folks,

I am trying to calculate the 3-parameter Gamma distribution in R using several methods. I am able to produce results, but the estimated shape, scale, and threshold (location) values do not match the results provided by Minitab. Can anyone help me understand why this is happening and how to align the R results with Minitab?

here is my dataset

x <- c(1.2, 1.6, 1.7, 1.8, 1.9, 2.0, 2.2, 2.6, 3.0, 3.5, 4.0, 4.8, 5.6, 6.6, 7.6)

minitab Result
#-----------------------
#Shape: 2.37046
#scale: 1.19117
#Threshold: 0.51634

#-----------------------------------------------------
# Method-1                                          
# 3-parameter Gamma Negative Log-Likelihood         
#-----------------------------------------------------

negloglik_gamma3 <- function(params, data) {
  shape <- params[1]
  scale <- params[2]
  shift <- params[3]
  
  y <- data - shift
  
  if (shape <= 0 || scale <= 0 || any(y <= 0)) return(1e10)
  
  -sum(dgamma(y, shape = shape, scale = scale, log = TRUE))
}


x_moments <- x - min(x) + 0.01
init_shape <- mean(x_moments)^2 / var(x_moments)
init_scale <- var(x_moments) / mean(x_moments)
init_shift <- min(x) - 0.5  

init_params <- c(shape = init_shape, scale = init_scale, shift = init_shift)


lower_bounds <- c(1e-4, 1e-4, -Inf)
upper_bounds <- c(Inf, Inf, min(x) - 1e-6)


fit <- optim(par = init_params,
             fn = negloglik_gamma3,
             data = x,
             method = "L-BFGS-B",
             lower = lower_bounds,
             upper = upper_bounds,
             control = list(fnscale = 1))


shape_est <- fit$par[1]
scale_est <- fit$par[2]
shift_est <- fit$par[3]

cat("Estimated shape:", round(shape_est, 5), "\n")
cat("Estimated scale:", round(scale_est, 5), "\n")
cat("Estimated shift:", round(shift_est, 5), "\n")

#-----------------------------------------------------
#Result of method 1
#-----------------------------------------------------
# Estimated shape: 0.50503 
# Estimated scale: 4.23743 
# Estimated shift: 1.2  
#------------------------------------------------------
#  Method 2
#-----------------------------------------------------

library(FAdist); 
library(fitdistrplus)
fit<-fitdist(x,dgamma3,start=list(shape=1,scale=1,thres=0.5),
             upper=c(Inf,Inf,min(x)))

print(fit)
#-----------------------------------------------------------
# Result of method 2
#----------------------------------------------------------
#shape : 0.2922451 
#scale : 7.3236 
#thres : 1.2 
#---------------------------------------------------------
  
  
#-----------------------------------------------------------  
#---Method 3
#-----------------------------------------------------------


theta_vals <- seq(min(x) - 1, min(x) - 0.01, length.out = 100)

best_ll <- -Inf
best_params <- NULL

for (theta in theta_vals) {
  shifted_data <- x - theta
  if (any(shifted_data <= 0)) next
  
  sample_mean <- mean(shifted_data)
  sample_var <- var(shifted_data)
  
  shape_est <- sample_mean^2 / sample_var
  scale_est <- sample_var / sample_mean
  
  loglik <- sum(dgamma(shifted_data, shape = shape_est, scale = scale_est, log = TRUE))
  
  if (loglik > best_ll) {
    best_ll <- loglik
    best_params <- c(shape = shape_est, scale = scale_est, threshold = theta)
  }
}

cat("Best Estimates:\n")
cat("Shape (alpha):", best_params[1], "\n")
cat("Scale (beta):", best_params[2], "\n")
cat("Threshold (theta):", best_params[3], "\n")

#---------------------------------------------------------
# Result of method 3
#---------------------------------------------------------
# Shape (alpha): 1.199152 
# Scale (beta): 1.809612 
# Threshold (theta): 1.17 
#---------------------------------------------------------

thanks in advance.

2

Hi @Mahendar_Vanaparthi ,

if you look closely in the R results you will realize that all values from the different packages for the Threshold/Shift parameter are near the minimum value of your data.

While fitting parametric distributions which include a Threshold/Shift it can happen(which is happening in your example) that the loglikelihood is unbounded when the numerical optimization is getting close to the minimum value of your dataset.

The following was generated by Microsoft copilot since i am lazy to write this on my own:
BEGIN COPILOT ( Microsoft. (2025). Copilot (GPT-4) [Large language model] . https://copilot.microsoft.com)

The log-likelihood for a single observation ( x ) from a 3-parameter gamma distribution is:

\log \left( \frac{(x - \theta)^{\alpha - 1} e^{-(x - \theta)/\beta}}{\beta^\alpha \Gamma(\alpha)} \right)

This simplifies to:

(\alpha - 1) \log(x - \theta) - \frac{x - \theta}{\beta} - \alpha \log(\beta) - \log(\Gamma(\alpha))
  • α (alpha): Shape — controls the skewness and tail behavior.
  • β (beta): Scale — stretches or compresses the distribution.
  • θ (theta): Location (threshold) — shifts the distribution along the x-axis.

:warning: Behavior Near the Lower Bound

As \theta \to x^- , we have x - \theta \to 0^+ , and:

  • \log(x - \theta) \to -\infty
  • The term (\alpha - 1) \log(x - \theta) dominates the behavior

:white_check_mark: Cases:

  • If \alpha > 1 :

    • (\alpha - 1) \log(x - \theta) \to -\infty
    • The log-likelihood is bounded below

  • If \alpha < 1 :

    • (\alpha - 1) \log(x - \theta) \to +\infty
    • The log-likelihood becomes unbounded

  • If \alpha = 1 :

    • The log term disappears
    • The likelihood is bounded, but may still be steep

:white_check_mark: Conclusion:

The likelihood becomes unbounded as \theta \to \min(x_i) if and only if:

\boxed{\alpha < 1}

This is why penalization or bounding \theta is essential when estimating the location parameter in such cases.

END COPILOT
So basically all the algorithms have a problem converging due to this unbounded nature of the likelihood.

Using the Minitab's threshold and calling the base R optim function will give you (up to numerical error) the same answer

# your data
x <- c(1.2,
       1.6,
       1.7,
       1.8,
       1.9,
       2.0,
       2.2,
       2.6,
       3.0,
       3.5,
       4.0,
       4.8,
       5.6,
       6.6,
       7.6)

# use minitad treshold
shifted_data <- x - 0.51634

# setup fit 
init_shape <- mean(shifted_data) ^ 2 / var(shifted_data)
init_scale <- var(shifted_data) / mean(shifted_data)
init_params <- c(shape = init_shape, scale = init_scale)
lower_bounds <- c(1e-4, 1e-4)
upper_bounds <- c(Inf, Inf)

# fit
fit <- optim(
  par = init_params,
  fn = \(params, data) {
    shape <- params[1]
    scale <- params[2]
    - sum(dgamma(
      data,
      shape = shape,
      scale = scale,
      log = TRUE
    ))
  },
  data = shifted_data,
  method = "L-BFGS-B",
  lower = lower_bounds,
  upper = upper_bounds,
  control = list(fnscale = 1)
)

fit$par
#  shape    scale 
#2.370547 1.191143

Unfortunately there is no do this and it is solved answer. You need to take care of the problematic numerical optimization for the threshold near the minimum of the data.
Then the problem when the optimization hits a shape value smaller or equal 1 has to be addressed in your fitting routine.

I don't know if the Minitab support is willing to give you more information on how they do the optimization and what steps they perform before estimating all 3 parameters (if they estimate them simultaneously and not using multi step procedures).

Hope it gives you a good starting point for further investigation.