Maximum Probability - Incomplete Gamma Law
07-27-2019, 01:32 PM
Post: #1 Eddie W. Shore Senior Member Posts: 956 Joined: Dec 2013
Maximum Probability - Incomplete Gamma Law
Blog Post: http://edspi31415.blogspot.com/2019/07/h...ximum.html

The program IGLMAX calculates four calculation points of the Incomplete Gamma Law:

Three parameters of A, γ, and β of the IGL (Incomplete Gamma Law):

g(x) = 1 / (β^γ * gamma(γ)) * x^(γ - 1) * e^(x/β)

where:

X = list of data points, where x_i ≥ 0
s = number of data points
r = number of points where x_i = 0 (zero points)

A = ln(mean(X)) - (Σ ln(X_i))/(s - r)
γ = 1/(4*A) * (1 + √(1 + 4*A/3))
β = mean(X)/γ

And

p = probability that x is not exceeded
p = r/s + (1 - r/s) * (1 - uigf(γ, x)/gamma(γ))

Gamma Function:
gamma(γ) = ∫( t^(γ - 1) * e^(-t) dt, 0, ∞ )

Upper Incomplete Gamma Function:
uigf(γ, x) = ∫( t^(γ - 1) * e^(-t) dt, x/β, ∞ )

One particular application is determining the maximum limit that rainfall exceeds x (inches or mm). The book "Pocket Computers in Agrometeorology" introduces this concept and provides a program for the classic TI-59 (see source below).

HP Prime Program IGLMAX

Code:
EXPORT IGLMAX(L1,X) BEGIN // list of data ≥0, X // 2019-06-16 EWS LOCAL S,R,M,L,I; LOCAL A,Y,B,N,G,P; // set up S:=SIZE(L1); R:=0; // count zeros M:=0; // ΣX L:=0; // Σ(LN(X)) // counting loop FOR I FROM 1 TO S DO IF L1(I)==0 THEN R:=R+1; ELSE M:=M+L1(I); L:=L+LN(L1(I)); END; END; // parameters A:=LN(M/(S-R))-L/(S-R); Y:=(4*A)^(−1)*(1+√(1+4*A/3)); B:=M/(S-R)*1/Y; // gamma G:=CAS.Gamma(Y); // upper incomplete gamma N:=∫(T^(Y-1)*e^(−T),T,X/B,∞); // maximum probability P:=R/S+(1-R/S)*(1-N/G); // results RETURN {"A=",A,"γ=",Y, "β=",B,"Max Prob=",P}; END;

Example:

Data from a city of the rainfall in 2017 and 2018 (in inches):

2017
January: 3.90
February: 2.84
March: 2.31
April: 0.98
May: 0.64
June: 0.05
July: 0.00
August: 0.01
September: 0.00
October: 0.33
November: 0.72
December: 1.08

2018
January: 2.49
February: 2.66
March: 3.06
April: 2.94
May: 2.33
June: 0.81
July: 0.05
August: 0.00
September: 0.00
October: 0.14
November: 0.50
December: 2.24

Parameters:
A: 0.7237035089
γ (Gamma): 0.8296776362
β (Beta): 1.812752248

Probability that X inches of rainfall will not exceed:
X = 1 in: 0.593857
X = 2 in: 0.781173
X = 3 in: 0.879613

Source:

R.A. Gommes "Pocket Computers In Agrometeorology" Food and Agriculture Organization of the United Nations. FAO PLANT PRODUCTION AND PROTECTION PAPER. Rome, 1983. ISBN 92-5-101336-5
 « Next Oldest | Next Newest »

User(s) browsing this thread: 1 Guest(s)