Fit data to the two parameter gamma distribution.
Available in version 6.1.0 and later.
function dim_gamfit_n ( x : numeric, ; float, double optgam : logical, dims [*] : integer ) return_val : float or double
A variable of type 'float' or 'double' and any dimensionality. Since a distribution is being fit, there should be a 'reasonably' large sample size (say, at least 30 values).optgam
If optgam=False, an estimate of the returned parameters will be provided regardless of the number of non-missing (_FillValue) observations.
If optgam=True, the following attributes can be associated with optgam:
- optgam@inv_scale=True, the returned 'scale' value will be (1/scale). The default is optgam@inv_scale=False.
- optgam@pcrit specifies the minimum percent of observations
required for parameters to be returned. The default is pcrit=0.
So, if it is desired that (say) 75% of the observation must be nonmissing, then
optgam = True optgam@pcrit = 75
The dimension(s) of x on which to estimate the distribution parameters. Usually, this is the record ('time') dimension.
The estimated parameters returned will be the same type as x. The leftmost dimension will be size 3. The order will be  shape;  scale;  unconditional probability ( 0.0 <=> 1.0 ) of a 0.0 observation).
The dim_gamfit_n will return estimates of the shape and scale of the 2-parameter gamma distribution. These estimates are derived using maximum likelihood as described in
A Note on the Gamma Distribution Thom (1958): Monthly Weather Review, pp 117-122. specifically: eqn 22 for gamma; just above eqn 21 http://water.usgs.gov/osw/bulletin17b/MWR_Thom_1958.pdfIn addition, a third value, the unconditional probability of encountering an observation equal to zero is returned. The latter parameter has nothing to do with the gamma distribution. However, this function was included to analyze precipitation observations and this 3rd parameter is useful for certain applications.
The 2-parameter gamma distribution is modeled via:
gamma2_distribution = (x^(shape-1) * exp( -x/scale )) / (scale^shape * gamma(shape))Note_1: The scale is returned as [1/scale] in some applications. This is called the rate.
Note_2: As an aside, the scale and shape parameters for the classic two-parameter gamma distribution can be estimated via the "method of moments".
m1 = avg(x) ; 1st moment m2 = avg(x^2) ; 2nd moment or, if (say) x(time,lat,lon) m1 = dim_avg_n(x,0) ; 1st moment m2 = dim_avg_n(x^2,0) ; 2nd moment then shape = m1^2 / (m2 - m1^2) ; gamma scale = (m2 - m1^2) / m1 ; beta ---------- Zhou Xin (11 Aug 2015) provided the following information: I looked at the source codes used by dim_gamfit_n() and cdfgam_p(). If I understood it correctly, the "SCALE" parameter in dim_gamfit_n() means "beta" in the Thom (1958) paper Eqn (1), while the "SCALE" parameter in cdfgam_p() means "1/beta". So when using the output of dim_gamfit_n() for cdfgam_p(), the "SCALE" parameter has to be inverted.
; The following data were from: http://www.wessa.net/rwasp_fitdistrgamma.wasp g = (/112,118,132,129,121,135,148,148,136,119,104,118 \ ,115,126,141,135,125,149,170,170,158,133,114,140 \ ,145,150,178,163,172,178,199,199,184,162,146,166 \ ,171,180,193,181,183,218,230,242,209,191,172,194 \ ,196,196,236,235,229,243,264,272,237,211,180,201 \ ,204,188,235,227,234,264,302,293,259,229,203,229 \ ,242,233,267,269,270,315,364,347,312,274,237,278 \ ,284,277,317,313,318,374,413,405,355,306,271,306 \ ,315,301,356,348,355,422,465,467,404,347,305,336 \ ,340,318,362,348,363,435,491,505,404,359,310,337 \ ,360,342,406,396,420,472,548,559,463,407,362,405 \ ,417,391,419,461,472,535,622,606,508,461,390,432 /) ; Fit the 2-parameter (shape and scale) gamma distribution via the "method of moments" g1 = avg(g) g2 = avg(g^2) gshape = g1^2 / (g2 - g1^2) ; alpha gscale = (g2 - g1^2) / g1 ; beta print("gshape = "+gshape) print("gscale = "+gscale) ; Estimate the 2-parameter gamma distribution via maximum likelihood gpar = dim_gamfit_n(g, False, 0) ; gpar(3) print(gpar)The online results, which returned rate (=1/scale) were:
shape=gamma=5.459, rate=0.019476 ==> scale= (1/rate) = 51.345NCL's (edited) results from dim_gamfit_n are
; Method of Moments (0) gshape = 5.4973 ; shape (0) gscale = 50.9884 ; scale ; dim_gamfit_n (0) gpar = 5.49911 ; shape (0) gpar = 50.9716 ; scale (0) gpar = 0.0 ; pzeroExample 2
If x(time,lat,lon) and it is desired that 60% of the values be present at each grid point before the parameters are calculated. Also, return the inverse scale (1/scale).
gopt = True gopt@pcrit = 60.0 gopt@inv_scale = True xpar = dim_gamfit_n(x, gopt, 0) ; xpar(3,nlat,mlon)