analysis.jl

using Distributions

# Komolgorov - Smirnov 
# ( adapted from HypothesisTests.jl )
function distance_ks(x::Vector{<:Real}, d::UnivariateDistribution)
    n = length(x)
    cdfs = cdf.(d, sort(x))
    δp = maximum((1:n) / n - cdfs)
    δn = -minimum((0:n-1) / n - cdfs)
    δ = max(δn, δp)
    return δ
end

function exact_distance_ks(f::Function, d::Distribution, h = 1e-4)
    x = [-1:h:1 ... ]
    F = cdf.(d, x)
    y = f.(x)
    p = exp.(y); p = p / (2 * mean(p))
    Δx = x[2:end].-x[1:end-1]
    P = cumsum(p .* h)
    return maximum(abs.(P .- F))
end 


# Cramér - von Mises
function distance_cm(x::Vector{<:Real}, d::UnivariateDistribution)
    n = length(x)
    x = sort(x)
    cdfs = cdf.(d, x)
    ω² = mean( ((1:n) / n .+ .5 - cdfs ).^2 ) + 1/(12*n^2)
    return ω²
end

function exact_distance_cm(f::Function, d::Distribution, h = 1e-4)
    x = [-1:h:1 ... ]
    y = f.(x)
    F = cdf.(d, x)
    f = pdf.(d, x)
    p = exp.(y); p = p / (2 * mean(p))
    P = cumsum(p .* h)
    return 2 * mean( (P .- F).^2 .* p )   
end 

# Quantify convergence rate
function asymptotic_rate(x::Vector, y::Vector)
    log_x = log.(x)
    log_y = log.(y) 
    A = hcat( ones( length(x)), log_x )

    (a,b) = (A'A) \ (A' * log_y)
    return b
end