04-11-2021, 08:49 AM
Goodmorning everyone. Can any of you point me to an algorithm written in HP PPL for the calculation of the "studentized range q distribution" CDF?
Thanks, Roberto.
Thanks, Roberto.
(04-13-2021 08:18 PM)compaqdrew Wrote: [ -> ]The trick will be the infinite integration limits.
Usually you can just use large/small values if the contribution outside the range will be small.
Studentized_Range_q_P_Value(q,k,ν):=
BEGIN
return 1 -
√(2/π) * k * (ν/2)^(ν/2) / Gamma(ν/2)
* ∫(s^(ν-1)*exp(-ν*s*s/2)
* ∫(exp(-z*z/2) * normald_cdf(z,z+q*s)^(k-1)
, z = -∞ .. +∞)
, s = 0. .. +∞);
END;
do
mathx = require'mathx'
local erf, tgamma = mathx.erf, mathx.tgamma
local exp, sqrt, pi, huge = math.exp, math.sqrt, math.pi, math.huge
local quad = require'quad'.quad
local function cdf(a, b) -- area from a to b
a = erf(0.7071067811865476*a)
b = erf(0.7071067811865476*b)
return (b-a) * 0.5
end
function Studentized_Range_q_P_Value(q,k,v)
local c = sqrt(2/pi) * k * (v/2)^(v/2) / tgamma(v/2)
local k1, v1 = k-1, v-1
local function fs(s)
local function fz(z) return exp(-z*z/2) * cdf(z,z+q*s)^k1 end
return s^v1 * exp(-v*s*s/2) * quad(fz, -huge, huge)
end
return 1 - c * quad(fs, 0, huge)
end
end
#cas
quad6(f, a, b, d, n, eps)
BEGIN
LOCAL k, c, t0, t1, t, p, s, err;
LOCAL fp, fm, u, r, x, y, mode;
fp := fm := k := c := 0.;
t0 := 7.38905609893065; // e^2
a := approx(a);
b := approx(b);
mode := (a-a!=0) + 2*(b-b!=0);
IF mode == 0 THEN // tanh-sinh
c,d,s := d, (b-a)/2, f((a+b)/2)
ELSE
IF mode == 3 THEN // sinh-sinh
s := f(c)
ELSE // exp-sinh
c, s := a, b; // c = finite edge
IF mode == 1 THEN c, s := b, a END
IF s < 0 THEN d:=-d END
s := f(c+d)
END
END
REPEAT
p, t, t1 := 0, sqrt(t0), t0;
IF k==0 THEN t1 := t END
t0 := t;
IF mode == 0 THEN // tanh-sinh
REPEAT
u := exp(1/t-t);
r := 2*u/(1+u);
x := d*r;
IF a+x != a THEN
y := f(a+x);
IF y-y==0 THEN fp:=y ELSE fp*=c END
END
IF b-x != b THEN
y := f(b-x);
IF y-y==0 THEN fm:=y ELSE fm*=c END
END
y := (t+1/t)*r/(1+u) * (fp+fm);
p, t := p+y, t*t1;
UNTIL abs(y) <= eps*abs(p);
ELSE
t *= 0.5;
REPEAT
r := exp(t-0.25/t);
u := r; fp := 0;
IF mode == 3 THEN // sinh-sinh
r := 0.5*r - 0.5/r;
u := 0.5*u + 0.5/u;
y := f(c - d*r);
IF y-y==0 THEN fp := y*u END
ELSE // exp-sinh
x := c + d/r;
IF x == c THEN break END
y := f(x);
IF y-y==0 THEN fp := y/u END
END
y := f(c + d*r);
IF y-y == 0 THEN fp += y*u END
fp *= (t+0.25/t);
p, t := p+fp, t*t1;
UNTIL abs(fp) <= eps*abs(p);
END
s, err, k := s+p, abs(s-p), k+1;
// print(k, d*ldexp(s,1-k), err/abs(s));
UNTIL err <= 10*eps*abs(s) OR k>n;
IF mode==1 OR (mode==3 and a>b) THEN d:=-d END
return d*ldexp(s,1-k), err/abs(s);
END;
quad(f,a,b) := quad6(f,a,b,1,6,1e-9);
#end
#cas // Studentized P_Value, from q,k,v
SRP_qkv(q,k,v):=
1 - √(2/π)*k*(v/2)^(v/2)/Gamma(v/2)
* quad(s -> s^(v-1)*exp(-v*s*s/2) *
int(exp(-z*z/2) * normald_cdf(z,z+q*s)^(k-1)
, z, -inf, inf)
, 0., inf) (1)
#end
(05-19-2021 02:39 PM)Albert Chan Wrote: [ -> ]I was unable to get double integrals using only quad. Any ideas ?
#cas // Studentized P_Value, from q,k,v
SRP_qkv(q,k,v)
BEGIN
LOCAL s, z;
1 - √(2/π)*k*(v/2)^(v/2)/Gamma(v/2)
* quad(s -> s^(v-1)*exp(-v*s*s/2) *
quad(z -> exp(-z*z/2) * normald_cdf(z,z+q*s)^(k-1)
, -inf, inf) (1)
, 0, inf) (1)
END
#end
#cas // Studentized P_Value, from q,k,v
SRP_qkv(q,k,v) :=
1 - √(2/π)*k*(v/2)^(v/2)/Gamma(v/2)
* quad(s -> s^(v-1)*exp(-v*s*s/2) *
quad(z -> exp(-z*z/2) * normald_cdf(z,z+q*s)^(k-1)
, -inf, inf) (1)
, 0, inf) (1)
#end
#cas // Studentized P_Value, from q,k,v
SRP_qkv(q,k,v) :=
1 - √(2/π)*k*(v/2)^(v/2)/Gamma(v/2)
* quad (s -> s^(v-1) * exp(-v*s*s/2) *
quad6(z -> exp(-z*z/2) * normald_cdf(z,z+q*s)^(k-1)
, -inf, inf, 4., 6, 1e-9) (1)
, 0, inf) (1)
#end
XCas moyal.cc Wrote:static gen normal_cdf(const gen & g,GIAC_CONTEXT){
return rdiv(erf(ratnormal(plus_sqrt2_2*g),contextptr)+plus_one,2,contextptr)
}
#cas // Studentized P_Value, from q,k,v
SRP_qkv(q,k,v) :=
1 - k*2^(2-k) * v^(v/2)/(Gamma(v/2)*√π)
* quad6(y -> y^(v-1) * exp(-v*y*y) *
quad6(x -> exp(-x*x) * (erf(x+q*y)-erf(x))^(k-1)
, -inf, inf, 2.8, 6, 1e-9) (1)
, 0.0, inf, 0.71, 6, 1e-9) (1)
#end
#PYTHON EXPORT quadr6
from sys import *
from math import *
f=argv[0]
a=float(argv[1])
b=float(argv[2])
d=float(argv[3])
n=float(argv[4])
eps=float(argv[5])
fp = fm = k = c = 0.
t0 = 7.38905609893065 # e^2
mode = (a-a != 0) + 2* (b-b != 0)
if mode == 0: # tanh-sinh
c, d, s = d, (b-a)/2, f((a+b)/2)
else:
if mode == 3: # sinh-sinh
s = f(c)
else: # exp-sinh
if mode == 2:
c = a
else:
c =b
if bool(a>b) == bool(mode==2):
d = −d
s = f(c+d)
while err >=10*eps*abs(s) or k<n:
p, t, t1 = 0, sqrt(t0), t0
if k==0:
t1 = t
t0 = t
if mode == 0:
while abs(y) >= eps*abs(p):
u = exp(1/t-t)
r = 2*u/(1+u)
x = d*r
if a+x != a:
y = f(a+x)
if y-y == 0:
fp = y
else:
fp *= c
if b-x != b:
y = f(b-x)
if y-y==0:
fm =y
else:
fm *=c
y = (t+1/t)*r/(1+u) * (fp+fm)
p, t = p+y, t*t1
else:
t *= 0.5
while abs(fp) >= eps*abs(p):
r = exp(t-0.25/t)
u = r
fp =0
if mode == 3:
r = 0.5*r - 0.5/r
u = 0.5*u + 0.5/u
y = f(c-d*r)
if y-y == 0:
fp = y*u
else:
x = c + d/r
if x == c:
break
y = f(x)
if y-y == 0:
fp = y/u
y = f(c+d*r)
if y-y == 0:
fp += y*u
fp *= (t+0.25/t)
p, t = p+fp, t*t1
s, err, k = s+p, abs(s-p), k+1
if mode == 1 or (mode==3 and a>b):
d = −d
ris=[d*ldexp(s,1-k), err/abs(s)]
print(ris)
#end
EXPORT Quadratura_Python(f,a,b,d,n,eps)
BEGIN
PYTHON(quadr6,f,a,b,d,n,eps);
END;