estimating percentiles of gamma distribution

[Lea Nassal]

Hi,
I'm a new stata user and currently I'm trying to fit a gamma distribution to my data and compute the corresponding percentiles of the estimated distribution. I used the gammafit command to estimate the parameters of the distribution. I know that in R there exists an command "qgamma" which computes the percentiles of the gamma distribution. Now my question is if there exists a similar command in stata.

[Joseph Coveney]

Welcome to Statalist. You can start by typing the following at the command line in stata.

help gammap()

In general, statistical functions (densities, quantiles etc. and their inverses for various distributions) can be found by

help functions

and then clicking on the hyperlink for statistical functions in the screen that pops up.

[Nick Cox]

gammafit is from SSC, as you are asked to explain (FAQ Advice #12: explain where community-contributed commands you refer to come from). Also from SSC is qgamma, which will be found to use

`B' * invgammap(`A',`Psubi')

as the recipe for a predicted gamma quantile for cumulative probability p_i and estimated parameters A = alpha and B = beta

[Anh Trinh]

Dear Statalist,

I have a rainfall variable and I would like to fit a separate gamma distribution to each country’s time series; and within each country each year is assigned its corresponding percentile in its gamma distribution.

* Example generated by -dataex-. To install: ssc install dataex
clear
input byte country int year double rainfall
1 2005 1890.3680057525635
1 2006 1527.0653429031372
1 2007 1768.3753786087036
1 2008 2306.7355098724365
1 2009 1594.1028535366058
1 2010 1502.7944736480713
1 2011 1917.2996888160706
1 2012 2004.7223992347717
1 2013 1872.8261060714722
1 2014 1660.2155447006226
1 2015 1539.1846103668213
1 2016 1656.9323408603668
1 2017 1775.4609567523003
2 2005  1986.130392074585
2 2006 1900.2506830692291
2 2007 1841.2762060165405
2 2008   2161.82772731781
2 2009 1687.9958666563034
2 2010  1737.955777645111
2 2011 1773.8970966339111
2 2012 1850.1043663024902
2 2013 1876.7550468444824
2 2014   1925.95050573349
2 2015 2012.2574634552002
2 2016 1771.2601299285889
2 2017 2206.6164360046387
end

I would appreciate any help.
Thank you!

[Nick Cox]

Here's a sketch. Code not really tested. See gammafit (SSC) as mentioned in #1 and #3.

gen alpha = . 
gen beta = . 
gen P = . 

su country, meanonly 

quietly forval c = 1/`r(max)' { 
    capture gammafit rainfall if country == `c'
    if _rc = 0 { 
        replace alpha = e(alpha) if country == `c' 
        replace beta = e(beta) if country == `c'
        replace P = gammap(e(alpha), rainfall/e(beta)) if country == `c'
    }
    else di "country `c': no go " _rc 
}

[Anh Trinh]

Thank you, Nick! The code works well, just a minor typo that "if _rc==0".