Tiny Julia script to verify the the validity of the Wilk's theorem for a particular test statistics

stat-proj-INSS18.jl

using StatsBase, Distributions, Optim, ProgressMeter, Plots

nexp = 100000     # number of simulated experiments
nobs_mean = 1000  # number of expected counts in experiment

function logL(μ, h::Histogram)  # compute -log(L(μ))
    logL = 0.0
    # loop over bins
    for b in 1:length(h.edges[1])-1
        low  = h.edges[1][b]
        up   = h.edges[1][b+1]
        # cumulative of nobs_mean*(0.1μ(1-x)^2 + x^2) -> b+s hypothesis
        I(x) = nobs_mean*((1+0.1μ)x^3 - 0.3μ*x^2 + 0.3μ*x)
        pred = I(up) - I(low)  # expected counts in bin
        if pred < 0
            pred = 0
        end
        logL -= log(pdf(Poisson(pred), h.weights[b]))
    end
    return logL
end

hq = fit(Histogram, [], 0:0.01:5, closed = :left)  # histogram for q_0 values

# loop over experiments
@showprogress 0.1 "Generating experiments..." for n in 1:nexp
    # create bkg-only data with x^2, x ∈ [0,1]
    h = fit(Histogram,
            [cbrt(rand()) for ev in 1:rand(Poisson(nobs_mean))],
            0:0.1:1,
            closed = :left)
    # find likelihood minimum
    res = optimize(μ -> logL(μ,h), -1, 1)
    if res.converged
        if res.minimizer < 0
            push!(hq, 0)
        else
            # fill with q_0
            push!(hq, 2(logL(0,h) - logL(res.minimizer,h)))
        end
    else
        println("WARNING: convergence not reached")
    end
end

# plot!
f(x) = 1.71120e+02*(1/sqrt(x))*exp(-x/2)
plot(hq.edges[1][2:end],
     [hq.weights,f],
     xlabel="q_0",
     ylabel="\frac{dN_{exp}}{dq_0}",
     legend=false,
     xlims=(-0.1,5),
     yscale=:log10)