code 41000;
real procedure BIN(X, N, P);
value X, N, P; real X, N, P;
begin integer IX;
IX:= ENTIER(X);
if N < 0 ∨ N > ENTIER(N) then
BIN:= STATAL3 ERROR(“BIN”, 2, N) else
if P < 0 ∨ P > 1 then
BIN:= STATAL3 ERROR(“BIN”, 3, P) else
if IX < 0 then BIN:= 0 else
if IX ≥ N then BIN:= 1 else
if P = 0 then BIN:= 1 else
if P = 1 then BIN:= 0 else
if N > 1000 then
begin real B;
B:= 1 - INCOMPLETE BETA(P, IX + 1, N - IX, 10-12);
BIN:= if B < 0 then 0 else B;
end else
begin real procedure TAIL;
begin integer I; real PROB, CUM, LAST;
PROB:= CUM:= BINPROB(IX, N, P);
for I:= IX - 1, I - 1 while CUM > LAST do
begin PROB:=
PROB × (1 - P) / P × (I + 1) / (N - I);
LAST:= CUM; CUM:= CUM + PROB
end;
TAIL:= CUM
end;
if X > ENTIER(N / 2) then
begin IX:= N - IX - 1; P:= 1 - P;
BIN:= 1 - TAIL
end else BIN:= TAIL
end;
end BIN;
eop
code 41001;
real procedure BININV(PROB, N, P, LEFT);
value PROB, N, P, LEFT;
real PROB, N, P; Boolean LEFT;
begin integer X; real PX, PCUM;
if PROB ≤ 0 ∨ PROB ≥ 1 then
STATAL3 ERROR(“BININV”, 1, PROB) else
if N > ENTIER(N) ∨ N < 0 then
STATALS ERROR(“BININV”, 2, N) else
if P < 0 ∨ P > 1 then
STATAL3 ERROR(“BININV”, 3, P);
if P = 0 ∨ N = 0 then
BININV:= (if LEFT then -1 else 1)
else if P = 1 then
BININV:= (if LEFT then N - 1 else N + 1)
else if LEFT then
begin X:= PHINV(PROB) ×
SQRT(N × P × (1 - P)) - 0.5 + N × P;
if X < 0 then X:= 0
else if X > N then X:= N;
if PROB < (1 - P) ⭡ N then BININV:= -1 else
begin PX:= BINPROB(X, N, P); PCUM:= BIN(X, N, P);
if PCUM > PROB then
begin for PCUM:= PCUM - PX
while PCUM > PROB do
begin PX:= PX × X × (1 - P) /
(N - X + 1) / P;
X:= X - 1
end; X:= X - 1
end else
begin for PX:=
PX × (N - X) / (X + 1) × P / (1 - P)
while PCUM + PX < PROB do
begin X:= X + 1; PCUM:= PCUM + PX end
end;
BININV:= X
end
end else
begin X:= PHINV(1 - PROB) ×
SQRT(N × P × (1 - P)) + 0.5 + N × P;
if X < 0 then X:= 0 else
if X > N then X:= N;
if PROB < P ⭡ N then BININV:= N + 1 else
begin PCUM:= 1 - BIN(X - 1, N, P);
PX:= BINPROB(X, N, P);
if PCUM < PROB then
begin for PX:=
PX × X × (1 - P) / (N - X + 1) /P
while PCUM + PX < PROB do
begin X:= X - 1; PCUM:= PCUM + PX end
end else
begin for PCUM:= PCUM - PX
while PCUM > PROB do
begin PX:= PX × (N - X) × P /
(X + 1) / (1 - P);
X:= X + 1
end; X:= X + 1
end;
BININV:= X
end
end
end BININV;
eop
code 41004;
real procedure HYPERG(X, N, R, NN);
value X, N, R, NN; real X, N, R, NN;
begin integer I; real SUM, LAST, TERM; Boolean LEFT;
if N < 0 ∨ N > NN ∨ N - ENTIER(N) ≠ 0 then
STATAL3 ERROR(“HYPERG”, 2, N);
if R <0 ∨ R > NN ∨ R - ENTIER(R) ≠ 0 then
STATAL3 ERROR(“HYPERG”, 3, R);
if NN - ENTIER(NN) ≠ 0 then
STATAL3 ERROR(“HYPERG”, 4, NN);
LEFT:= true;
if N > NN / 2 then
begin LEFT:= false; N:= NN - N; X:= R - X - 1 end;
if R > NN / 2 then
begin LEFT:= ¬ LEFT; R:= NN - R;
X:= N - X - 1
end;
if N > R then begin I:= N; N:= R; R:= I end;
if X < 0 then HYPERG:= if LEFT then 0 else 1
else
if X ≥ N then HYPERG:= if LEFT then 1 else 0
else if NN > 105 then
begin real BETA, TAU, CHI;
TAU:= SQRT(R × N × (NN - N) × (NN - R) / NN) / NN;
CHI:= (X + .5 - N × R / NN) / TAU;
BETA:= (CHI × CHI + 2) / 12;
X:= if R ≤ NN / 4 then
2 × (SQRT((X + .5 + BETA)
× (NN - R - N + X + .5 + BETA))
- SQRT((N - X - .5 + BETA) ×
(R - X - .5 + BETA))) /
SQRT(NN + 1.5 - NN × NN / 2 / N / (NN - N))
else
CHI + (CHI × CHI - 1) ×
(2 × N - NN) × (NN - 2 × R)
/ 6 / TAU / NN / NN + CHI ×
(1 - 3 × (NN - N) × N / NN / NN)
/ 48 / TAU / TAU;
HYPERG:= PHI(if LEFT then X else -X)
end else
begin X:= ENTIER(X);
TERM:= SUM:= HYPERGPROB(X, N, R, NN);
if X > (N + 1) × (R + 1) / (NN + 2) then
begin LEFT:= ¬ LEFT; SUM:= 0;
for I:= X + 1, I + 1 while LAST < SUM do
begin TERM:= TERM × (N - I + 1) × (R - I + 1)
/ I / (NN - R - N + I);
LAST:= SUM; SUM:= SUM + TERM
end
end else
for I:= X, I - 1 while LAST < SUM do
begin TERM:= TERM × I × (NN - N - R + I)
/ (N - I + 1) / (R - I + 1);
LAST:= SUM; SUM:= SUM + TERM
end;
HYPERG:= if LEFT then SUM else 1 - SUM
end
end HYPERG;
eop
code 41005;
real procedure HYPERGINV(PROB, N, R, M, LEFT);
value PROB, N, R, M, LEFT;
real PROB, N, R, M; Boolean LEFT;
begin integer X; real PX, PCUM, LOW, UP;
if PROB ≤ 0 ∨ PROB ≥ 1 then
STATAL3 ERROR(“HYPERGINV”, 1, PROB) else
if N > ENTIER(N) ∨ N < 0 ∨ N > M then
STATAL3 ERROR(“HYPERGINV”, 2, N) else
if R > ENTIER(R) ∨ R < 0 ∨ R > M then
STATAL3 ERROR(“HYPERGINV”, 3, R) else
if M > ENTIER(M) ∨ M < 0 then
STATAL3 ERROR(“HYPERGINV”, 4, M);
LOW:= if N + R - M > 0 then N + R - M else 0;
UP:= if N < R then N else R;
if N = 0 ∨ R = 0 then
HYPERGINV:= (if LEFT then -1 else +1)
else if N = M ∨ R = M then
HYPERGINV:= (if LEFT then M - 1 else M + 1)
else if LEFT then
begin X:= PHINV(PROB) × SQRT((M - N) × N × R ×
(M - R) / (M × M × (M - 1))) + R × N / M + 0.5;
if X < LOW then X:= LOW else
if X > UP then X:= UP;
if PROB < HYPERGPROB(LOW, N, R, M) then
HYPERGINV:= LOW - 1
else
begin PX:= HYPERGPROB(X, N, R, M);
PCUM:= HYPERG(X, N, R, M);
if PCUM > PROB then
begin for PCUM:= PCUM - PX
while PCUM > PROB do
begin PX:= PX × X × (M - N - R + X) /
(N - X + 1) / (R - X + 1); X:= X - 1
end; X:= X - 1
end else
begin for PX:= PX × (N - X) × (R - X) /
(X + 1) / (R - X + 1)
while PCUM + PX < PROB do
begin X:= X + 1; PCUM:= PCUM + PX end
end;
HYPERGINV:= X
end
end else
begin X:= PHINV(1 - PROB) × SQRT((M - N) × N × R ×
(M - R) / (M × M × (M - 1))) + R × N / M - 0.5;
if X < LOW then X:= LOW else
if X > UP then X:= UP;
if PROB < HYPERGPROB(UP, N, R, M) then
HYPERGINV:= UP + 1
else
begin PCUM:= 1 - HYPERG(X - 1, N, R, M);
PX:= HYPERGPROB(X, N, R, M);
if PCUM < PROB then
begin for PX:= PX × X × (M - N - R + X) /
(N - X + 1) / (R - X + 1)
while PCUM + PX < PROB do
begin X:= X - 1; PCUM:= PCUM + PX end
end else
begin for PCUM:= PCUM - PX
while PCUM > PROB do
begin PX:= PX × (N - X) × (R - X) /
(X + 1) / (M - N - R + X + 1);
X:= X + 1
end; X:= X + 1
end;
HYPERGINV:= X
end
end
end HYPERGINV;
eop
code 41009;
real procedure NEGBIN(X, K, P);
value X, K, P; real X, K, P;
NEGBIN:= if K < 0 ∨ K > ENTIER(K)
then STATAL3 ERROR(“NEGBIN”, 2, K)
else if P ≤ 0 ∨ P > 1
then STATAL3 ERROR(“NEGBIN”, 3, P)
else if X ≥ K then
1 - BIN(K - 1, ENTIER(X), P)
else 0;
eop
code 41010;
real procedure NEGBININV(PROB, K, P, LEFT);
value PROB, K, P, LEFT;
real PROB, K, P; Boolean LEFT;
begin integer X; real PX, PCUM;
if PROB ≤ 0 ∨ PROB ≥ 1 then
STATAL3 ERROR(“NEGBININV”, 1, PROB) else
if K > ENTIER(K) ∨ K < 0 then
STATAL3 ERROR(“NEGBININV”, 2, K) else
if P ≤ 0 ∨ P > 1 then
STATAL3 ERROR(“NEGBININV”, 3, P);
if P = 1 ∨ K = 0 then
NEGBININV:= (if LEFT then K - 1 else K + 1)
else if LEFT then
begin X:= (PHINV(PROB) ×
SQRT(K × (1 - P)) + K - P / 2) / P;
if X < K then X:= K;
if PROB < P ⭡ K then NEGBININV:= K - 1 else
begin PX:= NEGBINPROB(X, K, P);
PCUM:= NEGBIN(X, K, P);
if PCUM > PROB then
begin for PCUM:= PCUM - PX
while PCUM > PROB do
begin PX:= PX × (X - K) / (1 - P)
/ (X - 1);
X:= X - 1
end; X:= X - 1
end else
begin for PX:= PX × (1- P) × X / (X - K + 1)
while PCUM + PX < PROB do
begin X:= X + 1; PCUM:= PCUM + PX end
end;
NEGBININV:= X
end
end else
begin X:= (PHINV(1 - PROB) ×
SQRT(K × (1 - P)) + K + P / 2) / P;
if X > K then X:= K;
PCUM:= 1 - NEGBIN(X - 1, K, P);
PX:= NEGBINPROB(X, K, P);
if PCUM < PROB then
begin for PX:= PX × (X - K) / (1 - P) / (X - 1)
while PCUM + PX < PROB do
begin X:= X - 1; PCUM:= PCUM + PX end
end else
begin for PCUM:= PCUM - PX
while PCUM > PROB do
begin PX:= PX × (1 - P) × X / (X - K + 1);
X:= X + 1
end; X:= X + 1
end
end;
NEGBININV:= X
end NEGBININV;
eop
code 41013;
real procedure POISSON(X, MU);
value X, MU; real X, MU;
begin integer IX;
real procedure KSI(K, L);
value K, L; real L; integer K;
begin real U, U2, W; W:= SQRT(L);
U:= 2 × (SQRT(K + 1) - W);
U2:= U × U; KSI:=
U + (U2 - 4) / 12 / W + (-U2 × U + 10 × U) / 72 /
L + (21 × U2 × U2 - 371 × U2 - 52) / 6480 / L / W
end;
IX:= ENTIER(X);
if IX < 0 then POISSON:= 0 else
if MU ≤ 0 then
POISSON:= STATAL3 ERROR(“POISSON”,2,MU)
else
if MU > 1000 then POISSON:= PHI(KSI(IX, MU))
else
begin integer I, MODUS; real MODUSPROB, PROB, CUM;
MODUS:= ENTIER(MU) + 1; if IX < MODUS then
begin PROB:= CUM:= POISSONPROB(IX, MU);
for I:= IX step -1 until 1 do
begin PROB:= PROB × I / MU;
CUM:= CUM + PROB
end
end else
begin MODUSPROB:= PROB:= CUM:=
POISSONPROB(MODUS, MU);
for I:= MODUS step -1 until 1 do
begin PROB:= PROB × I / MU;
CUM:= CUM + PROB
end; PROB:= MODUSPROB;
for I:= MODUS + 1 step 1 until IX do
begin PROB:= PROB × MU / I;
CUM:= CUM + PROB
end
end;
POISSON:= CUM
end
end POISSON;
eop
code 41014;
real procedure POISSONINV(PROB, MU, LEFT);
value PROB, MU, LEFT;
real PROB, MU; Boolean LEFT;
begin integer X; real PX, PCUM;
if PROB ≤ 0 ∨ PROB ≥ 1 then
STATAL3 ERROR(“POISSONINV”, 1, PROB) else
if MU ≤ 0 then
STATAL3 ERROR(“POISSONINV”, 2, MU);
if LEFT then
begin X:= (PHINV(PROB) / 2 + SQRT(MU)) ⭡ 2 - 1;
if X < 0 then X:= 0;
if PROB < EXP(-MU) then POISSONINV:= -1 else
begin PX:= POISSONPROB(X, MU);
PCUM:= POISSON(X, MU);
if PCUM > PROB then
begin for PCUM:= PCUM - PX
while PCUM > PROB do
begin PX:= PX × X / MU; X:= X - 1 end;
X:=X - 1
end else
begin for PX:= PX × MU / (X + 1)
while PCUM + PX < PROB do
begin X:= X + 1; PCUM:= PCUM + PX end
end;
POISSONINV:= X
end
end else
begin X:= (PHINV(1 - PROB) / 2 + SQRT(MU)) ⭡ 2 + 1;
if X < 0 then X:= 0;
PCUM:= 1 - POISSON(X - 1, MU);
PX:= POISSONPROB(X, MU);
if PCUM < PROB then
begin for PX:= PX × X / MU
while PCUM + PX < PROB do
begin X:= X - 1; PCUM:= PCUM + PX end
end else
begin for PCUM:= PCUM - PX
while PCUM > PROB do
begin PX:= PX × MU / (X + 1); X:= X + 1 end;
X:= X + 1
end;
POISSONINV:= X
end
end POISSONINV;
eop
code 41020;
real procedure WILCOX(X,M,N);
value X,M,N; real X,M,N;
begin
integer procedure MIN(A,B);
value A,B; integer A,B;
MIN := if A > B then B else A;
integer procedure MAX(A,B);
value A,B; integer A,B;
MAX := if A > B then A else B;
real WP1; Boolean X EVEN, RIGHT;
integer MN;
if M < 0 ∨ ENTIER(M) < M then
STATAL3 ERROR(“WILCOX”, 2, M);
if N < 0 ∨ ENTIER(N) < N then
STATAL3 ERROR(“WILCOX”, 3, N);
MN := M × N;
X := ENTIER(X/2);
if X < MN/2 then RIGHT := false else
begin RIGHT := true; X := MN-X-1 end;
X EVEN := ENTIER(X/2) × 2 = X;
M := MIN(M,N); N := MN/M;
if X < 0 then WP1 := 0 else
if M = 1 then WP1 := (X+1) / (N+1) else
if M = 2 then WP1 := (if X EVEN then
(X+2)×(X+2) / (2×(N+1)×(N+2)) else
(X+1)×(X+3) / (2×(N+1)×(N+2)) ) else
if 2×X = MN - 1 then WP1 := .5 else
if MN > 400 then
begin integer NOEM, N2, N3, N4, M2, M3, M4;
real F0Y, F3Y, F5Y, F7Y, T3, T5, T7, Y, Y2;
Y := (2×X + 1 - MN) / SQRT(MN × (M + N + 1) / 3);
Y2 := Y × Y;
NOEF := 10 × MN × (M + N + 1); F0Y := PHIDENS(Y);
N2 := N × N; N3 := N2 × N; N4 := N2 × N2;
M2 := M × M; M3 := M2 × M; M4 := M2 × M2;
F3Y := Y × (3 - Y2);
F5Y := Y × (-15 + Y2 × (10 - Y2));
F7Y := Y × (105 - Y2 × (105 - Y2 × (21 - Y2)));
T3 := (M2 + N2 + MN + M + N) / NOEM / 2;
T5 := ( 2 × (M4 + N4) + 4 +
(M3 × N + N3 × M + M3 + N3) +
6 × M2 × N2 + 7 × MN × (M + N) + M2 + N2 +
2 × MN - M - N) / (NOEM × NOEM × 2.1);
T7 := (M2 + N2 + MN + M + N) ⭡ 2
/ (NOEM × NOEM × 8);
WP1 := MAX(PHI(Y) - F0Y ×
(T3 × F3Y - T5 × F5Y - T7 × F7Y), 0);
end else
begin integer I,J,W,UP,UP1,UP2; real H1,H2;
if N × (X+1) ≤ 12000 then
begin M := N; N := MN / M end;
begin real array WP[0:X, 1:M];
UP2 := MIN(M, ENTIER((MN-X-1)/(N-1)));
for I := MAX(2, -ENTIER(X/2-M)) step 1
until UP2 do
begin UP:= X-(M-I)×2; UP1 := MIN(UP,I-1);
H1 := 1/(I+1);
for W:= MAX(0, X-(M-I)×N) step 1
until UP1 do
WP[W,I] := H1 × (W+1);
end;
for J := 2 step 1 until N do
begin UP:= MIN(X-(M-1)×J , J-1);
H2 := 1/(J+1);
for W:= MAX(0, X-(M-2)×N-J) step 1
until UP do
WP[W,1] := H2 × (W+1);
UP2 := (if J×M < X+1 then
ENTIER((MN-X-1)/(N-J)) else M);
for I := MAX(2,-ENTIER(X/J-M)) step 1
until UP2 do
begin UP:= X - (M-I)×J;
H1:= J/(I+J); H2 := I/(I+J);
UP1 := MIN(UP,J-1);
for W := MAX(0, X-(M-I)×N) step 1
until UP1 do
WP[W,I] := WP[W,I]×H1;
UP1 := MIN(UP,I×J-I-1);
for W := MAX(J, X-(M-I)×N) step 1
until UP1 do
WP[W,I] := WP[W,I]×H1 + WP[W-J,I-1]×H2;
UP1 := MIN(UP,I×J-1);
for W:= MAX(I×J-I, X-(M-I)×N) step 1
until UP1 do
WP[W,I] := H1 + WP[W-J,I-1]×H2;
end
end;
WP1 := WP[X,M];
end;
end;
WILCOX := if RIGHT then 1-WP1 else WP1;
end WILCOXCDF;
eop
code 41021;
real procedure WILCOXINV(PROB, M, N, LEFT);
value PROB, M, N, LEFT; real PROB, M, N; Boolean LEFT;
begin
integer X, W, WI, MN; real Z; Boolean MN EVEN;
MN := M × N;
MN EVEN := ENTIER(MN/2) × 2 = MN;
X := if MN EVEN then MN/2 - 1 else MN/2 - 1.5;
PROB := PROB + 10-13×(1-PROB);
if PROB ≤ 0 ∨ PROB ≥ 1 then
STATAL3 ERROR(“WILCOXINV”, 1, PROB) else
if M < 0 ∨ ENTIER(M) < M then
STATAL3 ERROR(“WILCOXINV”, 2, M) else
if N < 0 ∨ ENTIER(N) < N then
STATAL3 ERROR(“WILCOXINV”, 2, N) else
if MN = 0 then WI := -2 else
if PROB = .5 then WI := ENTIER((MN-1)/2) × 2 else
if M = 1 then WI := ENTIER(PROB×(N+1)) × 2 - 2
else
if N = 1 then WI := ENTIER(PROB×(M+1)) × 2 - 2
else if MN > 400 ∨ M=2 ∨ N=2 then
begin Z := PHINV(PROB) × SQRT(MN×(M+N+1)/3) + MN;
WI := W := ENTIER(Z/2) × 2;
if WI < 0 then WI := W := 0;
if WI > 2×MN then WI := W := 2×MN;
if WILCOX(W, M, N) ≤ PROB then
begin for W := W + 2
while WILCOX(W, M, N) ≤ PROB do WI := W
end else
begin for W:= W - 2
while WILCOX(W, M, N) > PROB do WI := W;
WI := WI - 2;
end;
end else
begin integer I,J,UP,UP1; real H1,H2;
Boolean RIGHT, EQUAL; real array WCMN[-1:X+2];
integer procedure MAX(A,B);
value A,B; integer A,B;
MAX := if A > B then A else B;
integer procedure MIN(A,B);
value A,B; integer A,B;
MIN := if A > B then B else A;
RIGHT := PROB > .5;
if RIGHT then PROB := 1 - PROB;
I := MAX(M,N); J := MIN(M,N);
if I × (X+1) > 12000 then
begin M := J; N := I end else
begin M := I; N := J end;
begin real array WP[0:X, 1:M];
for I := MAX(2,-ENTIER(X/2-M)) step 1
until M do
begin UP:= X-(M-I)×2; UP1 := MIN(UP,I);
H1 := 1/(I+1);
for W:= 0 step 1 until UP1 do
WP[W,I] := H1;
end;
for J := 2 step 1 until N do
begin UP:= MIN(X-(M-1)×J , J); H2 := 1/(J+1);
for W:= 0 step 1 until UP do
WP[W,1] := H2;
for I := MAX(2,-ENTIER(X/J-M)) step 1
until M do
begin UP:= X - (M-I)×J;
H1:= J/(I+J); H2 := I/(I+J);
UP1 := MIN(UP,J-1);
for W := 0 step 1 until UP1 do
WP[W,I] := WP[W,I]×H1;
UP1 := MIN(UP,I×J-I);
for W := J step 1 until UP1 do
WP[W,I] := WP[W,I]×H1 + WP[W-J,I-1]×H2;
UP1 := MIN(UP,I×J);
for W := I×J-I+1 step 1 until UP1
do WP[W,I] := WP[W-J,I-1]×H2;
end
end;
WCMN[-1] := 0;
WCMN[0] := WP[0,M];
for W := 1 step 1 until X do
WCMN[W] := WCMN[W-1] + WP[W,M];
if MN EVEN then
begin WCMN[X+1] := 1 - WCMN[X];
WCMN[X+2] := 1 - WCMN[X-1];
end
else
begin WCMN[X+1] := .5;
WCMN[X+2] := 1 - WCMN[X];
end;
end;
WI := PHINV(PROB) × SQRT(MN×(M+N+1)/3) + MN;
W := ENTIER(WI/2);
if W < 0 then WI := W := 0 else WI:= W;
if WCMN[W] ≤ PROB then
begin for W := W + 1
while WCMN[W] ≤ PROB do WI := W
end else
begin for W := W - 1
while WCMN[W] > PROB do WI := W;
WI := WI - 1;
end;
EQUAL := WCMN[WI] = PROB;
if RIGHT then
begin if EQUAL then WI := 2 × (MN - WI - 1)
else WI := 2 × (MN - WI - 2)
end
else WI := 2 × WI;
end;
WILCOXINV := if LEFT then WI else 2×MN - WI;
end WILCOXINV;
eop
code 41022;
real procedure WILCOXPROB(X,M,N);
value X,M,N; real X,M,N;
begin
integer procedure MIN(A,B);
value A,B; integer A,B;
MIN := if A> B then B else A;
integer procedure MAX(A,B);
value A,B; integer A,B;
MAX := if A> B then A else B;
real WP1;
integer MN;
if M < 0 ∨ ENTIER(M) < M then
STATAL3 ERROR(“WILCOXPROB”, 2, M);
if N < 0 ∨ ENTIER(N) < N then
STATAL3 ERROR(“WILCOXPROB”, 3, N);
MN := M × N;
X := MIN(X/2, MN - X/2);
M := MIN(M,N); N := MN/M;
if ENTIER(X) < X then WP1 := 0 else
if X < 0 then WP1 := 0 else
if MN = 0 then WP1 := 1 else
if M = 1 then WP1 := 1/(N+1) else
if M = 2 then
WP1 := ENTIER(X/2+1) / ((N+1)×(N+2)/2) else
if MN > 400 then
WP1 := WILCOX(2×X,M,N) - WILCOX(2×X-2,M,N) else
begin integer I,J,W,UP,UP1,UP2; real H1,H2;
if N × (X+1) ≤ 12000 then
begin M := N; N := MN/M end;
begin real array WP[0:X, 1:M];
UP2 := MIN(M, ENTIER((MN-X)/(N-1)));
for I := MAX(2,-ENTIER(X/2-M)) step 1
until UP2 do
begin UP:= X-(M-I)×2; UP1 := MIN(UP,I);
H1 := 1/(I+1);
for W:= MAX(0, X-(M-I)×N) step 1
until UP1 do
WP[W,I] := H1;
end;
for J := 2 step 1 until N do
begin UP:= MIN(X-(M-1)×J , J); H2 := 1/(J+1);
for W:= MAX(0, X-(M-2)×N-J) step 1
until UP do
WP[W,1] := H2;
UP2 := if J×M<X+1 then
ENTIER((MN-X)/(N-J)) else M;
for I := MAX(2,-ENTIER(X/2-M)) step 1
until UP2 do
begin UP:= X - (M-I)×J;
H1:= J/(I+J); H2 := I/(I+J);
UP1 := MIN(UP,J-1);
for W := MAX(0, X-(M-I)×N) step 1
until UP1 do
WP[W,I] := WP[W,1]×H1;
UP1 := MIN(UP,I×J-I);
for W := MAX(J, X-(M-I)×N) step 1
until UP1 do
WP[W,I] := WP[W,I]×H1 + WP[W-J,I-1]×H2;
UP1 := MIN(UP,I×J);
for W:=MAX(I×J-I+1,X-(M-I)×N) step 1
until UP1 do
WP[W,I] := WP[W-J,I-1]×H2;
end
end;
WP1 := WP[X,M];
end;
end;
WILCOXPROB := WP1
end WILCOXPROB;
eop
code 41023;
real procedure RUN(X, M, N);
value X, M, N; real X, M, N;
begin
real P, PCUM; integer IX, I, K, UP;
if M < 0 ∨ ENTIER(M) < M then
STATAL3 ERROR(“RUN”, 2, M) else
if N < 0 ∨ ENTIER(N) < N then
STATAL3 ERROR(“RUN”, 3, N) else
if M > N then begin K:= M; M:= N; N:= K end;
if M = 0 then RUN:= (if X < 1 then 0 else 1)
else if X < 2 then RUN:= 0 else
if M = N ∧ X ≥ M × 2 then RUN:= 1 else
if X > M × 2 then RUN:= 1 else
begin IX:= ENTIER(X); if IX ÷ 2 × 2 < IX then
begin PCUM:= P:= RUNPROB(IX, M, N);
K:= (IX - 1) / 2;
P:= P×K×2 / (M + N - K × 2); PCUM:= PCUM + P
end else
begin K:= IX / 2;
P:= PCUM:= RUNPROB(IX, M, N)
end;
for I:= K - 1 step -1 until 1 do
begin P:= P × (M + N - I × 2) × I / (N - I) /
(M - I) / 2; PCUM:= PCUM + P;
P:= P×2×I / (M + N - 2 × I); PCUM:= PCUM + P
end;
RUN:= PCUM
end
end RUN;
eop
code 41024;
real procedure RUNINV(PROB, M, N, LEFT);
value PROB, M, N, LEFT; real PROB, M, N; Boolean LEFT;
begin
integer H, R, MN; real P, PCUM;
if PROB ≤ 0 ∨ PROB ≥ 1 then
STATAL3 ERROR(“RUNINV”, 1, PROB) else
if M < 0 ∨ ENTIER(M) < M then
STATAL3 ERROR(“RUNINV”, 2, M) else
if N < 0 ∨ ENTIER(N) < N then
STATAL3 ERROR(“RUNINV”, 3, N);
MN:= M + N;
if M > N then begin H:= M; M:= N; N:= H end;
if M = 0 then
RUNINV:= (if LEFT then 0 else 2) else
if LEFT then
begin R:= PHINV(PROB) ×
SQRT((MN ⭡ 3 - MN) / (2 × M × N ×
(2 × M × N - MN))) + .5 + 2 × M × N / MN;
if R < 2 then R:= 2 else
if R > M × 2 then R:= M × 2;
if PROB < RUNPROB(2, M, N)
then RUNINV:= +1 else
begin PCUM:= RUN(R, M, N);
if PCUM ≤ PROB then
begin for P:= RUNPROB(R + 1, M, N)
while PCUM + P ≤ PROB do
begin R:= R + 1; PCUM:= PCUM + P end
end else
begin for P:= RUNPROB(R, M, N)
while PCUM - P > PROB do
begin R:= R - 1; PCUM:= PCUM - P end;
R:= R - 1
end; RUNINV:= R
end
end else
begin R:= PHINV(1 - PROB) ×
SQRT((MN ⭡ 3 - MN) / (2 × M × N ×
(2 × M × N - MN))) + 1.5 + 2 × M × N / MN;
if R < 2 then R:= 2 else
if R > M × 2 then R:= M × 2;
MN:= if M = N then 2 × M else 2 × M + 1;
if PROB < RUNPROB(MN, M, N)
then RUNINV:= MN + 1 else
begin PCUM:= 1 - RUN(R - 1, M, N);
if PCUM ≤ PROB then
begin for P:= RUNPROB(R - 1, M, N)
while PCUM + P ≤ PROB do
begin R:= R - 1; PCUM:= PCUM + P end
end else
begin for P:= RUNPROB(R, M, N)
while PCUM - P > PROB do
begin R:= R + 1; PCUM:= PCUM - P end;
R:= R + 1
end;
RUNINV:= R
end
end
end RUNINV;
eop
code 41025;
real procedure RUNPROB(X, M, N); value X, M, N;
real X, M, N;
begin integer K; real P; Boolean EVEN;
if M < 0 ∨ ENTIER(M) < M then
STATAL3 ERROR(“RUNPROB”, 2, M) else
if N < 0 ∨ ENTIER(N) < N then
STATAL3 ERROR(“RUNPROB”, 3, N);
if M > N then begin K:= N; N:= M; M:= K end;
EVEN:= ENTIER(X / 2) × 2 = X;
K:= if EVEN then X / 2 else (X + 1) / 2;
RUNPROB :=
if N = 0 then 0 else
if M = 0 then (if X = 1 then 1 else 0) else
if X < 2 ∨ X > 2 × M + 1 ∨ ENTIER(X) < X
then 0 else
if EVEN then
2 × M × N × EXP(2 × (LOGGAMMA(M) + LOGGAMMA(N) -
LOGGAMMA(K)) - LOGGAMMA(M - K + 1) - LOGGAMMA(N - K + 1)
- LOGGAMMA(M + N + 1)) else
(M + N - X + 1) × M × N / (K - 1) × EXP(2 × (LOGGAMMA(M)
+ LOGGAMMA(N) - LOGGAMMA(K - 1)) - LOGGAMMA(M + N + 1)
- LOGGAMMA(M - K + 2) - LOGGAMMA(N - K + 2))
end RUNPROB;
eop
code 41026;
real procedure KENDALL(X, N); value X, N; real X, N;
begin
integer I, G, IX; real P;
if N < 0 ∨ ENTIER(N) < N then
STATAL3 ERROR(“KENDALL”, 2, N);
G:= N × (N - 1) / 2; IX:= G + ENTIER(-(G - X) / 2) × 2;
if IX ≥ G then KENDALL:= 1 else
if IX < -G then KENDALL:= 0 else
if N > 9 then
KENDALL:= PHI(IX + 1 / SQRT(N × (N - 1) × (N+N+5) / 18))
else if IX > 0 then
begin P:= 0; for I:= G step -2 until IX + 2 do
P:= P + KENDALLPROB(I, N); KENDALL:= 1 - P
end else
begin P:= 0; for I:= -G step 2 until IX do
P:= P + KENDALLPROB(I, N); KENDALL:= P
end
end KENDALL;
eop
code 41027;
real procedure KENDALLINV(PROB, N, LEFT);
value PROB, N, LEFT; real PROB, N; Boolean LEFT;
begin integer S, G; real P, PCUM;
if PROB ≤ 0 ∨ PROB ≥ 1 then
STATAL3 ERROR(“KENDALLINV”, 1, PROB) else
if N < 0 ∨ ENTIER(N) < N then
STATAL3 ERROR(“KENDALLINV”, 2, N);
G:= N × (N - 1) / 2;
if N = 0 then
KENDALLINV:= (if LEFT then -1 else +1)
else if LEFT then
begin S:= PHINV(PROB) × SQRT(G × (N × 2 + 5) / 9);
S:= if ABS(S) > G then G × SIGN(S) else
G + ENTIER(-(G - S) / 2) × 2;
if PROB < KENDALLPROB(-G, N) then
KENDALLINV:= -G - 2
else
begin PCUM:= KENDALL(S, N);
if PCUM ≤ PROB then
begin for P:= KENDALLPROB(S + 2, N)
while PCUM + P ≤ PROB do
begin S:= S + 2; PCUM:= PCUM + P end
end else
begin for P:= KENDALLPROB(S, N)
while PCUM - P > PROB do
begin S:= S - 2; PCUM:= PCUM - P end;
S:=S-2
end; KENDALLINV:= S
end
end else
begin S:= PHINV(1 - PROB) × SQRT(G × (N × 2 + 5) / 9);
S:= if ABS(S) > G then G × SIGN(S) else
G - ENTIER((G - S) / 2) × 2;
if PROB < KENDALLPROB(G, N) then
KENDALLINV:= G + 2 else
begin PCUM:= 1 - KENDALL(S - 2, N);
if PCUM ≤ PROB then
begin for P:= KENDALLPROB(S - 2, N)
while PCUM + P ≤ PROB do
begin S:= S - 2; PCUM:= PCUM + P end
end else
begin for P:= KENDALLPROB(S, ND)
while PCUM - P > PROB do
begin S:= S + 2; PCUM:= PCUM - P end;
S:= S +2
end;
KENDALLINV:= S
end
end
end KENDALLINV;
eop
code 41028;
real procedure KENDALLPROB(X, N);
value X, N; real X, N;
begin integer G, IX;
real procedure PROB(S, N);
value S, N; integer S, N;
begin integer I; real P;
if N = 2 then PROB:= (if ABS(S) = 1 then .5
else 0) else
begin P:= 0;
for I:= 0 step 1 until N - 1 do
P:= P + PROB(S - N + 1 + I × 2, N - 1);
PROB:= P / N
end
end;
if N < 0 ∨ ENTIER(N) < N then
STATAL3 ERROR(“KENDALLPROB”, 2, N);
G:=N × (N - 1) / 2; IX:= ENTIER(X);
if N > 9 then
begin real S; S:= SQRT(N × (N-1) × (N+N+5) / 18);
KENDALLPROB:= PHI(IX + 1/S) - PHI(IX - 1/S)
end else
KENDALLPROB:= if IX < X ∨ ABS(IX) > G ∨
(G - IX) ÷ 2 × 2 < G - IX then 0 else PROB(IX, N)
end KENDALLPROB;
eop
code 41251;
real procedure BINPROB(X, N, P);
value X, N, P; real X, N, P;
BINPROB:= if N < 0 ∨ N > ENTIER(N)
then STATAL3 ERROR(“BINPROB”, 2, N)
else if P < 0 ∨ P > 1
then STATAL3 ERROR(“BINPROB”, 3, P)
else if X < 0 ∨ X > N ∨ X > ENTIER(X)
then 0 else if P = 0 ∨ N = 0
then (if X = 0 then 1 else 0)
else if P = 1
then (if X = N then 1 else 0)
else EXP(LOGGAMMA(N+1) - LOGGAMMA(X+1)
- LOGGAMMA(N-X+1) + X × LN(P) + (N-X) × LN(1-P));
eop
code 41252;
real procedure POISSONPROB(X, MU);
value X, MU; real X, MU;
POISSONPROB:= if MU ≤ 0
then STATAL3 ERROR(“POISSONPROB”, 2, MU)
else if X < 0 ∨ X > ENTIER(X) then 0
else EXP(-MU + X × LN(MU) - LOGGAMMA(X+1));
eop
code 41253;
real procedure HYPERGPROB(X, N, R, M);
value X, N, R, M; real X, N, R, M;
begin
integer procedure BINCOEF(N, K); value N, K;
integer N, K;
begin integer B, L, B1;
B1:= if K > N - K then N - K else K; B:= 1;
for L:= 1 step 1 until B1 do
B:= B × (N + 1 - L) ÷ L;
BINCOEF:= B
end;
if N < 0 ∨ N > M ∨ N > ENTIER(N)
then STATAL3 ERROR(“HYPERGPROB”, 2, N)
else if R < 0 ∨ R > M ∨ R > ENTIER(R)
then STATAL3 ERROR(“HYPERGPROB”, 3, R)
else if M > ENTIER(M)
then STATAL3 ERROR(“HYPERGPROB”, 4, M);
if X < 0 ∨ X < N + R - M ∨ X > N ∨ X > R
∨ X > ENTIER(X)
then HYPERGPROB:= 0 else
if N = 0 ∨ M = 0 then HYPERGPROB:=
(if X = 0 then 1 else 0) else
if N = M ∨ R = M then HYPERGPROB:=
(if X = M then 1 else 0) else
if M ≤ 51 then
HYPERGPROB:= (BINCOEF(N, X) × BINCOEF(M - N, R - X))
/ BINCOEF(M, R)
else
begin integer I; real PROB; PROB:= 0;
for I:= N, M-N, R, M-R do
PROB:= PROB + LOGGAMMA(I + 1);
for I:= N - X, X, M - N - R + X, R - X, M do
PROB:= PROB - LOGGAMMA(I + 1);
HYPERGPROB:= EXP(PROB)
end
end HYPERGPROB;
eop
code 41254;
real procedure NEGBINPROB(X, K, P);
value X, K, P; real X, K, P;
NEGBINPROB:= if K < 0 ∨ K > ENTIER(K)
then STATAL3 ERROR(“NEGBINPROB”, 2, K)
else if P ≤ 0 ∨ P > 1
then STATAL3 ERROR(“NEGBINPROB”, 3, P)
else if X < K ∨ X > ENTIER(X) then 0
else if P = 0 then 0
else if P = 1 ∨ K = 0
then (if X = K then 1 else 0)
else P × BINPROB(K - 1, X - 1, P);
eop
code 41255;
real procedure MULNOMPROB(XVEC, N, K, PVEC);
value N,K; real N,K; array XVEC, PVEC;
begin real XL, PL, LNPR, EPS, PSUM;
integer J, XSUM;
if N > ENTIER(N) ∨ N < 1 then
STATAL3 ERROR(“MULNOMPROB”,2,N) else
if K > ENTIER(K) ∨ K < 2 then
STATAL3 ERROR(“MULNOMPROB”,3,K) else
begin PSUM := 0; XSUM := 0; EPS := 10-14;
for J := 1 step 1 until K do
begin PL := PVEC[J];
if PL < 0 then
STATAL3 ERROR(“MULNOMPROB”,4,PL) else
PSUM := PSUM + PL
end;
if ABS(PSUM-1) > EPS then
STATAL3 ERROR(“MULNOMPROB”,4,PSUM) else
begin for J := 1 step 1 until K do
begin XL := XVEC[J];
if XL > ENTIER(XL) ∨ XL < 0 then
begin MULNOMPROB := 0; goto OUT end;
XSUM := XSUM + XL;
end;
if XSUM 10= N then
begin MULNOMPROB := 0; goto OUT end else
begin LNPR := LOGGAMMA(N+1);
for J := 1 step 1 until K do
if PVEC[J] = 0 then
begin if XVEC[J] ≠ 0 then
begin MULNOMPROB := 0;
goto OUT
end
end else
LNPR := LNPR - LOGGAMMA(XVEC[J]+1) +
XVEC[J] × LN(PVEC[J]);
MULNOMPROB := EXP(LNPR);
end;
end;
end; OUT:
end MULNOMPROB;
eop
code 41256;
real procedure MULHYPERGPROB(X,N,K,R);
value N,K; real N,K; array X,R;
begin integer I,J,L,SR,SX; real MHP,XJ,RJ;
if N < 1 ∨ N > ENTIER(N) then
STATAL3 ERROR(“MULHYPERGPROB”,2,N);
if K < 1 ∨ K > ENTIER(K) then
STATAL3 ERROR(“MULHYPERGPROB”,3,K);
SX := SR := 0;
MHP := 0;
for J := 1 step 1 until K do
begin RJ := R[J]; XJ := X[J];
if RJ < 0 ∨ RJ > ENTIER(RJ) then
STATAL3 ERROR(“MULHYPERGPROB”,4,RJ);
SR := SR + RJ; SX := SX + XJ;
if XJ > RJ ∨ XJ < 0 ∨ XJ > ENTIER(XJ) then
begin MULHYPERGPROB := 0; goto OUT end;
MHP := MHP + LOGGAMMA(RJ+1) - LOGGAMMA(XJ+1)
- LOGGAMMA(RJ-XJ+1);
end J ;
if SX ≠ N then
begin MULHYPERGPROB := 0; goto OUT end;
MHP := MHP + LOGGAMMA(N+1) + LOGGAMMA(SR-N+1)
- LOGGAMMA(SR+1);
MULHYPERGPROB := EXP(MHP);
OUT:
end MULHYPERGPROB;
eop
code 41500;
real procedure PHI(X); value X; real X;
begin real ABSX, ERF, ERFC, C, P, Q;
X:= X × .70710 67811 8655; ABSX:= ABS(X);
if X > 5.5 then PHI:= 1 else if X < -5.5
then PHI:= 0 else if ABSX ≤ 0.5 then
begin C:= X × X;
P:= ((-0.35609 84370 181510-1 × C +
0.69963 83488 619110+1) × C + 0.21979 26161 829410+2)
× C + 0.24266 79552 305310+3;
Q:= ((C +
0.15082 79763 040810+2) × C +
0.91164 90540 451510+2) × C +
0.21505 88758 698610+3;
PHI:= .5 × X × P / Q + .5
end else
begin if ABSX < 4 then
begin C:= ABSX;
P:= ((((((-0.13686 48573 827210-6 × C +
0.56419 55174 789710+0) × C +
0.72117 58250 883110+1) × C +
0.43162 22722 205710+2) × C +
0.15298 92850 469410+3) × C +
0.33932 08167 343410+3) × C +
0.45191 89537 118710+3) × C +
0.30045 92610 201610+3;
Q:= ((((((C +
0.12782 72731 962910+2) × C +
0.77000 15293 523010+2) × C +
0.27758 54447 439910+3) × C +
0.63898 02644 656310+3) × C +
0.93135 40948 306110+3) × C +
0.79095 09253 279010+3) × C +
0.30045 92609 369810+3;
C:= P/Q
end else
begin C:= 1 / X / X;
P:= (((0.22319 24597 341910-1 × C +
0.27866 13086 096510-0) × C +
0.22695 65935 396910-0) × C +
0.49473 09106 232510-1) × C +
0.29961 07077 035410-2;
Q:= (((C +
0.19873 32018 171410+1) × C +
0.10516 75107 067910+1) × C +
0.19130 89261 078310+0) × C +
0.10620 92305 284710-1;
C:= (C × (-P) / Q + 0.56418 95835 4776) / ABSX
end;
PHI:= .5 + .5 × SIGN(X) × (1 - C × EXP(-X × X))
end
end PHI;
eop
code 41501;
real procedure PHINV(PROB); value PROB; real PROB;
begin real EPS;
real procedure INVERF(X); value X; real X;
begin real ABSX, P, BETAX;
real array A[0 : 23];
real procedure CHEPOLSER(N, X, A);
value N, X; integer N; real X; array A;
begin integer K; real H, R, S, TX;
TX:= X + X; R:= A[N];
H:= A[N - 1] + R × TX;
for K:= N - 2 step -1 until 1 do
begin S:= R; R:= H;
H:= A[K] + R × TX - S
end;
CHEPOLSER:= A[0]- R + H × X
end CHEPOLSER;
ABSX:= ABS(X);
if ABSX ≤ 0.8 then
begin
A[ 0]:= 0.99288 53766 1894;
A[ 1]:= 0.12046 75161 4310;
A[ 2]:= 0.01607 81993 4210;
A[ 3]:= 0.00268 67044 3716;
A[ 4]:= 0.00049 96347 3024;
A[ 5]:= 0.00009 88982 1860;
A[ 6]:= 0.00002 03918 1276;
A[ 7]:= 0.00000 43272 7162;
A[ 8]:= 0.00000 09380 8141;
A[ 9]:= 0.00000 02067 3472;
A[10]:= 0.00000 00461 5970;
A[11]:= 0.00000 00104 1668;
A[12]:= 0.00000 00023 7150;
A[13]:= 0.00000 00005 4393;
A[14]:= 0.00000 00001 2555;
A[15]:= 0.00000 00000 2914;
A[16]:= 0.00000 00000 0680;
A[17]:= 0.00000 00000 0159;
A[18]:= 0.00000 00000 0037;
A[19]:= 0.00000 00000 0009;
A[20]:= 0.00000 00000 0002;
A[21]:= 0.00000 00000 0001;
INVERF:= CHEPOLSER(21, X × X / 0.32 - 1, A) × X
end else
if 1 - ABSX ≥ 2510-4 then
begin
A[ 0]:= 0.91215 88034 1755;
A[ 1]:= -0.01626 62818 6766;
A[ 2]:= 0.00043 35564 7295;
A[ 3]:= 0.00021 44385 7007;
A[ 4]:= 0.00000 26257 5108;
A[ 5]:= -0.00000 30210 9105;
A[ 6]:= -0.00000 00124 0606;
A[ 7]:= 0.00000 00624 0661;
A[ 8]:= -0.00000 00005 4015;
A[ 9]:= -0.00000 00014 2321;
A[10]:= 0.00000 00000 3438;
A[11]:= 0.00000 00000 3358;
A[12]:= -0.00000 00000 0146;
A[13]:= -0.00000 00000 0081;
A[14]:= 0.00000 00000 0005;
A[15]:= 0.00000 00000 0002;
BETAX:= SQRT(- LN((1 + ABSX) × (1 - ABSX)));
P:= -1.54881 30423 7326 × BETAX +
2.56549 01231 4782;
P:= CHEPOLSER(15, P, A);
INVERF:= if X < 0 then - BETAX × P
else BETAX × P
end else
begin
A[ 0]:= 0.95667 97090 2049;
A[ 1]:= -0.02310 70043 0907;
A[ 2]:= -0.00437 42360 9751;
A[ 3]:= -0.00057 65034 2265;
A[ 4]:= -0.00001 09610 2231;
A[ 5]:= 0.00002 51085 4703;
A[ 6]:= 0.00001 05623 3607;
A[ 7]:= 0.00000 27544 1233;
A[ 8]:= 0.00000 04324 8450;
A[ 9]:= -0.00000 00205 3034;
A[10]:= -0.00000 00438 9154;
A[11]:= -0.00000 00176 8401;
A[12]:= -0.00000 00039 9129;
A[13]:= -0.00000 00001 8693;
A[14]:= 0.00000 00002 7292;
A[15]:= 0.00000 00001 3282;
A[16]:= 0.00000 00000 3183;
A[17]:= 0.00000 00000 0167;
A[18]:= -0.00000 00000 0204;
A[19]:= -0.00000 00000 0097;
A[20]:= -0.00000 00000 0022;
A[21]:= -0.00000 00000 0001;
A[22]:= 0.00000 00000 0001;
A[23]:= 0.00000 00000 0001;
BETAX:= SQRT(- LN((1 + ABSX) × (1 - ABSX)));
P:= -0.55945 76313 29832 × BETAX +
2.28791 57162 6336;
P:= CHEPOLSER(23, P, A);
INVERF:= if X < 0 then - BETAX × P
else BETAX × P
end
end INVERSE ERROR FUNCTION;
EPS:= 10-14;
if PROB < EPS ∨ 1 - PROB < EPS then
STATAL3 ERROR(“PHINV”, 1, PROB);
PHINV:= INVERF(2 × PROB - 1) × 1.41421356237310
end PHINV;
eop
code 41502;
real procedure NORMAL(X, MU, SIGMA);
value X, MU, SIGMA; real X, MU, SIGMA;
NORMAL:= if SIGMA ≤ 0
then STATAL3 ERROR(“NORMAL”, 3, SIGMA)
else PHI((X - MU) / SIGMA);
eop
code 41503;
real procedure NORMALINV(PROB, MU, SIGMA);
value PROB, MU, SIGMA; real PROB, MU, SIGMA;
NORMALINV:= if SIGMA ≤ 0
then STATAL3 ERROR(“NORMALINV”, 3, SIGMA)
else if PROB < 10-14 ∨ PROB > 1 - 10-14
then STATAL3 ERROR(“NORMALINV”, 1, PROB)
else MU + PHINV(PROB) × SIGMA;
eop
code 41506;
real procedure CHISQ(X, DF);
value X, DF; real X, DF;
CHISQ:= if DF ≤ 0 then
STATAL3 ERROR(“CHISQ”, 2, DF)
else if X ≤ 0 then 0 else
GAMMA(X, DF / 2, 2);
eop
code 41507;
real procedure CHISQINV(PROB, DF);
value PROB, DF; real PROB, DF;
if PROB < 10-10 ∨ PROB > 1 - 10-10
then STATAL3 ERROR(“CHISQINV”, 1, PROB)
else if DF ≤ 0
then STATAL3 ERROR(“CHISQINV”, 2, DF)
else
begin real X;
X:= PHINV(PROB) × SQRT(2 × DF) + DF;
CHISQINV:= INVERSE(X, CHISQ(X, DF), PROB, 10-10)
end CHISQINV;
eop
code 41509;
real procedure NCCHISQ(X, DF, DELTA);
value X, DF, DELTA; real X, DF, DELTA;
begin real FACTOR1, FACTOR2, PROB, SUM, TERM;
integer M;
if DF < 1 ∨ DF > ENTIER(DF) then
STATAL3ERROR(“NCCHISQ”, 2, DF);
if DELTA < 0 then
STATAL3ERROR(“NCCHISQ”, 3, DELTA);
if X ≤ 0 then NCCHISQ:= 0 else
if DELTA = 0 then NCCHISQ:= CHISQ(X, DF) else
begin PROB:= CHISQ(X, DF); X:= X / 2; DF:= DF / 2;
DELTA:= DELTA / 2; FACTOR1:= EXP(-DELTA);
FACTOR2:= EXP(DF × LN(X) - X - LOGGAMMA(DF + 1));
TERM:= SUM:= PROB × FACTOR1; M:= 0;
for M:= M + 1
while ¬( TERM < 10-9 ∧ M > DELTA ) do
begin FACTOR1:= FACTOR1 × DELTA / M;
PROB:= PROB - FACTOR2;
FACTOR2:= FACTOR2 × X / (DF + M);
TERM:= PROB × FACTOR1; SUM:= SUM + TERM
end;
NCCHISQ:= SUM
end;
end NCCHISQ;
eop
code 41513;
real procedure GAMMA(X, ALPHA, SCALE);
value X, ALPHA, SCALE; real X, ALPHA, SCALE;
begin integer DELTA, UPP;
real BETA, START, SUM, TERM;
real procedure INCGAM(X, A, EPS);
value X, A, EPS; real X, A, EPS;
begin real C0, C1, C2, D0, D1, D2, X2, AX, P,
Q, R, S, R1, R2, SCF; integer N;
S:= EXP(-X + A × LN(X)); SCF:= 10+300;
if X ≤ 1 then
begin X2:= X × X; AX:= A × X;
D0:= 1; P:= A; C0:= S;
D1:= (A + 1) × (A + 2 - X);
C1:= (D1 + AX + 2 × X) × S;
R2:= C1 / D1;
for N:= 1, N + 1
while ABS((R2 - R1) / R2) > EPS do
begin P:= P + 2;
Q:= (P + 1) × (P × (P + 2) - AX);
R:= N × (N + A) × (P + 2) × X2;
C2:= (Q × C1 + R × C0) / P;
D2:= (Q × D1 + R × D0) / P;
R1:= R2; R2:= C2 / D2;
C0:= C1; C1:= C2; D0:= D1; D1:= D2;
if ABS(C1) > SCF ∨ ABS(D1) > SCF then
begin C0:= C0 / SCF; C1:= C1 / SCF;
D0:= D0 / SCF; D1:= D1 / SCF
end
end; INCGAM:= R2 / A / EXP(LOGGAMMA(A))
end else
begin C0:= A × S; C1:= (1 + X) × C0;
Q:= X + 2 - A;
D0:= X; D1:= X × Q; R2:= C1 / D1;
for N:= 1, N + 1
while ABS((R2 - R1) / R2) > EPS do
begin Q:= Q + 2; R:=N × (N + 1 - A);
C2:= Q × C1 - R × C0; D2:= Q × D1 - R × D0;
R1:= R2; R2:= C2 / D2;
C0:= C1; C1:= C2; D0:= D1; D1:= D2;
if ABS(C1) > SCF ∨ ABS(D1) > SCF then
begin C0:= C0 / SCF; C1:= C1 / SCF;
D0:= D0 / SCF; D1:= D1 / SCF
end
end; INCGAM:= 1 - R2 / A / EXP(LOGGAMMA(A))
end
end INCGAM;
if ALPHA ≤ 0 then
STATAL3 ERROR(“GAMMA”, 2, ALPHA)
else if SCALE ≤ 0 then
STATAL3 ERROR(“GAMMA”, 3, SCALE)
else if X ≤ 0 then GAMMA:= 0 else
if ALPHA ≥ 500 then
begin
GAMMA:= PHI(((X / SCALE / ALPHA) ⭡ .33333333353333
- 1 + 1 / (9 × ALPHA)) × 3 × SQRT(ALPHA))
end else
begin X:= X / SCALE; BETA:= ALPHA - ENTIER(ALPHA) + 1;
START:= if X ≥ 40 then 1 else
INCGAM(X, BETA, 10-12);
if ALPHA < 1 then
GAMMA:= START + EXP(-X + ALPHA × LN(X)
- LOGGAMMA(ALPHA + 1))
else if ALPHA < 2 then GAMMA:= START
else if X > 700 then GAMMA:= 1
else
begin UPP:= ENTIER(ALPHA) - 2; SUM:= TERM:=
EXP(-X + (ALPHA - 1) × LN(X) - LOGGAMMA(ALPHA));
for DELTA:= 1 step 1 until UPP do
begin TERM:= TERM × (ALPHA - DELTA) / X;
SUM:= SUM + TERM
end;
GAMMA:= START - SUM
end
end
end GAMMA;
eop
code 41514;
real procedure GAMMAINV(PROB,ALPHA,SCALE);
value PROB,ALPHA,SCALE; real PROB,ALPHA,SCALE;
begin real X,TOL;
TOL:= 10-10;
if ALPHA ≤ 0 then
STATAL3 ERROR(“GAMMAINV”,2,ALPHA) else
if SCALE ≤ 0 then
STATAL3 ERROR(“GAMMAINV”,3,SCALE) else
if PROB ≤ TOL ∨ PROB ≥ 1 - TOL then
STATAL3 ERROR(“GAMMAINV”,1,PROB);
X:= ALPHA × SCALE;
GAMMAINV:= INVERSE(X,GAMMA(X,ALPHA,SCALE),PROB, TOL)
end GAMMAINV;
eop
code 41517;
real procedure BETA(X, ALPHA1, ALPHA2);
value X, ALPHA1, ALPHA2; real X, ALPHA1, ALPHA2;
BETA:= if ALPHA1 ≤ 0
then STATAL3 ERROR(“BETA”, 2, ALPHA1)
else if ALPHA2 ≤ 0
then STATAL3 ERROR(“BETA”, 3, ALPHA2)
else if X ≤ 0 then 0
else if X ≥ 1 then 1
else INCOMPLETE BETA(X, ALPHA1, ALPHA2, 10-12);
eop
code 41518;
real procedure BETAINV(PROB,ALPHA1,ALPHAZ2);
value PROB,ALPHA1,ALPHA2; real PROB,ALPHA1,ALPHA2;
begin real X,Y,TOL;
comment DEFINE ACCURACY;
TOL:= 10-10;
comment TEST FOR ADMISSIBILITY OF PARAMETERS;
if ALPHA1 ≤ 0 then
STATAL3 ERROR(“BETAINV”,2,ALPHA1)
else
if ALPHA2 ≤ 0 then
STATAL3 ERROR(“BETAINV”,3,ALPHA2)
else
if PROB ≤ TOL ∨ PROB ≥ 1 - TOL then
STATAL3 ERROR(“BETAINV”,1,PROB);
X:= 0;
BETAINV:= INVERSE(X,BETA(X,ALPHA1,ALPHA2),PROB,TOL)
end BETAINV;
eop
code 41521;
real procedure FISHER(X, DF1, DF2); value X, DF1, DF2;
real X, DF1, DF2;
begin real IB;
if DF1 ≤ 0 then
FISHER:= STATAL3 ERROR(“FISHER”,2,DF1)
else
if DF2 ≤ 0 then
FISHER:= STATAL3 ERROR(“FISHER”,3,DF2)
else
if X ≤ 0 then FISHER:= 0
else
begin IB:= INCOMPLETE BETA(DF2/(DF2 + DF1 × X)
, DF2/2,DF1/2,10-12);
if IB < 0 then IB:= 0
else if IB > 1 then IB:= 1;
FISHER:= 1 - IB
end
end FISHER;
eop
code 41522;
real procedure FISHERINV(PROB, DF1, DF2);
value PROB, DF1, DF2; real PROB, DF1, DF2;
begin
if PROB < 10-10 ∨ PROB > 1 - 10-10
then STATAL3 ERROR(“FISHERINV”, 1, PROB)
else if DF1 ≤ 0
then STATAL3 ERROR(“FISHERINV”, 2, DF1)
else if DF2 ≤ 0
then STATAL3 ERROR(“FISHERINV”, 3, DF2)
else
begin real X;
X:= if PROB ≤ .5 then .5 else
if DF2 ≤ 4 then 1 else
DF2 / (DF2 - 2) + PHINV(PROB) ×
SQRT(2 × DF2 × DF2 × (DF1 + DF2 - 2) /
(DF1 × (DF2 - 4) × (DF2- 2) × (DF2 - 2)));
FISHERINV:=
INVERSE(X, FISHER(X, DF1, DF2), PROB, 10-10)
end
end FISHERINV;
eop
code 41525;
real procedure NCFISHER(X,DF1,DF2,DELTA);
value X,DF1,DF2,DELTA; real X,DF1,DF2,DELTA;
begin integer J; real XX,FAKTOR1,FAKTOR2,EPS,SUM;
if DF1 ≤ 0 then STATAL3 ERROR(“NCFISHER”,2,DF1)
else
if DF2 ≤ 0 then STATAL3 ERROR(“NCFISHER”,3,DF2)
else
if DELTA < 0 then
STATAL3 ERROR(“NCFISHER”,4,DELTA)
else if X ≤ 0 then NCFISHER:= 0 else
begin XX:= (DF1 × X) / (DF1 × X + DF2); EPS:= 10-12;
DF1:= DF1 / 2; DF2:= DF2 / 2; DELTA:= DELTA / 2;
FAKTOR1:= 1;
FAKTOR2:= SUM:= INCOMPLETE BETA(XX,DF1,DF2,EPS);
if DELTA = 0 then goto UIT;
J:= 0; for J:= J + 1 while FAKTOR2 > EPS do
begin FAKTOR1:= FAKTOR1 × DELTA / J;
FAKTOR2:= FAKTOR1 ×
INCOMPLETE BETA(XX,DF1 + J,DF2,EPS);
SUM:= SUM + FAKTOR2
end;
UIT: NCFISHER:= EXP(-DELTA) × SUM
end
end NCFISHER;
eop
code 41530;
real procedure STUDENT(X, DF); value X, DF;
real X, DF;
begin realIB;
if DF ≤ 0 then
STUDENT:= STATAL3 ERROR(“STUDENT”,2,DF)
else
begin IB:=
INCOMPLETE BETA(DF/(DF + X × X),DF/2,0.5,10-12);
if IB < 0 then IB:= 0
else if IB > 2 then IB:= 2;
STUDENT:= if X < 0 then IB / 2 else 1 - IB / 2
end
end STUDENT;
eop
code 41531;
real procedure STUDENTINV(PROB, DF);
value PROB, DF; real PROB, DF;
begin
if PROB < 10-10 ∨ PROB > 1 - 10-10
then STATAL3 ERROR(“STUDENTINV”, 1, PROB)
else if DF ≤ 0
then STATAL3 ERROR(“STUDENTINV”, 2, DF)
else
begin real X, U, U2;
U:= PHINV(PROB); U2:= U × U;
X:= U × (1 + (U2 + 1) / 4 / DF +
(3 + U2 × (U2 × 5 + 16)) / 96 / DF / DF);
STUDENTINV:= INVERSE(X, STUDENT(X, DF), PROB, 10-10)
end
end STUDENTINV;
eop
code 41533;
real procedure NCSTUDENT(X, DF, DELTA);
value X, DF, DELTA; real X, DF, DELTA;
begin real A, B, A2, WB, D2, TOL, TOLI, H, HELP, RESULT;
Boolean DFEVEN;
real procedure INTEGRATE(Y0, Y4, F0, F2, F4);
value Y0, Y4, F0, F2, F4; real Y0, Y4, F0, F2, F4;
begin real F1, F3, Y2, TEE, Y;
Y2:=(Y0 + Y4)/2;
Y :=(Y0 + Y2)/2; F1:= EXP(H×(1 + Y×Y))/(1 + Y×Y);
Y :=(Y2 + Y4)/2; F3:= EXP(H×(1 + Y×Y))/(1 + Y×Y);
TEE:=6×F2 - 4×(F1 + F3) + F0 + F4;
INTEGRATE:=if ABS(TEE) < TOLI
then (Y4 - Y0)×(4×(F1 + F3) + 2×F2 +
F0 + F4 - TEE/15)
else INTEGRATE(Y0, Y2, F0, F1, F2) +
INTEGRATE(Y2, Y4, F2, F3, F4);
end INTEGRATE;
real procedure SUMMATION OF FACTORS M;
begin integer I;
real MSUM, COEF, MIMIN2, MIMIN1, MI;
Boolean ADD;
MSUM:= 0;
if DF > 1 then
begin
MIMIN2:= A×WB × EXP(H) × PHI(HELP×WB) ×
.3989422804;
if DFEVEN then MSUM:= MSUM + MIMIN2;
if DF > 2 then
begin COEF:= 1;
MIMIN1:= B×(HELP×MIMIN2 +
A×.1591549431×EXP(-.5×D2));
if ¬ DFEVEN then MSUM:= MSUM + MIMIN1;
ADD:= DFEVEN;
for I:= 2 step 1 until DF - 2 do
begin MI:= (I - 1)/I×B×
(COEF×HELP×MIMIN1 + MIMIN2);
if ADD then MSUM:= MSUM + MI;
ADD:=¬ ADD; COEF:= 1/(I - 1)/COEF;
MIMIN2:= MIMIN1; MIMIN1:= MI;
end I;
end DF>2;
end DF>1;
SUMMATION OF FACTORS M:= MSUM;
end SUMMATION OF FACTORS M;
procedure INITIALISATION;
begin TOL:= 10-8;
if DF < 1 ∨ ENTIER(DF) ≠ DF then
STATAL3 ERROR(“NCSTUDENT”, 2, DF);
DFEVEN:= ENTIER(DF/2) = DF/2;
A:= X/SQRT(DF); A2:= A×A; D2:= DELTA×DELTA;
HELP:= DELTA×A;
B:=DF/(DF + X×X); WB:= SQRT(B); H:=-D2×B×.5;
if ABS(A) > TOL then TOLI:= 180 × TOL / ABS(A);
end INITIALISATION;
INITIALISATION;
RESULT :=
if DFEVEN then
PHI(-DELTA) + SUMMATION OF FACTORS M × 2.5066282746
else
PHI(-DELTA×WB) + SUMMATION OF FACTORS M × 2 +
(if ABS(A) ≤ TOL then 0 else
.31830 98862 × INTEGRATE(0, A, EXP(H),
EXP(H×(1 + A2/4))/(1 + A2/4),
EXP(H×(1 + A2))/(1 + A2)) / 12);
NCSTUDENT:=
if TOL ≤ RESULT ∧ RESULT ≤ 1 - TOL then RESULT
else
if ABS(RESULT) < TOL then 0 else
if ABS(RESULT - 1) < TOL then 1 else
STATAL3 ERROR(“NCSTUDENT”, 0, RESULT);
end NCSTUDENT;
eop
code 41534;
real procedure NCSTUDENTINVC(PROB, DF, DELTA);
value PROB, DF, DELTA; real PROB, DF, DELTA;
begin real X, Y, TOL;
procedure CORNISH FISHER EXPANSION;
begin real UA, UA2, UA3, UA4, UA5;
integer DF4, DFDF;
UA:=PHINV(PROB);
UA2:=UA×UA; UA3:=UA2×UA; UA4:=UA2×UA2;
UA5:=UA4×UA; DF4:=DF×4; DFDF:=DF×DF;
X:=-UA/DFDF/32;
X:=X×DELTA - (UA2 - 1)/DFDF/24;
X:=X×DELTA + UA/DF4 + (UA3 + UA×4)/DFDF/16;
X:=X×DELTA + 1 + (UA2×2 + 1)/DF4 +
(UA4×4 + UA2×12 + 1)/DFDF/32;
X:=X×DELTA + UA + (UA3 + UA)/DF4 + (UA5×5 + UA3×16
+ UA×3)/DFDF/96;
end INITIAL APPROXIMATION BY CORNISH-FISHER METHOD;
TOL:= 10-7;
if PROB < TOL ∨ PROB > 1 - TOL then
STATAL3 ERROR(“NCSTUDENTINV”, 1, PROB);
if DF < 1 ∨ ENTIER(DF) ≠ DF then
STATAL3 ERROR(“NCSTUDENTINV”, 2, DF);
CORNISH FISHER EXPANSION;
NCSTUDENTINV:=
INVERSE(X, NCSTUDENT(X, DF, DELTA), PROB, TOL)
end NCSTUDENTINV;
eop
code 41539;
real procedure LOGNORMAL(X, MU, SIGMA);
value X, MU, SIGMA; real X, MU, SIGMA;
LOGNORMAL:= if SIGMA ≤ 0
then STATAL3 ERROR(“LOGNORMAL”, 3, SIGMA)
else if X ≤ 0 then 0
else PHI((LN(X) - MU) / SIGMA);
eop
code 41540;
real procedure LOGNORMALINV(PROB, MU, SIGMA);
value PROB, MU, SIGMA; real PROB, MU, SIGMA;
LOGNORMALINV:= if SIGMA ≤ 0
then STATAL3 ERROR(“LOGNORMALINV”, 3, SIGMA)
else if PROB < 10-14 ∨ PROB > 1 - 10-14
then STATAL3 ERROR(“LOGNORMALINV”, 1, PROB)
else EXP(PHINV(PROB) × SIGMA + MU);
eop
code 41541;
real procedure CAUCHY(X, LOC, SCALE);
value X, LOC, SCALE; real X, LOC, SCALE;
CAUCHY:= if SCALE ≤ 0
then STATAL3 ERROR(“CAUCHY”, 3, SCALE)
else ARCTAN((X - LOC) / SCALE) × .31830988618379 + .5;
eop
code 41542;
real procedure CAUCHYINV (PROB, LOC, SCALE);
value PROB, LOC, SCALE; real PROB, LOC, SCALE;
begin real ARG;
ARG:= 3.1415 92653 5898 × PROB;
CAUCHYINV:= if PROB ≤ 0 ∨ PROB ≥ 1 then
STATAL3 ERROR (“CAUCHYINV”, 1, PROB)
else if SCALE ≤ 0 then
STATAL3 ERROR (“CAUCHYINV”, 3, SCALE)
else
-SCALE × COS (ARG) / SIN (ARG) + LOC
end CAUCHYINV;
eop
code 41545;
real procedure WEIBULL(X, LOC, SCALE, ALPHA);
value X, LOC, SCALE, ALPHA; real X, LOC, SCALE, ALPHA;
WEIBULL:= if SCALE ≤ 0
then STATAL3 ERROR(“WEIBULL”, 3, SCALE)
else if ALPHA ≤ 0
then STATAL3 ERROR(“WEIBULL”, 4, ALPHA)
else if X ≤ LOC then 0
else 1 - EXP(-((X - LOC) / SCALE) ⭡ ALPHA);
eop
code 41546;
real procedure WEIBULLINV(PROB,LOC,SCALE,ALPHA);
value PROB,LOC,SCALE,ALPHA; real PROB,LOC,SCALE,ALPHA;
begin
if PROB ≤ 10-10 ∨ PROB ≥ 1 - 10-10 then
STATAL3 ERROR(“WEIBULLINV”,1,PROB) else
if SCALE ≤ 0 then
STATAL3 ERROR(“WEIBULLINV”,3,SCALE)
else
if ALPHA ≤ 0 then
STATAL3 ERROR(“WEIBULLINV”,4,ALPHA);
WEIBULLINV:= LOC + SCALE × (-LN(1 - PROB)) ⭡ (1 / ALPHA)
end WEIBULLINV;
eop
code 41550;
real procedure LOGISTIC(X, MU, SIGMA);
value X, MU, SIGMA; real X, MU, SIGMA;
LOGISTIC:= if SIGMA ≤ 0
then STATAL3 ERROR(“LOGISTIC”, 3, SIGMA)
else 1 / (1 + EXP(-(X - MU) / SIGMA));
eop
code 41551;
real procedure LOGISTICINV (PROB, MU, SIGMA);
value PROB, MU, SIGMA; real PROB, MU, SIGMA;
LOGISTICINV:= if SIGMA ≤ 0
then STATAL3 ERROR(“LOGISTICINV”, 3, SIGMA)
else if PROB ≤ 0 ∨ PROB ≥ 1
then STATAL3 ERROR(“LOGISTICINV”, 1, PROB)
else - SIGMA × LN((1 - PROB) / PROB) + MU;
eop
code 41556;
real procedure KOLSMIR(D, XSIZE, YSIZE, EPS);
value D, XSIZE, YSIZE, EPS; real D, XSIZE, YSIZE, EPS;
begin integer I, KGV, M, N;
integer procedure GGD(M, N);
value M, N; integer M, N;
GGD:= if N = 0 then M else GGD(N, M - M ÷ N × N);
procedure APPROX;
begin integer K; real SUM, X, TERM, THETA;
SUM:= .5; THETA:= (1 + (M / KGV) ⭡ 1.2) / (M + N);
X:= (I / KGV + THETA) ⭡ 2 × 2 / (1 / M + 1 / N);
for K:= 1, K + 2 while TERM > EPS do
begin TERM:= EXP(-X × K × K);
SUM:= SUM - TERM × (1 - EXP(-X × (2 × K + 1)))
end;
KOLSMIR:= 2 × SUM
end;
procedure EXACT;
begin integer array LOW[0:N]; array H[0:M];
integer DMN, MN1, X, Y, UPP; real SUM, BINCOEF;
BINCOEF:= 1; LOW[0]:= 0; MN1:= M + N + 1;
DMN:= I × M × N / KGV;
for X:= 1 step 1 until N do
begin integer T, TN;
T:= M × X - DMN; TN:= T ÷ N;
LOW[X]:= if T ≤ 0 then 0 else
if TN = T / N then TN else TN + 1;
BINCOEF:= BINCOEF × (MN1 - X) / X;
if BINCOEF > 10318 then
begin EPS:= 10-4;goto L end
end;
H[0]:= 1;
for Y:= 1 step 1 until M do H[Y]:= 0;
for X:= 0 step 1 until N do
begin Y:= LOW[X]; SUM:= H[Y];
UPP:= M - LOW[N - X];
for Y:= Y + 1 step 1 until UPP do
H[Y]:= SUM:= SUM + H[Y]
end;
KOLSMIR:= SUM / BINCOEF
end;
if XSIZE ≤ 0 ∨ XSIZE > ENTIER(XSIZE) then
STATAL3 ERROR(“KOLSMIR”, 2, XSIZE);
if YSIZE ≤ 0 ∨ YSIZE > ENTIER(YSIZE) then
STATAL3 ERROR(“KOLSMIR”, 3, YSIZE);
if XSIZE < YSIZE then
begin N:= XSIZE; M:= YSIZE end
else begin M:= XSIZE; N:= YSIZE end;
if EPS < 0 ∨ EPS > 10-2 then EPS:= 10-3;
if XSIZE < YSIZE then
begin N:= XSIZE; M:= YSIZE end
else begin M:= XSIZE; N:= YSIZE end;
KGV:= M × N / GGD(M, N);
I:= ENTIER((1 + 10-14) × D × KGV);
if EPS ≥ 10-3 then L: APPROX else EXACT;
end KOLSMIR;
eop
code 41558;
real procedure BIVANORM(H, K, RHO); value H, K, RHO;
real H, K, RHO;
begin real B;
real procedure V(H, K, EPS); value H, K, EPS;
real H, K, EPS;
if H = 0 ∨ K = 0 then V:= 0 else
if ABS(H) < ABS(K) then
V:= (PHI(H) - .5) × (PHI(K) - .5) - V(K, H, EPS)
else
if ABS(K)> 8 then
V:=.15915 49430 9189 × ARCTAN (K/H)
else
begin real M, L, L2, S, R, T, SS, TSN; integer N;
L:= K / H; M:= H × H / 2; L2:= L × L; R:= EXP(-M);
S:= 1 - R; T:= L; SS:= T × S;
for N:= 1, N + 1 while ABS(TSN) ≥ EPS do
begin R:= R × M / N; S:= S - R; T:= -T × L2;
TSN:=S × T / (2 × N + 1);
SS:= SS + TSN
end;
V:= SS × .15915 49430 9189
end V;
if H < -8 ∨ K < -8 then B:=0 else
if H > 8 ∧ K > 8 then B:=1 else
B:= if ABS(RHO) > 1 then
STATAL3 ERROR(“BIVANORM”, 3, RHO)
else if ABS(RHO) = 1 then
(if RHO = 1 then (if K ≤ H then PHI(K)
else PHI(H))
else (if H ≤ -K then 0
else PHI(K) - PHI(H)))
else V(H,(K - RHO × H)/ SQRT(1 - RHO ⭡ 2), 10-14)
+ V(K,(H - RHO × K)/ SQRT(1 - RHO ⭡ 2), 10-14)
+ .5 × (PHI(H) + PHI(K))
- .15915 49430 9189 × ARCCOS(RHO);
if B < 0 then BIVANORM:=0 else BIVANORM:=B;
end BIVANORM;
eop
code 41560;
real procedure STUDRANGE(Q,N,NU);
value Q,N,NU; real Q,N,NU;
begin real X, PI, LN4, LNSQRT2PI, LNSQRTPI4;
array E[1 : 3];
real procedure POWER(X)TO:(N); value X, N;
real X; integer N;
begin integer N2; real Y;
Y:= 1;
WHILE POS N:
if N ≤ 0 then goto END WHILE POS N;
N2:= N ÷ 2;
WHILE EVEN N:
if N2 × 2 ≠ N then goto END WHILE EVEN N;
N:= N2; X:= X × X; N2:= N ÷ 2;
goto WHILE EVEN N;
END WHILE EVEN N:
N:= N - 1; Y:= Y × X;
goto WHILE POS N;
END WHILE POS N:
POWER:= Y
end POWER;
real procedure RANGE(T,N); value T,N; real T,N;
begin real U; real array E[1:3];
E[1]:= E[2]:= 10-7;
RANGE:= N × QADRAT(U, -5, +5,
PHIDENS(U) × POWER(PHI(U + T) - PHI(U), N - 1), E);
end RANGE;
real procedure INTEGRAND(X); value X; real X;
begin real XQ;
XQ:= X / Q;
INTEGRAND:= EXP(NU × (LN4 + LN(XQ) - LNSQRT2PI
- XQ × XQ / 2)) × (1 - RANGE(X, N)) / X;
end INTEGRAND;
if N < 2 ∨ ENTIER(N) < N then
STATAL ERROR(“STUDRANGE”, 2, N);
if NU < 1 ∨ ENTIER(NU) < NU then
STATAL ERROR(“STUDRANGE”, 3, NU);
E[1]:= E[2]:= 10-6; PI:= ARCTAN(1) × 4;
LNSQRT2PI:= .5 × LN(2 × PI); LN4:= LN(4);
LNSQRTPI4:= .5 × LN(PI) - LN4;
STUDRANGE:= if Q ≤ 0 then 0 else
1 - 2 × EXP(NU × (LN(NU) / 2 + LNSQRTPI4)
- LOGGAMMA(NU / 2)) ×
QADRAT(X, 10-6, Q × 7, INTEGRAND(X), E);
end STUDRANGE;
eop
code 41561;
real procedure EXPON (X, LAMBDA);
value X, LAMBDA; real X, LAMBDA;
EXPON:= if LAMBDA ≤ 0
then STATAL3 ERROR(“EXPON”, 2, LAMBDA)
else if X ≤ 0 then 0
else 1 - EXP(-LAMBDA × X);
eop
code 41562;
real procedure EXPONINV(PROB, LAMBDA);
value PROB, LAMBDA; real PROB, LAMBDA;
EXPONINV:= if LAMBDA ≤ 0
then STATAL3 ERROR(“EXPONINV”,2,LAMBDA)
else if PROB ≤ 0 ∨ PROB ≥ 1
then STATAL3 ERROR(“EXPONINV”,1,PROB)
else - LN(1 - PROB) / LAMBDA;
eop
code 41563;
real procedure ERLANG(X, N, SCALE);
value X, N, SCALE; real X, N, SCALE;
ERLANG:= if SCALE ≤ 0
then STATAL3 ERROR(“ERLANG”, 3, SCALE)
else if N ≤ 0 ∨ ENTIER(N) < N
then STATAL3 ERROR(“ERLANG”, 2, N)
else GAMMA(X, N, SCALE);
eop
code 41564;
real procedure ERLANGINV(PROB,ALPHA,SCALE);
value PROB,ALPHA,SCALE; real PROB,ALPHA,SCALE;
if PROB ≤ 10-10 ∨ PROB ≥ 1 - 10-10 then
STATAL3 ERROR(“ERLANGINV”,1,PROB)
else if ALPHA ≤ 0 ∨ ENTIER(ALPHA) < ALPHA then
STATAL3 ERROR(“ERLANGINV”,2,ALPHA)
else if SCALE ≤ 0 then
STATAL3 ERROR(“ERLANGINV”,3,SCALE)
else ERLANGINV:= GAMMAINV(PROB,ALPHA, SCALE);
eop
code 41565;
real procedure LAPLACE(X, MU, SIGMA);
value X, MU, SIGMA; real X, MU, SIGMA;
begin
if SIGMA ≤ 0 then
STATAL3 ERROR(“LAPLACE”, 3, SIGMA);
X:= (X - MU) / SIGMA;
LAPLACE:= .5 × (1 + (1 - EXP(-ABS(X))) × SIGN(X))
end LAPLACE;
eop
code 41566;
real procedure LAPLACEINV(PROB,LOC,SCALE);
value PROB, LOC, SCALE; real PROB, LOC, SCALE;
begin
if SCALE ≤ 0 then
STATAL3 ERROR(“LAPLACEINV”,3,SCALE);
if PROB ≤ 0 ∨ PROB ≥ 1 then
STATAL3 ERROR(“LAPLACEINV”,1,PROB);
if PROB ≤ .5 then
LAPLACEINV := LOC + LN(2 × PROB) × SCALE
else
LAPLACEINV := LOC - LN(2 × (1 - PROB)) × SCALE;
end LAPLACEINV;
eop
code 41567;
real procedure UNIFORM(X, A, B);
value X, A, B; real X, A, B;
begin
if B ≤ A then
STATAL3ERROR(“UNIFORM”, 2, B);
UNIFORM:= if X ≤ A then 0 else
if X ≥ B then 1 else (X - A) / (B - A)
end UNIFORM;
eop
code 41568;
real procedure UNIFORMINV(PROB, A, B);
value PROB, A, B; real PROB, A, B;
begin
if B ≤ A then STATAL3ERROR(“UNIFORMINV”, 2, B);
if PROB ≤ 0 ∨ PROB ≥ 1 then
STATAL3ERROR(“UNIFORMINV”, 1, PROB);
UNIFORMINV:= (B - A) × PROB + A
end UNIFORMINV;
eop
code 41569;
real procedure BINORCOR(R,RHO,N);
value R, RHO, N; real R, RHO, N;
begin
real procedure SAMCORBIVNORDEN(R,RHO,N);
value R, RHO, N; real R, RHO, N;
begin
if ABS(R) ≥ 1 then SAMCORBIVNORDEN := 0
else
begin
real W1, W3, Y1, Y2, Y3, Y4, N1, R2, RRHO,
R2RHO2, W2, PI, RHO2;
R2 := R × R;
RHO2 := RHO × RHO;
RRHO := R × RHO;
R2RHO2 := R2 × RHO2;
W1 := SQRT(1 - R2);
W2 := SQRT(1 - RHO2);
W3 := SQRT(1 - RHO2 × R2);
PI := ARCTAN(1) × 4;
N1 := N - 1;
if N < 15 then
begin
real array SB[3:N]; integer I;
SB[3] :=
(1-RHO2) / PI / W1 × (1 + RRHO ×
ARCCOS(-RRHO) / W3) / (1 - R2RHO2);
if N =3 then SAMCORBIVNORDEN := SB[3]
else
begin
SB[4] :=
(1 - RHO2) × W2 / PI × (3 × RRHO +
(1 + 2 × R2RHO2) × ARCCOS(-RRHO) / W3)
/ (1 - R2RHO2) / (1 - R2RHO2);
if N = 4 then SAMCORBIVNORDEN := SB[4]
else
begin for I := 5 step 1 until N do
SB[I] :=
(2 × I - 5) / (I - 3) × RRHO × W1 × W2
/ (1 - R2RHO2) × SB[I-1] +
(I - 3) / (I - 4) × (1 - RHO2) ×
(1 - R2) / (1 - R2RHO2) × SB[I-2];
SAMCORBIVNORDEN := SB[N];
end;
end;
end
else
begin
Y1 := (RRHO + 2) / 8;
Y2 := (3 × RRHO + 2) × (3 × RRHO + 2) / 128;
Y3 := (((15 × RRHO + 18) × RRHO - 4)
× RRHO - 8) × 5 / 1024;
Y4 := ((((3675 × RRHO + 4200) × RRHO - 2520)
× RRHO - 3360) × RRHO - 336) / 32768;
SAMCORBIVNORDEN :=
(N - 2) / SQRT(N - 1) × (1 - RHO2) × W2
× (W1 × W2 / (1 - RRHO)) ⭡ N1 ×
SQRT((1 - RRHO) / 2 / PI)
/ ((1 - R2) × W1 × (1 - RHO2) × W2)
× ((((Y4 / N1 + Y3) / N1 + Y2) / N1 + Y1)
/ N1 + 1);
end;
end;
end SAMCORBIVNORDEN;
if ABS(RHO) ≥ 1 then
STATAL3 ERROR(“BINORCOR”, 2, RHO)
else if N > ENTIER(N) ∨ N < 3 then
STATAL3 ERROR(“BINORCOR”, 3, N)
else if R ≤ -1 then BINORCOR := 0
else if R ≥ 1 then BINORCOR := 1
else if RHO = 0 then
BINORCOR :=
STUDENT(R × SQRT((N - 2) / (1 - R × R)), N - 2)
else
if N ≤ 500 then
begin real W1, W2, W3, PI, R2, RHO2, RHO3,
RHO4, RRHO, R2RHO2;
R2 := R × R;
RHO2 := RHO × RHO;
RHO3 := RHO2 × RHO;
RHO4 := RHO2 × RHO2;
RRHO := R × RHO;
R2RHO2 := R2 × RHO2;
W1 := SQRT(1 - R2);
W2 := SQRT(1 - RHO2);
W3 := SQRT(1 - RHO2 × R2);
PI := ARCTAN(1) × 4;
if N = 3 then
BINORCOR :=
(ARCCOS(-R) - RHO × W1 / W3 × ARCCOS(-RRHO)) / PI
else if N = 4 then
BINORCOR :=
W1 × W2 × SAMCORBIVNORDEN(R, RHO, 3) / RHO -
(W2 / RHO - ARCCOS(RHO)) / PI
else if N = 5 then
BINORCOR :=
W1 × W2 × SAMCORBIVNORDEN(R, RHO, 4) / 2 / RHO
- R × (1 - R2) / 2 × SAMCORBIVNORDENC(R, RHO, 3)
- W1 × (1 + RHO2) / (2 × PI × RHO) × ARCCOS(-RRHO)
/ W3 + ARCCOS(-R) / PI
else
if N = 6 then
BINORCOR :=
W2 × (1 - 4 × RHO2) / (3 × PI × RHO3)
+ ARCCOS(RHO) / PI
- (1 - RHO2) × W1 × W2 / 3 / RHO3 ×
SAMCORBIVNORDEN(R, RHO, 3)
+ (1 - RHO2) × R / 3 / RHO2 ×
SAMCORBIVNORDEN(R, RHO, 4)
+ W1 × W2 / 3 / RHO × SAMCORBIVNORDEN(R, RHO, 5)
else if N = 7 then
BINORCOR :=
ARCCOS(-R) / PI - (3 + 6 × RHO2 - RHO4) ×
ARCCOS(-RRHO) / W3 × W1 / (8 × PI × RHO)
-R × (1 - R2) × (4 - 3 × RHO2 + 3 × RHO4) / 8 /
RHO2 × SAMCORBIVNORDEN(R, RHO, 3)
- R2 × W1 × W2 × (2 - RHO2) / 8 / RHO ×
SAMCORBIVNORDEN(R, RHO, 4)
+ (1 - RHO2) × R / 4 / RHO2 ×
SAMCORBIVNORDEN(R, RHO, 5)
+ W1 × W2 / 4 / RHO × SAMCORBIVNORDEN(R, RHO, 6)
else if N = 8 then
BINORCOR := ARCCOS(RHO)
/ PI - W2 × (3 - 11 × RHO2 + 23 × RHO4)
/ 15 / PI / RHO4 / RHO
+ (W2 / RHO) ⭡ 5 × W1 / 5 ×
SAMCORBIVNORDEN(R, RHO, 3)
- R × (1-RHO2) × (1-RHO2) / 5 / RHO4 ×
SAMCORBIVNORDEN(R, RHO, 4)
+ (3 × R2 - 1) × (1 - RHO2) × W2 / W1 / 15 / RHO3
× SAMCORBIVNORDEN(R, RHO, 5)
+ (1 - RHO2) × R / 5 / RHO2 ×
SAMCORBIVNORDEN(R, RHO, 6)
+ W1 × W2 / 5 / RHO × SAMCORBIVNORDENC(R, RHO, 7)
else
begin real array E[1:3]; real X;
E[1] := E[2] := 10-6;
BINORCOR :=
STUDENT(-RHO × SQRT(N - 1) / W2, N - 1)
+ QADRAT(X, 0, R, SAMCORBIVNORDEN(X, RHO, N), E)
end;
end
else
begin real R2, RHO2, RHO3; integer N1, N2, N3;
N1 := N - 1; N2 := N1 × N1; N3 := N1 × N2;
R2 := R × R;
RHO2 := RHO × RHO;
RHO3 := RHO2 × RHO;
BINORCOR :=
if ABS(RHO) ≤ .7 then
PHI((R × SQRT((N - 2.5) / (1 - R2)) -
RHO × SQRT((N - 1.5) / (1 - RHO2)))
/ SQRT(1 + RHO2 / 2 / (1 - RHO2) + R2 /
2 / (1 - R2)))
else
PHI((.5 × LN((1 + R) × (1 - RHO) / (1 - R) /
(1 + RHO)) - RHO / 2/ N1 -
(5 × RHO + RHO3) / 8 / N2 - (11 × RHO +
(2 + RHO2) × 3 × RHO3) / 16 / N3)
/ SQRT(1 / N1 + (4 - RHO2) / 2 / N2 +
(22 - (6 - 3 × RHO2) × RHO2) / 6 / N3));
end;
end BINORCOR;
eop
code 41571;
real procedure EXTVAL(X, LOC, SCALE);
value X, LOC, SCALE; real X, LOC, SCALE;
begin
if SCALE ≤ 0 then
STATAL3 ERROR(“EXTVAL”, 3, SCALE);
X:= -(X - LOC) / SCALE;
EXTVAL:= if X > LN(-LN(MINREAL))
then 0
else EXP(-EXP(X));
end EXTVAL;
eop
code 41572;
real procedure EXTVALINV(PROB, LOC, SCALE);
value PROB, LOC, SCALE; real PROB, LOC, SCALE;
begin
if SCALE ≤ 0 then
STATAL3 ERROR(“EXTVALINV”, 3, SCALE);
if PROB ≤ 0 ∨ PROB ≥ 1 then
STATAL3 ERROR(“EXTVALINV”, 1, PROB);
EXTVALINV:= -SCALE × LN(-LN(PROB)) + LOC
end EXTVALINV;
eop
code 41751;
real procedure UNIFORMDENS(X, A, B);
value X, A, B; real X,A,B;
begin
if B ≤ A
then STATAL3ERROR(“UNIFORMDENS”, 2, B);
UNIFORMDENS:= if X ≤ A ∨ X > B
then 0 else 1 / (B - A)
end UNIFORMDENS;
eop
code 41752;
real procedure PHIDENS(X); value X;
PHIDENS:= .39894228040143 × EXP(- X × X / 2);
eop
code 41753;
real procedure NORMALDENS(X, MU, SIGMA);
value X, MU, SIGMA; real X, MU, SIGMA;
NORMALDENS:= if SIGMA ≤ 0
then
STATAL3 ERROR(“NORMALDENS”, 3, SIGMA)
else EXP(-(((X - MU) / SIGMA) ⭡ 2) / 2)
× .39894228040143 / SIGMA;
eop
code 41754;
real procedure LOGNORMALDENS(X, MU, SIGMA);
value X, MU, SIGMA; real X, MU, SIGMA;
LOGNORMALDENS:= if SIGMA ≤ 0
then
STATAL3 ERROR(“LOGNORMALDENS”, 3, SIGMA)
else if X ≤ 0 then 0
else
EXP(-(((LN(X) - MU) / SIGMA) ⭡ 2) / 2)
× .39894228040143 / X / SIGMA;
eop
code 41755;
real procedure EXPONDENS(X, LAMBDA); value X, LAMBDA;
real X, LAMBDA;
EXPONDENS:= if LAMBDA ≤ 0
then
STATAL3 ERROR(“EXPONDENS”, 2, LAMBDA)
else if X ≤ 0 then 0
else LAMBDA × EXP(- LAMBDA × X);
eop
code 41756;
real procedure GAMMADENS(X,ALPHA, SCALE);
value X,ALPHA,SCALE; real X,ALPHA,SCALE;
if X ≤ 0 then GAMMADENS:= 0 else
begin
if ALPHA ≤ 0
then STATAL3 ERROR(“GAMMADENS”,2,ALPHA)
else if SCALE ≤ 0
then STATAL3 ERROR(“GAMMADENS”, 35,SCALE);
GAMMADENS :=
EXP(- ALPHA × LN(SCALE) - LOGGAMMA(ALPHA) -
X / SCALE + (ALPHA - 1) × LN(X))
end GAMMADENS;
eop
code 41757;
real procedure ERLANGDENS(X,ALPHA,SCALE);
value X,ALPHA,SCALE; real X,ALPHA,SCALE;
if X < 0
then STATAL3 ERROR(“ERLANGDENS”,1,X) else
if ALPHA ≤ 0 ∨ ENTIER(ALPHA) < ALPHA then
STATAL3 ERROR(“ERLANGDENS”,2,ALPHA) else
if SCALE ≤ 0
then STATAL3 ERROR(“ERLANGDENS”,3,SCALE)
else ERLANGDENS:= if X = 0 then 0 else
EXP(-ALPHA × LN(SCALE) - LOGGAMMA(ALPHA) - X / SCALE +
(ALPHA - 1) × LN(X));
eop
code 41758;
real procedure CHISQDENS(X, DF);
value X, DF; real X, DF;
CHISQDENS:= if DF ≤ 0
then STATAL3ERROR(“CHISQDENS”, 2, DF)
else if X ≤ 0 then 0 else
EXP((DF / 2 - 1) × LN(X) - X / 2 -
DF × LN(2) / 2 - LOGGAMMA(DF / 2));
eop
code 41759;
real procedure WEIBULLDENS(X,LOC,SCALE,ALPHA);
value X,LOC,SCALE,ALPHA; real X,LOC,SCALE,ALPHA;
begin
if SCALE ≤ 0
then STATAL3 ERROR(“WEIBULLDENS”,3,SCALE) else
if ALPHA ≤ 0
then STATAL3 ERROR(“WEIBULLDENS”,4,ALPHA);
WEIBULLDENS:= if X ≤ LOC then 0 else
(ALPHA / SCALE) × EXP((ALPHA - 1) ×
LN((X - LOC) / SCALE) - ((X - LOC) / SCALE) ⭡ ALPHA)
end WEIBULLDENS;
eop
code 41760;
real procedure BETADENS(X,ALPHA1,ALPHA2);
value X,ALPHA1,ALPHA2; real X,ALPHA1,ALPHA2;
begin real BET;
if ALPHA1 ≤ 0
then STATAL3 ERROR(“BETADENS”,2,ALPHA1);
if ALPHA2 ≤ 0
then STATAL3 ERROR(“BETADENS”,3,ALPHA2);
if X ≤ 0 ∨ X ≥ 1 then BETADENS:= 0 else
begin BET:= EXP(LOGGAMMA(ALPHA1 + ALPHA2) -
LOGGAMMA(ALPHA1) - LOGGAMMA(ALPHA2));
BETADENS:= BET × EXP((ALPHA1 - 1) × LN(X) +
(ALPHA2 - 1 ) × LN(1 - X))
end;
end BETADENS;
eop
code 41761;
real procedure FISHERDENS(X, DF1, DF2);
value X, DF1, DF2;
real X, DF1, DF2;
FISHERDENS:= if DF1 ≤ 0
then STATAL3 ERROR(“FISHERDENS”, 2, DF1)
else
if DF2 ≤ 0
then STATAL3 ERROR(“FISHERDENS”, 3, DF2)
else
if X ≤ 0 then 0 else
EXP(LOGGAMMA((DF1 + DF2) / 2) -
LOGGAMMA(DF1 / 2) - LOGGAMMA(DF2 / 2) +
(DF1 × LN(DF1) + DF2 × LN(DF2)) / 2 +
(DF1 / 2 - 1) × LN(X) - (DF1 + DF2) / 2 ×
LN(DF1 × X + DF2));
eop
code 41762;
real procedure STUDENTDENS(X, DF);
value X, DF; real X, DF;
STUDENTDENS:= if DF ≤ 0
then STATAL3 ERROR(“STUDENTDENS”, 2, DF)
else
EXP(LOGGAMMA((DF + 1) / 2) - LOGGAMMA(DF / 2) -
(DF + 1) / 2 × LN(1 + X × X / DF) -
LN(DF) / 2 - .57236494299247);
eop
code 41763;
real procedure CAUCHYDENS(X,LOC,SCALE);
value X,LOC,SCALE; real X,LOC,SCALE;
begin real PI,Q;
if SCALE ≤ 0
then STATAL3 ERROR(“CAUCHYDENS”,3,SCALE);
PI:= 3.1415926535898;
Q:= (X - LOC) / SCALE;
CAUCHYDENS:= 1 / (PI × SCALE × (1 + Q × Q))
end CAUCHYDENS;
eop
code 41764;
real procedure LAPLACEDENS(X,LOC,SCALE);
value X, LOC, SCALE; real X, LOC, SCALE;
begin
if SCALE ≤ 0
then STATAL3 ERROR(“LAPLACEDENS”,3,SCALE);
X := (X - LOC) / SCALE;
LAPLACEDENS := .5 / SCALE × EXP(-ABS(X));
end LAPLACEDENS;
eop
code 41765;
real procedure LOGISTICDENS(X, MU, SIGMA);
value X, MU, SIGMA; real X, MU, SIGMA;
begin
if SIGMA ≤ 0
then STATAL3 ERROR(“LOGISTICDENS”, 3, SIGMA);
X:= EXP(-(X - MU) / SIGMA);
LOGISTICDENS:= X / ((1 + X) × (1 + X) × SIGMA)
end LOGISTICDENS;
eop
code 41766;
real procedure EXTVALDENS(X, LOC, SCALE);
value X, LOC, SCALE; real X, LOC, SCALE;
begin
if SCALE ≤ 0 then
STATAL3 ERROR(“EXTVALDENS”, 3, SCALE);
X:= (X - LOC) / SCALE;
EXTVALDENS:= EXP(-(X + EXP(-X))) / SCALE
end EXTVALDENS;
eop
code 49999;
real procedure ANDERSON DARLING(X, L, U, SORTED);
value L, U, SORTED;
array X; integer L, U; Boolean SORTED;
begin integer I, N;
real MU, SIGMA, XI, FACTOR, SUM, ESTIMATE;
N:= U - L +1;
if N ≤ 1 then
STATAL3 ERROR(“ANDERSON-DARLING”, 3, U);
if ¬ SORTED then VECQSORT(X, L, U);
if X[L] = X[U] then
STATAL3 ERROR(“ANDERSON-DARLING”, 1, X[L]);
comment FIRST THE ESTIMATION (IN THE USUAL WAY) OF
EXPECTATION AND STANDARD DEVIATION OF THE NORMAL
DISTRIBUTION. ;
MU:= SIGMA:= 0; ESTIMATE:= (X[L] + X[U]) / 2;
for I:= L step 1 until U do
begin XI:= X[I]:= X[I] - ESTIMATE;
MU:= MU + XI; SIGMA:= SIGMA + XI × XI;
end;
MU:= MU / N;
SIGMA:= SQRT((SIGMA - N × MU × MU) / (N - 1));
comment TRANSFORMATION OF THE OBSERVATIONS TO
UNIFORM(0, 1)-DISTRIBUTED QUANTITIES.;
for I:= L step 1 until U do
X[I]:=PHI((X[I] - MU) / SIGMA);
comment ANDERSON-DARLING TEST QUANTITY;
SUM:= 0; FACTOR:= -1;
for I:= L step 1 until U do
begin FACTOR:= FACTOR + 2;
SUM:= SUM +
FACTOR × (LN(X[I]) + LN(1 - X[L + U - I]))
end;
ANDERSON DARLING:=
(1 + 4 / N - 25 / N / N) × (-SUM / N - N)
end OF ANDERSON DARLING;
eop