**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.

## No comments:

## Post a Comment