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, ⏨-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 > ⏨5 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 ≔ - 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 + ⏨-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 ≔ AX(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 ≔ -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 - 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 ≔ ⏨-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 ≠ 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 1815⏨-1 × C +
0·69963 83488 6191⏨+1) × C + 0·21979 26161 8294⏨+2)× C + 0·24266 79552 3053⏨+3;
Q ≔ ((C +
0·15082 79763 0408⏨+2) × C +
0·91164 90540 4515⏨+2) × C +
0·21505 88758 6986⏨+3;
PHI ≔ ·5 × X × P / Q + ·5
end else
begin if ABSX < 4 then
begin C ≔ ABSX;
P ≔ ((((((-0·13686 48573 8272⏨-6 × C +
0·56419 55174 7897⏨+0) × C +
0·72117 58250 8831⏨+1) × C +
0·43162 22722 2057⏨+2) × C +
0·15298 92850 4694⏨+3) × C +
0·33932 08167 3434⏨+3) × C +
0·45191 89537 1187⏨+3) × C +
0·30045 92610 2016⏨+3;
Q ≔ ((((((C +
0·12782 72731 9629⏨+2) × C +
0·77000 15293 5230⏨+2) × C +
0·27758 54447 4399⏨+3) × C +
0·63898 02644 6563⏨+3) × C +
0·93135 40948 3061⏨+3) × C +
0·79095 09253 2790⏨+3) × C +
0·30045 92609 3698⏨+3;
C ≔ P/Q
end else
begin C ≔ 1 / X / X;
P ≔ (((0·22319 24597 3419⏨-1 × C +
0·27866 13086 0965⏨-0) × C +
0·22695 65935 3969⏨-0) × C +
0·49473 09106 2325⏨-1) × C +
0·29961 07077 0354⏨-2;
Q ≔ (((C +
0·19873 32018 1714⏨+1) × C +
0·10516 75107 0679⏨+1) × C +
0·19130 89261 0783⏨+0) × C +
0·10620 92305 2847⏨-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 ⩾ 25⏨-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 ≔ ⏨-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 < ⏨-14 ∨ PROB > 1 - ⏨-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 ∨ PROB > 1 - ⏨-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)
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 < ⏨-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 ≔ ⏨+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 + 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, ⏨-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;
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, ⏨-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;
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,⏨-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 ∨ PROB > 1 - ⏨-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)
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 ≔ ⏨-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,⏨-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 ∨ PROB > 1 - ⏨-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)
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 ≔ ×F2 - 4×(F1 + F3) + F0 + F4;
INTEGRATE ≔ ̲f 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 ≔ ⏨-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 ≔ F/(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 ≔ HINV(PROB);
UA2 ≔ A×UA; UA3 ≔ A2×UA; UA4 ≔ A2×UA2;
UA5 ≔ A4×UA; DF4 ≔ F×4; DFDF ≔ F×DF;
X ≔ UA/DFDF/32;
X ≔ ×DELTA - (UA2 - 1)/DFDF/24;
X ≔ ×DELTA + UA/DF4 + (UA3 + UA×4)/DFDF/16;
X ≔ ×DELTA + 1 + (UA2×2 + 1)/DF4 +
(UA4×4 + UA2×12 + 1)/DFDF/32;
X ≔ ×DELTA + UA + (UA3 + UA)/DF4 + (UA5×5 + UA3×16
+ UA×3)/DFDF/96;
end INITIAL APPROXIMATION BY CORNISH-FISHER METHOD;
TOL ≔ ⏨-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 < ⏨-14 ∨ PROB > 1 - ⏨-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 ∨ PROB ⩾ 1 - ⏨-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 > ⏨318 then
begin EPS ≔ ⏨-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 > ⏨-2 then EPS ≔ ⏨-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 + ⏨-14) × D × KGV);
if EPS ⩾ ⏨-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 ≔ × T / (2 × N + 1);
SS ≔ SS + TSN
end;
V ≔ SS × ·15915 49430 9189
end V;
if H < -8 ∨ K < -8 then B ≔ else
if H > 8 ∧ K > 8 then B ≔ 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), ⏨-14)
+ V(K,(H - RHO × K)/ SQRT(1 - RHO ⭡ 2), ⏨-14)
+ ·5 × (PHI(H) + PHI(K))
- ·15915 49430 9189 × ARCCOS(RHO);
if B < 0 then BIVANORM ≔ else BIVANORM ≔ ;
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] ≔ ⏨-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] ≔ ⏨-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, ⏨-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 ∨ PROB ⩾ 1 - ⏨-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] ≔ ⏨-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] ≔ HI((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