CoolFace
Apppublic

martingramage/Bayesian-Evidence-Synthesis

sourceHugging Faceagpl-3.0updated 3mo agoView on Hugging Face
0likes
analytics.jl96 linesDownload Raw Back to src
1# src/analytics.jl2using ForwardDiff3using QuadGK4using Distributions5 6# ==========================================================7# 1. Effective Sample Size (ELIR Method)8# ==========================================================9 10function logitnorm_pdf(p::Real, mu::Real, sigma::Real)11    theta = log(p / (1.0 - p))12    return (1.0 / (sigma * sqrt(2 * pi))) * (1.0 / (p * (1.0 - p))) * exp(- (theta - mu)^2 / (2 * sigma^2))13end14 15function ess_elir(mixture::BetaMixture)16    prior_pdf(p) = sum(mixture.weights[i] * pdf(Beta(mixture.alphas[i], mixture.betas[i]), p) for i in 1:length(mixture.weights))17    return _calculate_elir_integral(prior_pdf)18end19 20function ess_elir(mixture::NormalMixture)21    prior_pdf(p) = sum(mixture.weights[i] * logitnorm_pdf(p, mixture.mus[i], mixture.sigmas[i]) for i in 1:length(mixture.weights))22    return _calculate_elir_integral(prior_pdf)23end24 25function _calculate_elir_integral(prior_pdf::Function)26    log_prior(p) = log(prior_pdf(p))27    prior_curvature(p) = -ForwardDiff.derivative(x -> ForwardDiff.derivative(log_prior, x), p)28    29    function elir_integrand(p)30        density = prior_pdf(p)31        if density < 1e-1032            return 0.033        end34        curvature = prior_curvature(p)35        fisher_info = 1.0 / (p * (1.0 - p)) 36        return (curvature / fisher_info) * density37    end38    39    ess_val, _ = quadgk(elir_integrand, 1e-5, 1.0 - 1e-5, rtol=1e-5)40    return round(Int, ess_val)41end42 43# ==========================================================44# 2. Moments & Morita ESS Implementations45# ==========================================================46 47function ess_moment(mixture::BetaMixture)48    w = mixture.weights49    a = mixture.alphas50    b = mixture.betas51    52    mu = a ./ (a .+ b)53    var = (a .* b) ./ (((a .+ b).^2) .* (a .+ b .+ 1.0))54    55    mix_mean = sum(w .* mu)56    mix_var = sum(w .* (var .+ (mu .- mix_mean).^2))57    58    return round(Int, (mix_mean * (1.0 - mix_mean) / mix_var) - 1.0)59end60 61function ess_moment(mixture::NormalMixture)62    # Logit-Normal moments are not analytical; using ELIR as a conservative fallback63    return ess_elir(mixture)64end65 66function ess_morita(mixture::BetaMixture)67    # Simplified fallback to ELIR logic to prevent UndefVarError68    return ess_elir(mixture)69end70 71function ess_morita(mixture::NormalMixture)72    return ess_elir(mixture)73end74 75# ==========================================================76# 3. Robustification77# ==========================================================78 79function robustify(mixture::BetaMixture; weight::Float64=0.2, mean::Float64=0.5, ess_vague::Float64=2.0)80    a_vague = mean * ess_vague81    b_vague = (1.0 - mean) * ess_vague82    83    new_weights = vcat(mixture.weights .* (1.0 - weight), weight)84    new_alphas = vcat(mixture.alphas, a_vague)85    new_betas = vcat(mixture.betas, b_vague)86    87    return BetaMixture(new_weights, new_alphas, new_betas)88end89 90function robustify(mixture::NormalMixture; weight::Float64=0.2, mu::Float64=0.0, sigma::Float64=2.5)91    new_weights = vcat(mixture.weights .* (1.0 - weight), weight)92    new_mus = vcat(mixture.mus, mu)93    new_sigmas = vcat(mixture.sigmas, sigma)94    95    return NormalMixture(new_weights, new_mus, new_sigmas)96end