Maximum Probability - Incomplete Gamma Law

[Eddie W. Shore]

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