HP PPL and Statistical distributions

[Albert Chan]

A nice thing about infinity, you can compress it, shift it, and it is still infinity.

cdf function, Φ(z) = 1/2 * (1+erf(z/√2)).
We can substitute x=√2*z, and avoid calls to Φ(z), which internally called erf.

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)
}

The substitution also simplified exp(-z*z/2) to exp(-x*x)

Do the same way for outer integral, letting y = √2*s, and pull out the constants, we have:
(d adjusted to offset effect of scaling: 4/√2 ≈ 2.8, 1/√2 ≈ 0.71)

Code:

#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

Cas> SRP_qkv(4.511, 8, 12)       → 0.100017379208, finished in 0.085 second