fit gamma distribution based on median interquartile range

Question

I would like to create a gamma probability distribution in R with a median of 3 and interquartile range of [2,5]. I am familiar with the methods of moments to estimate the scale and shape of a gamma distribution but this is based on the mean and standard deviation and both are not available. Is there a way to estimate a gamma distribution in R with a specific median and interquartile range?

Answer 1

A gamma distribution with given median and IQR

It is not easy because there is no closed form expression for the median and quartiles. For the median you can find several approximate formula's. For the quartiles you could use R's approximation of the quantile function, qgamma, and then use an iterative algorithm that finds $k$ and $\theta$ that suit your target.

One way to make such algorithm can be by optimizing a cost function

f = function(par, targetmedian, targetIQR) {
  k = par[1]
  theta = par[2]

  median = theta*k*(1-1/9/k)^3
  ### alternatively compute the median like below
  #median = qgamma(0.5,shape = k,scale = theta)

  IQR = qgamma(0.75,shape = k,scale = theta) - qgamma(0.25,shape = k,scale = theta)

  error = (median-targetmedian)^2+(IQR-targetIQR)^2
  return(error)
}

### optim function below gives result k = 2.905083 theta = 1.160847
p2 = optim(par  = c(1,1), f, targetmedian = 3, targetIQR = 2.5)$par

### check median and IQR
qgamma(0.5, shape = p2[1], scale = p2[2]) ### 2.994302
qgamma(0.75,shape = p2[1], scale = p2[2])-qgamma(0.25,shape = p1[1],scale = p1[2]) # 2.494646

A gamma distribution with given quartiles

In this case you do not have two, but three values, whereas the gamma distribution only has two parameters. It might be possible that an exact fit is not possible. Now the cost function for the optimization becomes more important.

• You could use the Kolmogorov statistic (the distance between the cdf and your quartiles that define a empirical distribution)

• You can consider data bins $x\leq 2$, $2<x\leq 3$, $3 < x \leq 5$ and $5<x$ with each 25% probability and optimize the cross-entropy.

example code:

f = function(par, target_quartiles, fit_method = "Kolmogorov") {
  k = par[1]
  theta = par[2]

  if (fit_method %in% c("Kolmogorov", "cross-entropy") == FALSE) { 
    stop("use method Kolmogorov or cross-entropy")
  }

  p_quartiles = pgamma(target_quartiles, shape = k, scale = theta)
  
  if (fit_method == "Kolmogorov") {
     statistic = max(abs(p_quartiles - c(0.25,0.5,0.75)))
  }
  if (fit_method == "cross-entropy") {
     p1 = diff(c(0,p_quartiles,1))
     p2 = c(0.25,0.25,0.25,0.25)
     statistic = -sum(p2*log(p1))
  }     
  return(statistic)
}   

### optim function below gives result k = 2.399316 theta = 1.502534
par1 = optim(par  = c(1,1), f, target_quartiles = c(2,3,5))$par

### 0.2724194 0.4775805 0.7724194
pgamma(c(2,3,5), shape = par1[1], scale = par1[2])

### optim function below gives result k = 2.408474 theta = 1.530167
par2 = optim(par  = c(1,1), f, target_quartiles = c(2,3,5), fit_method = "cross-entropy")$par

### 0.2625420 0.4645755 0.7610788
pgamma(c(2,3,5), shape = par2[1], scale = par2[2])