CODEWITHSUNDEEP
rstatisticsgamma-distribution
user18441
user18441

Reputation: 663

How can I estimate the shape and scale of a gamma dist. with a particular mean and a 95% quantile?

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?

Upvotes: 5

Views: 2449

Answers (2)

Pariksheet Nanda
Pariksheet Nanda

Reputation: 75

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

Upvotes: 2

Greg Snow
Greg Snow

Reputation: 49640

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.

Upvotes: 3

Related Questions