Saturday, July 27, 2019

HP Prime and TI-84 Plus CE: Maximum Probability - Incomplete Gamma Law

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.