1
$\begingroup$

I have this dataframe below:

structure(list(x = c(0, 1, 2, 3, 8, 0, 1, 0, 1, 2, 48, 49, 50, 
51, 52, 53, 54, 55, 56, 1, 2, 3, 4, 5, 13, 14, 15, 16, 17, 0, 
71, 72, 73, 74, 75, 37, 38, 39, 40, 41, 20, 21, 22, 23, 24, 12, 
13, 14, 15, 16, 4, 5, 6, 7, 80, 81, 82, 83, 84, 85, 86, 1, 2, 
3, 4, 5, 6, 88, 89, 90, 91, 26, 27, 28, 14, 15, 16, 5, 6, 7, 
7, 8, 9, 10, 11, 1, 3, 29, 30, 0, 1, 20, 21, 22, 2, 3, 0, 1, 
2, 3, 4, 0, 1, 2, 3, 4, 0, 1, 2, 3, 1, 2, 1, 2, 97, 98, 99, 100, 
101, 102, 103, 104, 105, 106, 83, 16, 18, 19, 86, 87, 91, 92, 
93, 20, 21, 99, 100, 102, 103, 9, 10, 14, 0, 1, 2, 3, 4, 5, 6, 
9, 10, 11, 12, 13, 14, 79, 80, 81, 82, 83, 0, 1, 2, 3, 4, 5, 
7, 11, 21, 63, 64, 65, 66, 3, 18, 19, 20, 21), y = c(1665.05355793238, 
2290.46808761358, 2479.44521586597, 1027.63470042497, 2382.83166897297, 
599.608228149824, 200.079963575583, 2715.53603386879, 2567.91178971529, 
2709.88102674484, 454.382090631872, 2268.14043972082, 2147.62290922552, 
2269.1387550775, 2247.31983098201, 1903.39138268307, 2174.78291538358, 
2359.51909126411, 2488.39004804939, 212.851575751527, 461.398994384333, 
567.150629704352, 781.775113821961, 918.303706148872, 659.580393221229, 
1107.37695799186, 1160.80594193377, 1412.61328924168, 1689.48879626486, 
650.35140038206, 280.949480787385, 641.231373456365, 412.116433330579, 
416.824746264203, 614.763061881065, 415.905685925856, 494.374217984441, 
201.745910386788, 486.030122926459, 647.782697262242, 269.898762604455, 
469.544556583467, 605.790549736819, 483.045965372879, 668.017897433514, 
35.2706101682852, 265.693628564011, 285.116345260642, 291.023782284086, 
357.428790589795, 205.920375034591, 229.606221692753, 230.952761338012, 
241.641164634028, 494.819531862799, 611.564698476112, 603.238720579422, 
524.322001078981, 565.296378873638, 532.431853589369, 597.174114277044, 
140.829281608807, 260.737164468854, 306.72700499362, 283.410379620422, 
366.813913489692, 387.570173754128, 648.038833109662, 606.075737104722, 
686.408686154056, 705.914347674276, 814.318649679422, 720.467501226114, 
897.123468788341, 1015.63389401929, 892.041182056069, 1587.7192387376, 
469.889315964654, 371.252178900875, 516.214449755847, 223.317826933191, 
233.990588575922, 254.111292146219, 388.602676983443, 477.858510450125, 
53.9083766029653, 128.198042456082, 521.591475207009, 535.519377609133, 
560.44362318865, 396.430340241175, 654.226210802794, 1303.48454248847, 
1274.19490102748, 1615.11958998814, 1445.12280311435, 1579.06517848652, 
1591.83073895803, 1926.41506405128, 2029.74201928871, 1762.10456025088, 
1093.34139243141, 998.179578598123, 980.315594917862, 1016.85559533851, 
996.357576390554, 1465.61544544157, 1639.27601695352, 3099.03079020977, 
3250.79935610294, 613.15010408836, 794.168323367026, 953.81404011324, 
1942.00127470493, 750.459735572338, 875.899619173259, 1166.28903575242, 
788.564458085224, 1324.95006324351, 1111.05266567506, 1020.29043261707, 
1130.40108603239, 1529.81020075083, 1480.88505926495, 524.398260936374, 
512.986586335697, 833.124012635089, 957.75556050241, 220.735156724928, 
584.620856728405, 1243.39788535796, 1578.20351953804, 1386.65165398503, 
906.813018662913, 798.186995701282, 831.365377249207, 764.519073183124, 
672.076289062505, 669.879217186302, 1388.61752325296, 1104.50366285443, 
1722.3172301054, 313.881701797247, 1217.21950829029, 1325.65616729856, 
2137.73167350888, 1958.17336059175, 2355.88909001648, 2261.16462984681, 
2499.84794190526, 908.743901610374, 1927.14038875327, 1880.97180615366, 
1977.65100651979, 1714.75837931037, 1265.48484068702, 1193.29000986667, 
1156.81486114406, 1199.7373066445, 1116.24029749935, 1141.89299393632, 
1545.56326286064, 1401.44881976186, 1346.97588634375, 1640.27575962036, 
1458.74134225518, 1102.85304791946, 1617.68811492418, 1380.66626644135, 
1614.92309916019, 1590.29306322336, 1516.55950581282, 1616.58492165804, 
785.721945669502, 1857.46475183964, 3546.82112717628, 3451.61194658279, 
2874.73860728741)), .Names = c("x", "y"), class = "data.frame", row.names = c(186L, 
187L, 188L, 189L, 190L, 191L, 192L, 193L, 194L, 195L, 196L, 197L, 
198L, 199L, 200L, 201L, 202L, 203L, 204L, 206L, 207L, 208L, 209L, 
210L, 211L, 212L, 213L, 214L, 215L, 216L, 217L, 218L, 219L, 220L, 
221L, 222L, 223L, 224L, 225L, 226L, 227L, 228L, 229L, 230L, 231L, 
232L, 233L, 234L, 235L, 236L, 237L, 238L, 239L, 240L, 241L, 242L, 
243L, 244L, 245L, 246L, 247L, 248L, 249L, 250L, 251L, 252L, 253L, 
254L, 255L, 256L, 257L, 258L, 259L, 260L, 261L, 262L, 263L, 264L, 
265L, 266L, 267L, 268L, 269L, 270L, 271L, 272L, 273L, 274L, 275L, 
276L, 277L, 278L, 279L, 280L, 281L, 282L, 283L, 284L, 285L, 286L, 
287L, 288L, 289L, 290L, 291L, 292L, 293L, 294L, 295L, 296L, 298L, 
299L, 301L, 302L, 303L, 304L, 305L, 306L, 307L, 308L, 309L, 310L, 
311L, 312L, 313L, 314L, 315L, 316L, 317L, 318L, 319L, 320L, 321L, 
322L, 323L, 324L, 325L, 327L, 328L, 329L, 330L, 334L, 335L, 336L, 
337L, 338L, 339L, 340L, 341L, 342L, 343L, 344L, 345L, 346L, 347L, 
348L, 349L, 350L, 351L, 352L, 353L, 354L, 355L, 356L, 357L, 358L, 
359L, 360L, 361L, 362L, 363L, 364L, 365L, 366L, 367L, 368L, 369L, 
370L))

