4

Is there any way, in R, to calculate the scale and shape of a gamma distribution, given a particular value of mean (or median) and a particular quantile (the 95% quantile)?

So for example I have a mean = 130

and a 95% quantile = 300

with an offset of the distribution at 80

is there any way to obtain the scale and shape of a gamma that meet these criteria?

mnel
  • 113,303
  • 27
  • 265
  • 254
user18441
  • 643
  • 1
  • 7
  • 15

2 Answers2

3

Here is one approach:

myfun <- function(shape) {
    scale <- 130/shape
    pgamma(300, shape, scale=scale) - 0.95
}

tmp <- uniroot( myfun, lower=2, upper=10 )

myshape <- tmp$root
myscale <- 130/tmp$root

qgamma(0.95, shape=myshape, scale=myscale)
integrate( function(x) x*dgamma(x,shape=myshape,scale=myscale), 
    lower=0, upper=Inf )

I am not sure what you mean by offset of 80, if that is just where the gamma becomes non-zero then subtract 80 from 130 and 300 and do the same as above.

Greg Snow
  • 48,497
  • 6
  • 83
  • 110
  • thank you, it seems though that for some combinations of numbers (for example: mean =8, quartile 95% = 32) this isnt working, it gives an error such as: Error in uniroot(myfun, lower = 2, upper = 10) : f() values at end points not of opposite sign – user18441 Jan 11 '13 at 22:54
  • @user18441, you need to give it starting values that bracket the answer, meaning one of them returns a postive value from `myfun` and the other a negative. Try a few values and see what you get, then choose 2 for the starting values. See `?uniroot` for more detail. – Greg Snow Jan 12 '13 at 18:47
1

Use the optim() function to simultaneously fit the two parameters of the distribution:

offset <- 80
x <- c(q50 = 130, q95 = 300) - offset # Subtract offset and re-add later.
expected <- c(0.5, 0.95)

gamma_fit_error <- function(params) {
  with(as.list(params), {
    mat <- t(matrix(c(x,
                      qgamma(expected, shape = shape, scale = scale)),
                    nrow = length(x)))
    c(dist(mat))
  })
}

opt <- optim(par = c(shape = 1, scale = 10),
             fn = gamma_fit_error,
             method = "L-BFGS-B",
             lower = c(shape = 0.1, scale = 1),
             upper = c(shape = 100, scale = 100))

The fn input to optim() requires a single params input from which we need to unpack the individual parameters; therefore with() unpacks shape and scale from the named vector params input.

gamma_fit_error() returns the distance between the data x at the 50% and 95% quantiles and the calculated values at those quantiles using the gamma distribution parameters. dist() finds the distance between the rows, which is why one needs to t() transpose the values from column-wise matrix to row-wise. Because dist() returns a matrix, we use the c() shorthand to convert it into a vector.

Outputs:

> opt
$par
     shape      scale 
 0.9765883 74.5730641 

$value
[1] 0.007629429

$counts
function gradient 
     250      250 

$convergence
[1] 0

$message
[1] "CONVERGENCE: REL_REDUCTION_OF_F <= FACTR*EPSMCH"

> with(as.list(opt$par), pgamma(x, shape = shape, scale = scale))
      q50       q95 
0.4999896 0.9499950

> with(as.list(opt$par), offset + qgamma(expected, shape = shape, scale = scale))
[1] 130.0015 300.0075

> x + offset
q50 q95 
130 300

> hist(with(as.list(opt$par), offset + rgamma(1e4, shape = shape, scale = scale)))

histogram of the gamma distribution

Pariksheet Nanda
  • 65
  • 8