HP Prime and TI-84 Plus CE: Maximum Probability - Incomplete Gamma Law
Introduction
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
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;
TI-84 Plus Program IGLMAX
(Text to enter)
Note: probability is rounded to six decimal places
"EWS 2019-06-16"
Input "DATA (X≥0):",L1
Input "X:",X
"INITIALIZE"
dim(L1)→S
0→R
0→M
0→L
"COUNTING LOOP"
For(I,1,S)
If L1(I)=0
Then
R+1→R
Else
M+L1(I)→M
L+ln(L1(I))→L
End
End
"PARAMETERS"
ln(M/(S-R))-L/(S-R)→A
(4*A)^(1)*(1+√(1+4*A/3))→Y
M/(S-R)*1/Y→B
"GAMMA"
fnInt(T^(Y-1)*e^(T),T,0,500)→G
"INCOMPLETE GAMMA"
fnInt(T^(Y-1)*e^(T),T,X/B,500)→N
"PROB"
R/S+(1-R/S)*(1-N/G)→P
round(P,6)→P
ClrHome
Disp "A=",A,"GAMMA=",Y,"BETA=",B
Pause
Disp "MAX PROB (FIX 6)=",P
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
Eddie
All original content copyright, © 2011-2019. Edward Shore. Unauthorized use and/or unauthorized distribution for commercial purposes without express and written permission from the author is strictly prohibited. This blog entry may be distributed for noncommercial purposes, provided that full credit is given to the author.
Introduction
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
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;
TI-84 Plus Program IGLMAX
(Text to enter)
Note: probability is rounded to six decimal places
"EWS 2019-06-16"
Input "DATA (X≥0):",L1
Input "X:",X
"INITIALIZE"
dim(L1)→S
0→R
0→M
0→L
"COUNTING LOOP"
For(I,1,S)
If L1(I)=0
Then
R+1→R
Else
M+L1(I)→M
L+ln(L1(I))→L
End
End
"PARAMETERS"
ln(M/(S-R))-L/(S-R)→A
(4*A)^(1)*(1+√(1+4*A/3))→Y
M/(S-R)*1/Y→B
"GAMMA"
fnInt(T^(Y-1)*e^(T),T,0,500)→G
"INCOMPLETE GAMMA"
fnInt(T^(Y-1)*e^(T),T,X/B,500)→N
"PROB"
R/S+(1-R/S)*(1-N/G)→P
round(P,6)→P
ClrHome
Disp "A=",A,"GAMMA=",Y,"BETA=",B
Pause
Disp "MAX PROB (FIX 6)=",P
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
Eddie
All original content copyright, © 2011-2019. Edward Shore. Unauthorized use and/or unauthorized distribution for commercial purposes without express and written permission from the author is strictly prohibited. This blog entry may be distributed for noncommercial purposes, provided that full credit is given to the author.