A gamma distribution with given median and IQR
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])