I want to fit a gamma distribution to this dataset. I am trying to use this function rgamma(n=, shape =, scale= ). However I do not quite understand the way to compute the different parameters of the gamma distribution. I am quite a dummy in this and the examples I have seen so far do not really explain how to get the shape and scale values. Can someone give me some insights on it?

$\endgroup$
7
  • 1
    $\begingroup$ This is a bivariate dataset, whereas a Gamma distribution is univariate only. What, then, do you mean by "fit a gamma distribution" to it? And why do you want to do that? $\endgroup$
    – whuber
    Jul 7, 2015 at 14:42
  • $\begingroup$ +1 to the comment above. Also the rgamma function generates random values from the gamma distribution so that would not be the function you would want to "fit" the data. $\endgroup$ Jul 7, 2015 at 14:44
  • 1
    $\begingroup$ What you are really asking, then, is how to fit a function of the form $y = A x^k\exp(-x)$ to your $(x,y)$ data. The Gamma distribution appears to have nothing to do with it. What is more relevant is how you need to assess the goodness of fit. If these points don't exactly line up along one of those curves, why would that occur? What is the source and nature of the variability that could produce such discrepancies? $\endgroup$
    – whuber
    Jul 7, 2015 at 15:03
  • 1
    $\begingroup$ I plotted the data. There is no aspect of the $(x,y)$ data that looks remotely like a Gamma-shaped curve. The distribution of $y$ might possibly have a Gamma shape. If that's what you want to fit, then it would be best to eliminate the superfluous $x$ variable from your example. In any event, it's essential that you edit this question to clarify what you're trying to do, so that people don't get frustrated trying to figure out what you need. $\endgroup$
    – whuber
    Jul 7, 2015 at 15:23
  • 1
    $\begingroup$ A simple plot(A$x, A$y) (where A is your structure) shows that you may have some left-tail dependance, but nothing about the shape of the bivariate distribution. Yes, the very "uppermost" points in that plot seem to "trace" out the classic shape of a right-skewed density, but that doesn't mean anything; it's an optical illusion. $y$ itself may be gamma (see below in the answer and comments). Marginal $x$ certainly is not (it's bimodal). $\endgroup$
    – Avraham
    Jul 7, 2015 at 15:24

1 Answer 1

9
$\begingroup$

Please see update section below for clarification based on OP comments made whilst this answer was being posted.

In general, the standard gamma distribution has two parameters, a shape paramemeter and a scale or rate paramemeter, where rate = 1 / scale. R help for rgamma uses shape as the canonical form; using that, with a gamma density defined as:

$$ \large f(x) = \frac{x^{\alpha - 1}e^{-\frac{x}{\theta}}}{\theta^\alpha \Gamma(\alpha)} $$

$\alpha$ is the shape parameter and $\theta$ is the scale parameter. The mean of the distribution is $\alpha\theta$ and the variance is $\alpha\theta^2$.

The simplest way to fit would be using the Method of Moments. Calculate the mean and variance of your data. Now you have two equations ($\hat{\mu} = \alpha\theta, \hat{\sigma}^2 = \alpha\theta^2$) so you can solve for $\alpha$ and $\theta$.

The other most common technique is tp use maximum likelihood estimation. That would require using a non-linear optimizer (like optim or nloptr) to solve for the parameters which minimimze the negative log-likelihood, which you could get using dgamma with log = TRUE.

You could set up a function like so:

GammaNLL <- function(pars, data){
  alpha <- pars[[1]]
  theta <- pars[[2]]
  return (-sum(dgamma(data, shape = alpha, scale = theta, log = TRUE)))
}

and pass that into nloptr as

Fit <- nloptr(x0 = c(guess1, guess2), eval_f = GammaNLL, lb = c(0,0), data = data,
              opts = list(algorithm = "NLOPT_LN_SBPLX", maxeval = 1e5))

Fit$solution will have your MLE estimates.

In your particular case, however, your observations are pairs, so either you need to fit a multivariate distribution, which is less simple, or you want to fit the gamma to $x$ and $y$ independantly, which can be done using either method above.

Hope that helps.

Update based on clarifications in question.

Plotting $y$ against $x$ looks as such plot(A$x, A$y, pch =16):

enter image description here

The fact that the "top" lookes vaguely like a gamma is irrelevant; this is a scatter plot, not a density.

A kernel-smoothed density estimates of the marginal of $y$ may indicate a gamma (or some other right-skewed distribution) plot(density(A$y, kernel = 'e')):

enter image description here

$x$ is bimodal, so it's marginal is certainly not gamma plot(density(A$x, kernel = 'e')): enter image description here

So I think you need to re-formulate your question more like "is there a reasonable way to fit a relationship between $x$ and $y$"

$\endgroup$
6
  • $\begingroup$ This approach is hopeless: because the integral is much greater than $1$, the data do not represent points on a Gamma distribution at all and no Gamma density function comes anywhere close to such a set of values. $\endgroup$
    – whuber
    Jul 7, 2015 at 15:05
  • $\begingroup$ @whuber The marginal $y$ can be fit using either method (MoM gives $\theta \approx 508.9, \alpha \approx 2.17$; MLE: $\theta \approx 546.2, \alpha \approx 2.02$). The $x$ values are clearly bimodal though, so any fit wouldn't be good, although MoM will give an answer, simply due to algebra. I did say the method wouldn't work if there was a need for a multivariate fit. $\endgroup$
    – Avraham
    Jul 7, 2015 at 15:13
  • $\begingroup$ For clarification, please read the OP's comments to the question. $\endgroup$
    – whuber
    Jul 7, 2015 at 15:14
  • $\begingroup$ @whuber They must have been added as I was answering, as the question was pristine when I started. Do you suggest deleting this answer? $\endgroup$
    – Avraham
    Jul 7, 2015 at 15:15
  • 1
    $\begingroup$ I think some kind of modification would be needed to get the answer to match the question and the data. It is indeed frustrating to see a question changed while you're writing out an answer. That's why it's always a good idea to ask for clarification if you're uncertain about the interpretation or to wait for a response after a clarification has been requested. One thing you might do is delete the answer so that people don't downvote it--you can undelete it later if it turns out to be exactly what the OP was looking for. $\endgroup$
    – whuber
    Jul 7, 2015 at 15:18

Not the answer you're looking for? Browse other questions tagged or ask your own question.