Skip to content

Instantly share code, notes, and snippets.

@bicycle1885
Created November 6, 2021 19:24
Show Gist options
  • Save bicycle1885/9b470af6ec9856a4e0401bc94798f2a9 to your computer and use it in GitHub Desktop.
Save bicycle1885/9b470af6ec9856a4e0401bc94798f2a9 to your computer and use it in GitHub Desktop.
Approximation of Generalized Chi-Squared Distribution
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
@bicycle1885
Copy link
Author

bicycle1885 commented Nov 6, 2021

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

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment