Probabilistic measures with Monte Carlo draws

[sethaxen]

Sometimes we don't have a density function for a predictive distribution but instead have IID draws from that distribution. e.g. in Bayesian inference, we often only have access to the posterior predictive distribution through Monte Carlo draws. It would be useful if this package supported statistical estimation of the measure in such cases.

Proposal

I propose the API for probabilistic measures be expanded so that the first argument to a probabilistic measure can be an array of predictions representing these draws instead of just a distribution, where possible (e.g. this would not be possible for LogScore). When the first argument is an array of draws, then weights must be an array with the same size.

Example: CRPS

For example, the continuously ranked probability score (CRPS) can be defined for a predictive distribution q as
CRPS(q, y) = E[|ŷ-y|] - E[|ŷ-ŷ'|], where ŷ' and ŷ are independent random variables drawn from q.
For many named distributions, this can be evaluated exactly. But it can also be estimated from draws.

The naive estimator of CRPS for an IID sample of size N is.
mean(abs.(ŷ .- y)) - mean(abs.(ŷ .- ŷ')) * (N // (2 * (N-1))). For continuous predictions, this can be reformulated using order statistics to reduce the cost to O(Nlog(N)) (see https://hal.science/hal-02976423/document)

[ablaom]

Without having given this much thought, is the existing API not sufficient? If you have a number of draws y1, ..., yN of some distribution d, then your proxy for a prediction of y, is the corresponding empirical distribution d = Σ δ_yi, no? You could represent this by a raw array of draws, but it would be more consistent to use something that at least implements Base.rand.

At the moment LearnAPI.jl is not specific about the interface each kind of target proxy should support but i'd suggest Sampleable at least require Base.rand.

[sethaxen]

I think using some representation of an empirical distribution makes sense and isn't too cumbersome if a fit method exists for that object. For existing Monte Carlo draws, I don't think it makes sense for the proxy to be a Sampleable, since we won't and shouldn't sample from our draws; we already have the draws. Something like Distributions.DiscreteNonParametric could work, but here we need to be careful, since the exact estimator of $$E[|X-X'|]$$ for Distributions.DiscreteNonParametric will be biased if we apply it to a Monte Carlo sample, and a correction is needed. (Distributions currently has no type for an empirical distribution JuliaStats/Distributions.jl#1337).

But going with Distributions.DiscreteNonParametric, perhaps this is sensible?

using Distributions, StatisticalMeasuresBase, ScientificTypesBase, LearnAPI, IrrationalConstants

_pos_or_one(x) = ifelse(x ≤ 0, oneunit(x), x)

# mean absolute difference from self
# use corrected=true when probabilities are empirical and false otherwise
function mean_abs_diff(ŷ::Distributions.DiscreteNonParametric; corrected::Bool=false)
    p = Distributions.probs(ŷ)
    x = Distributions.support(ŷ)
    psum = pi = first(p)
    EXX = -2 * first(x) * pi * max(0, 1 - (1 - corrected) * pi)
    for (xi, pi) in Iterators.drop(zip(x, p), 1)
        psum += pi
        if corrected
            EXX += 2 * xi * pi * (1 - 2 * (1 - psum) / _pos_or_one(1 - pi))
        else
            EXX += 2 * xi * pi * (1 - 2 * (1 - psum) - pi)
        end
    end
    return EXX
end
mean_abs_diff(ŷ::Distributions.Normal; kwargs...) = 2 * Distributions.scale(ŷ) / sqrtπ

# mean absolute deviation from y
function mean_abs_dev(ŷ::Distributions.DiscreteNonParametric, y)
    p = Distributions.probs(ŷ)
    x = Distributions.support(ŷ)
    return sum(zip(x, p)) do (xi, pi)
        return pi * abs(xi - y)
    end
end
function mean_abs_dev(ŷ::Distributions.Normal, y)
    μ = Distributions.mean(ŷ)
    v = Distributions.var(ŷ)
    return 2 * v * Distributions.pdf(ŷ, y) + (μ - y) * (1 - 2 * Distributions.cdf(ŷ, y))
end

struct CRPS end

const crps = CRPS()

StatisticalMeasuresBase.@trait(
    CRPS,
    orientation=StatisticalMeasuresBase.Score(),
    kind_of_proxy=LearnAPI.Distribution(),
    observation_scitype=ScientificTypesBase.Infinite,
)

function (::CRPS)(ŷ::Distributions.UnivariateDistribution, y)
    EXX = mean_abs_diff(ŷ; corrected=false)
    EXy = mean_abs_dev(ŷ, y)
    return EXX / 2 - EXy
end
function (::CRPS)(ŷ::Distributions.DiscreteNonParametric, y; empirical::Bool=false)
    EXX = mean_abs_diff(ŷ; corrected=empirical)
    EXy = mean_abs_dev(ŷ, y)
    return EXX / 2 - EXy
end

Example:

julia> ŷ = randn(1_000);

julia> y = randn()
-1.1524247322809593

julia> crps(Distributions.fit(DiscreteNonParametric, ŷ), y; empirical=true)
-0.722180048862493

julia> crps(Normal(), y)
-0.7118369119295415

[sethaxen]

@ablaom Does the above look like a sensible approach? Also, would a contribution of CRPS to this package be welcome?

[ablaom]

Apologies for not following through earlier. I appreciate the time you've put into this.

Without having reviewed this in detail, this looks like you are fitting CRPS into the same scheme as existing proper scoring rules, LogLoss, BrierLoss, etc. The wrinkle is that while those other rules support UnivariateFinite (also "discrete", but where the sample elements are labelled and do not have the empirical interpretation) here we want to support a discrete distribution where the sample space is really continuous but has discrete support.

I agree, it looks it might be hard to get Distributions.jl to add an empirical distribution, given the past effort, unless perhaps you want to do that yourself. I agree using DiscreteNonParametric makes sense, since we already have UnivariateFinite from CategoricalDistributions to handle the bona fide discrete case.

Can you say how much of this can be made generic, ie apply to other distributions other than "empirical" and normal? I mean, ideally, those methods mean_abs_diff and mean_abs_dev should live at Distributions.jl, right? (I haven't; checked). I remember that when we implemented BrierScore we had to add some methods to Distributions.jl, which was quite a project.

In any, case, if only Normal and DiscreteNonParametric were supported, this would be a valuable contrbution I would support. Lead time for review would be a week or so.

You can check out the BrierScore doc string to get an idea of what is required here.

[sethaxen]

here we want to support a discrete distribution where the sample space is really continuous but has discrete support

For this specific example, yes, though when CRPS is formulated as E[|ŷ-y|] - E[|ŷ-ŷ'|] it also makes sense for discrete numeric sample spaces. e.g. if one has a posterior distribution and a Poisson observational model, CRPS could be used to evaluate the model via the posterior predictive distribution. It never really makes sense for a categorical distribution though.

Can you say how much of this can be made generic, ie apply to other distributions other than "empirical" and normal? I mean, ideally, those methods mean_abs_diff and mean_abs_dev should live at Distributions.jl, right? (I haven't; checked). I remember that when we implemented BrierScore we had to add some methods to Distributions.jl, which was quite a project.

I'm not certain these methods should live in Distributions, since they're quite specific. That being said, mean_abs_dev can be implemented quite generically using functions from Distributions:

function mean_abs_dev(ŷ::Distributions.UnivariateDistribution, y)
    prob_lower = cdf(ŷ, y)
    prob_higher = ccdf(ŷ, y)
    if Distributions.value_support(typeof(ŷ)) === Distributions.Discrete
        prob_higher += pdf(ŷ, y)
    end
    μ_lower = mean(truncated(ŷ; upper=y))
    μ_higher = mean(truncated(ŷ; lower=y))
    return (y - μ_lower) * prob_lower + (μ_higher - y) * prob_higher
end

However, this requires that mean is implemented for Truncated distributions, which currently it is only for Normal and Exponential, or that truncated returns something other than Truncated, which it does for DiscreteUniform, Uniform, and LogUniform. e.g.

julia> ŷ = Exponential(); y = 2;

julia> mean_abs_dev(ŷ, y)
1.2706705664732256

julia> mean(abs(rand(ŷ) - y) for _ in 1:100_000_000)
1.2705551367700967

julia> ŷ = DiscreteUniform(-5, 28); y = 8;

julia> mean_abs_dev(ŷ, y)
8.852941176470589

julia> mean(abs(rand(ŷ) - y) for _ in 1:100_000_000)
8.85276041

mean_abs_diff is a little trickier but can be implemented as

function mean_abs_diff(ŷ::Distributions.UnivariateDistribution, y)
    μ_lower, μ_higher = mean(JointOrderStatistics(ŷ, 2, 1:2))
    return μ_higher - μ_lower
end

I've started work on a PR that, if accepted, would implement this mean for most distributions in Distributions.jl: JuliaStats/Distributions.jl#1950

In any, case, if only Normal and DiscreteNonParametric were supported, this would be a valuable contrbution I would support. Lead time for review would be a week or so.

You can check out the BrierScore doc string to get an idea of what is required here.

Great! I'll probably take a few weeks to finish up the Distributions PR, and I might also add some more Truncated means while I'm at it so that it's easier to test the implementations. Then I'll tackle a PR here.

[ablaom]

For this specific example, yes, though when CRPS is formulated as E[|ŷ-y|] - E[|ŷ-ŷ'|] it also makes sense for discrete numeric sample spaces. e.g. if one has a posterior distribution and a Poisson observational model, CRPS could be used to evaluate the model via the posterior predictive distribution. It never really makes sense for a categorical distribution though.

Yes, of course you are right. It makes sense also for distributions with integer support.

I am happy for you to decide of the scope of this, ie how generic to make it, and also how much is practical to make a responsibility of Distributions, etc. Good luck.