Approximation of Generalized Chi-Squared Distribution

genx2cdf.jl

using LinearAlgebra
using Distributions

# Liu, Huan, Yongqiang Tang, and Hao Helen Zhang. 2009. “A New Chi-Square
# Approximation to the Distribution of Non-Negative Definite Quadratic Forms in
# Non-Central Normal Variables.” Computational Statistics & Data Analysis 53
# (4): 853–56.
function genx2cdf(A, μ, Σ, t)
    # k = 1
    AΣᵏ = A
    c₁ = 1dot(μ, AΣᵏ, μ)
    AΣᵏ = AΣᵏ*Σ
    c₁ += tr(AΣᵏ)
    # k = 2
    AΣᵏ = AΣᵏ*A
    c₂ = 2dot(μ, AΣᵏ, μ)
    AΣᵏ = AΣᵏ*Σ
    c₂ += tr(AΣᵏ)
    # k = 3
    AΣᵏ = AΣᵏ*A
    c₃ = 3dot(μ, AΣᵏ, μ)
    AΣᵏ = AΣᵏ*Σ
    c₃ += tr(AΣᵏ)
    # k = 4
    AΣᵏ = AΣᵏ*A
    c₄ = 4dot(μ, AΣᵏ, μ)
    AΣᵏ = AΣᵏ*Σ
    c₄ += tr(AΣᵏ)
    s₁ = c₃/c₂^(3/2)
    s₂ = c₄/c₂^2
    μ_Q = c₁
    σ_Q = √(2c₂)
    β₁ = √8*s₁
    β₂ = 12*s₂
    u = (t - μ_Q)/σ_Q
    if s₂^2 > s₂
        a = 1/(s₁ - √(s₁^2 - s₂))
    else
        a = 1/s₁
    end
    δ = s₁*a^3 - a^2
    l = a^2 - 2δ
    μᵪ = l + δ
    σᵪ = √2*a
    cdf(NoncentralChisq(l, δ), u*σᵪ + μᵪ)
end

Monte-Carlo estimation:

function monte_carlo_cdf(A, μ, Σ, t; n = 100_000)
    L = cholesky(Σ).L
    d = size(μ, 1)
    count = 0
    for _ in 1:n
        X = L*randn(d) .+ μ
        count += dot(X, A, X) < t
    end
    return count / n
end