HP PPL and Statistical distributions

[Albert Chan]

(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.

This is indeed the case.

Recently, we discussed tanh-sinh, exp-sinh, sinh-sinh algorithm here
Combined as 1, we have quad (my Lua implementation here, emace67 double_exponential.py here).

First, I clean-up CAS code, and run in XCas (s lower limit changed to 0. produced some speed-up)

Code:

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;

Picked an example, from studentized-range-q-table
From first page, P-value = 0.10, k=8, df=12, we have q = 4.511

XCas> Studentized_Range_q_P_Value(4.511, 8, 12)

Evaluation time: 6.24, returned 0.100017379207

Now, rewrite code in lua, using quad to integrate.

Code:

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

lua> Studentized_Range_q_P_Value(4.511, 8, 12)
0.10001737920797527

We got the same answer ... But, LuaJIT take 0.003 second