HP PPL and Statistical distributions

[Albert Chan]

Now that we have it running, we may optimize a bit.

Inner integral, \(\int_{-∞}^∞\), work best if bulk of area inside \(\int_{-d}^d\) (quad defaulted d = 1)

Inner integrand is bell-shaped, with center shifted left, we like to expand d = 4 (*)

One way is simple substitution: \(\int_{-∞}^∞ f(x)\,dx = \int_{-∞}^∞ f(4u)(4\,du) \)
Now, LHS \(\int_{-4}^4\) = RHS \(\int_{-1}^1\)

This is the reason quad6() exposed other parameters. We can simply set d = 4.

Code:

#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

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

---

(*) without shift, quad(z->exp(-z^2/2), -inf, inf) work better with d=8
However, if shifted, it might get worse, as d=9 example shows.

lua> Q = require'quad'
lua> f = function(z) return exp(-z*z/2) end
lua> Q.quad(Q.count(f), -huge, huge, 1), Q.n
2.506628274630998       119
lua> Q.quad(Q.count(f), -huge, huge, 2), Q.n
2.506628274630838       101
lua> Q.quad(Q.count(f), -huge, huge, 4), Q.n
2.506628274630999       43
lua> Q.quad(Q.count(f), -huge, huge, 8), Q.n
2.506628274630997       33
lua> Q.quad(Q.count(f), -huge, huge, 9), Q.n
2.506628274630931       55