Interoperability with MeasureTheory

[cscherrer]

For use with MeasureTheory, it would be helpful if we could write logfuncdensity(log_f, base), with expected behavior

logdensityof(logfuncdensity(log_f, base), x) == log_f(x) - logdensityof(base, x)

This would allow us to write

DensityInterface.logfuncdensity(log_f, base::Measure) = MeasureTheory.∫exp(log_f, base)

I know this package doesn't want base measures to be a requirement, so maybe we could find a way to make base optional.

∫exp is also defined for Distributions; it's not clear to me yet how to extend the interface to cover that case.

[oschulz]

Hm, do we want logfuncdensity to act an "alias" (for two arguments) if MeasureTheory (resp. MeasureBase) already provides ∫exp?

If the second argument will be a Measure, then MeasureBase seems to be the natural place to define the function (after all, it already does). Why not have different names for the two functions? The the case with the measure is more specific (which we want, if course, if a Measure is available) so maybe it's good if it has a less generic name.

It would be great if the return value of MeasureTheory.∫exp would support DensityInterface.hasdensity and DensityInterface.logdensityof, though. I think that would already give us interoperability.

[cscherrer]

The issue is that logfuncdensity(log_f) doesn't make sense at all for MeasureTheory, because a function alone isn't enough to specify a measure. I guess we could just not use that function, but I think a desirable quality of an interface is for it to make sense across packages.

[oschulz]

Yes, but shouldn't that function be called logfuncmeasure(log_f, measure) instead of logfuncdensity, since the return value is a measure, not a density (whereas logfuncdensity(log_f) will really return "just" a density)?

[oschulz]

logdensityof is flexible on purpose, regarding it's input (can be a density or something that has a base measure), to help with generic programming. But for construction of densities and measure, shouldn't we be more specific?

[cscherrer]

I thought it went

• logdensityof(d) returns a (a density) function
• logfuncdensity(f) takes a function, treats is as a density, and returns a measure (even if you don't call it that)

If not, then what does logfuncdensity do?

[oschulz]

In general logfuncdensity just returns a density that is defined by a log-function. Like a likelihood density, defined by a log-likelihood function.

Lets say we have a distribution d. logdensityof(d) will return a function that computes the log of the probability density. But given an arbitrary "black-box" function log_f that computes a the log of a probability density at given points, logdensityof can only create a density object (the probability density), not a distribution - not enough information.

Except if log_f was the result of calling logdensityof(d), which returns Base.Fix1(logdensityof, d), so the information is available.

I guess we should say logfuncdensity can return a measure if there's an implied base measure for log_f - but in most cases, it will just return a density.

[cscherrer]

I thought these were inverses, how can that be if they both return densities?

[oschulz]

How about this - after we add this to Distributions (JuliaStats/Distributions.jl#1416):

@inline DensityInterface.hasdensity(::Distribution) = true

DensityInterface.logdensityof(d::Distribution, x) = logpdf(d, x)

DensityInterface.densityof(d::Distribution, x) = pdf(d, x)

we add something like this to MeasureBase (caveat - I didn't test this, and some of it may already be there):

const MeasureLike = Union{MeasureBase.Measure, Distributions.Distribution}

# ParameterFreeMeasure means a Measure subtype that's singleton type:
abstract type ParameterFreeMeasure <: MeasureBase.Measure end

struct LebesgueMeasure <: MeasureBase.ParameterFreeMeasure end
struct CountingMeasure <: MeasureBase.ParameterFreeMeasure end


MeasureBase.logdensity(
    μ::Distributions.ContinuousDistribution, ::LebesgueMeasure, x
) = logpdf(μ, x)

MeasureBase.logdensity(
    μ::Distributions.DiscreteDistribution, ::CountingMeasure, x
) = logpdf(μ, x)


struct LogFuncRDIntegral{LogF,M<:MeasureLike} <: MeasureBase.Measure
    log_f::LogF
    μ::M
end

# Do whatever MeasureBase/MeasureTheory currently does here:
MeasureBase.logdensity(μ::LogFuncRDIntegral, ν::Measure, x) =
    MeasureBase.logdensity(MeasureBase.∫exp(μ.log_f, μ.μ), ν, x)


function logfuncmeasure end

MeasureBase.logfuncmeasure(log_f, μ::MeasureLike) =
    MeasureBase.LogFuncRDIntegral(log_f, μ)

MeasureBase.logfuncmeasure(
    log_f::Base.Fix1{
        Union{typeof(logdensityof),typeof(logpdf)},
        <:Distributions.ContinuousDistribution
    },
    ::LebesgueMeasure
) = log_f.x

MeasureBase.logfuncmeasure(
    log_f::Base.Fix1{
        Union{typeof(logdensityof),typeof(logpdf)},
        <:Distributions.DiscreteDistribution
    },
    ::CountingMeasure
) = log_f.x


function basemeasure end

MeasureBase.basemeasure(m::MeasureBase.LogFuncRDIntegral(m)) = m.μ

MeasureBase.basemeasure(m::Distributions.ContinuousDistribution) =
    MeasureBase.LebesgueMeasure()

MeasureBase.basemeasure(m::Distributions.DiscreteDistribution) =
    MeasureBase.CountingMeasure()


struct RadonNikodymDerivative{Mu,Nu}
    μ::Mu
    ν::Nu
end

DensityInterface.hasdensity(::RadonNikodymDerivative) = true

# If the μ is a singleton, we can do logfuncmeasure based on type:
MeasureBase.logfuncmeasure(
    log_f::Base.Fix1{typeof(logdensity),RadonNikodymDerivative{Mu,Nu}},
    μ::Nu
) where {Mu,Nu<:ParameterFreeMeasure} = log_f.x


function rd_derivative end

MeasureBase.rd_derivative(ν::MeasureLike, μ::MeasureLike) =
    MeasureBase.RadonNikodymDerivative(ν, μ)

# If the ν is a singleton, we can do rd_derivative based on type:
MeasureBase.rd_derivative(
    μ::LogFuncRDIntegral{LogF,M}, ν::M
) where {LogF,M<:ParameterFreeMeasure} = DensityInterface.logfuncdensity(ν.log_f)



DensityInterface.logdensityof(m::LogFuncRDIntegral, x) =
    m.log_f(x) + DensityInterface.logdensityof(m.μ, x)

DensityInterface.logdensityof(m::MeasureBase.RadonNikodymDerivative, x) =
    MeasureBase.logdensity(m.μ, m.ν, x) 
    # Maybe rename `MeasureBase.logdensity` to `MeasureBase.eval_log_rd_derivative` or so
    # to make the function names less similar?


DensityInterface.logdensityof(m::MeasureBase.Measure) =
    DensityInterface.logdensityof(
        MeasureBase.rd_derivative(m, MeasureBase.basemeasure(m))
    )

DensityInterface.logdensityof(m::MeasureBase.Measure, x) =
    DensityInterface.logdensityof(m)(x)


DensityInterface.logdensityof(
    m::RadonNikodymDerivative{<:Distributions.ContinuousDistribution,LebesgueMeasure},
    x...
) = logdensityof(m.μ, x...)

DensityInterface.logdensityof(
    m::RadonNikodymDerivative{<:Distributions.DiscreteDistribution,CountingMeasure},
    x...
) = logdensityof(m.μ, x...)

(I'm prefixing functions and types with package name a lot here for clarity).

Then this should work:

using DensityInterface, MeasureBase

μ = LebesgueMeasure() # Or another ParameterFreeMeasure compatible with log_f
m = logfuncmeasure(log_f, μ)

logdensityof(rd_derivative(m, μ), x) == log_f(x)
logdensityof(rd_derivative(m, μ)) === log_f

logdensityof(m, x) == log_f(x)
logdensityof(m) === log_f

It would the measure-theoretically explicit analog of (this would also work)

d = logfuncdensity(log_f)

logdensityof(d, x) == log_f(x)
logdensityof(d) === log_f

And this should also work:

d = SomeContinousDistribution()

log_f = logdensityof(rd_derivative(d, LebesgueMeasure())
logfuncmeasure(log_f, LebesgueMeasure()) === d

log_f = Base.Fix1(logpdf, d)
logfuncmeasure(log_f, LebesgueMeasure()) === d

It would be the measure-theoretically explicit analog of (this would also work)

log_f = logdensityof(d)
logfuncdensity(log_f) === d

log_f = logpdf(d)
logfuncdensity(log_f) === d

One could also think about making the measure theory API at least partially trait-based and replace MeasureLike by Any. Then I could define something like this in BAT:

MeasureBase.basemeasure(m::BAT.NonNormalizedPosterior) = m.prior

function MeasureBase.rd_derivative(μ::BAT.NonNormalizedPosterior{Li,Pr}, ν::Pr) where {Li,Pr}
    ν == ν.prior || throw(ArgumentError("ν must equal prior in rd_derivative of posterior"))
    ν.likelihood
end

(Note: The type is actually called BAT.PosteriorDensity, I use BAT.NonNormalizedPosterior here for clarity - though maybe I should rename it to that at some point).

Maybe similar things would be possible in Turing/PPL?

Sorry, I know this is a lengthy post. On the plus side, I think the code above may be quite complete already.

@cscherrer, what do you think of this?

@devmotion and @phipsgabler, how do you see this from a Distributions and Turing/PPL point of view?

[cscherrer]

Thanks for working through this @oschulz . I'll have a look and discuss with @mschauer, and get back to you soon.

[oschulz]

Thanks for working through this @oschulz . I'll have a look and discuss with @mschauer, and get back to you soon.

Thanks!

Regarding

[logdensityof/logfuncdensity] I thought these were inverses, how can that be if they both return densities?

logdensityof returns a function that computes the log-value of something that is or has a density. So it does not return a "density object", it returns a log-function (and logdensityof(logdensityof(d)) would not make sense).

logfuncdensity, on the other hand, takes a log-function (assumed to be the log-function of a density) and returns a corresponding "density object" (which is usually not a function, and even if it is, that function should certainly not compute the log-density value).

Actually, to make my proposal above work nicely, we should say that logfuncdensity returns a density or something that has a density (e.g. a measure). The latter will only happen in special(ized) cases, though (e.g. when it's input log-density function originated from a distribution or measure and still carries that information in a dispatchable way).

A major aim of DensityInterface is to let people work with "density objects" that do not have to be named log_... in the code to make clear if it's (already) a log or not:

It's so awkward to write log_likelihood = some_func; prior = SomeDistribution(). I had that in examples during the early days of BAT.jl and several users asked "why one with log_ and not the other?". On the other hand, likelihood = some_func; prior = SomeDistribution() is misleading, dangerous and wrong, if some_func returns the log of the likelihood.

But likelihood = logfuncdensity(some_func); prior = SomeDistribution() is clear to the reader. And it's safe, the code using likelihood and prior (and the posterior, and so on) can use logdensityof (resp. densityof) and be sure to get the log (resp. not the log).

[oschulz]

Maybe we actually should add something like

radon_nikodym_integral(d, μ) = logfuncmeasure(logdensityof(d), μ)

to my proposal above, to take advantage of the fact that things that implement the DensityInterface API (d in this case) don't need to have log in the name (on the LHS here).

To make things even nicer, we could replace the log_f-based struct LogFuncRDIntegral in my proposal with

struct RadonNikodymIntegral{F,M<:MeasureLike} <: MeasureBase.Measure
    d::D
    μ::M
end

requiring d to support logdensityof (can e.g. be the result of logfuncdensity(log_f)).

RadonNikodymIntegral would form an elegant symmetry with RadonNikodymDerivative, name-wise.

[oschulz]

I updated the proposal above a bit, renaming LogFuncMeasure to LogFuncRDIntegral and adding

DensityInterface.logdensityof(m::LogFuncRDIntegral, x) =
    m.log_f(x) + DensityInterface.logdensityof(m.μ, x)

Would indeed by nicer with a struct RadonNikodymIntegral, I think, then it could be

DensityInterface.logdensityof(m::RadonNikodymIntegral, x) =
    DensityInterface.logdensityof(m.d, x) + DensityInterface.logdensityof(m.μ, x)

[cscherrer]

Hi @oschulz , just going through this...

In a few places, you define a logdensity in terms of the three-argument method. This was how we started out, but it doesn't work well because

1. It forces the base measure to always be explicit. Especially for new users, it's much simpler if the base measure can be a separate function associated with a measure
2. This requires three-argument methods for all pairs of measures, which is intractable.

Instead we define

logdensity(::SomeMeasure, x)
basemeasure(::SomeMeasure)

The three-argument methods are then derived from this. If you take every possible measure and repeatedly apply basemeasure, you get a tree. The three-arugment logdensity is in terms of walking this tree. This approach allows us to limit ourselves to methods between "primitive" measures for definition of three-argument logdensity methods.

---

DensityInterface.logdensityof(
    m::RadonNikodymDerivative{<:Distributions.ContinuousDistribution,LebesgueMeasure},
    x...
) = logdensityof(m.μ, x...)

We can't just ignore the base measure here, that will lead to incorrect results. Maybe we can talk about the x... in a call, I don't understand what's intended here.

A ContinuousDistribution doesn't necessarily have Lebesgue as its base measure. For example, a Dirichlet in three dimensions is a two-dimensional manifold (with boundary), and is singular with respect to Lebesgue measure.

---

logdensityof returns a function that computes the log-value of something that is or has a density. So it does not return a "density object", it returns a log-function (and logdensityof(logdensityof(d)) would not make sense).

logfuncdensity, on the other hand, takes a log-function (assumed to be the log-function of a density) and returns a corresponding "density object" (which is usually not a function, and even if it is, that function should certainly not compute the log-density value).

Actually, to make my proposal above work nicely, we should say that logfuncdensity returns a density or something that has a density (e.g. a measure). The latter will only happen in special(ized) cases, though (e.g. when it's input log-density function originated from a distribution or measure and still carries that information in a dispatchable way).

We can't define a logfuncdensity like this without knowing a base measure to use for it. That's not just a MeasureTheory thing. In general, you either specify a base measure, or you sacrifice expressiveness (say, assuming everything is with respect to Lebesgue measure).

The function-like object passed to logfuncdensity might "know" the appropriate base measure. But it seems dangerous to not have that available as part of the interface.

---

A general observation... We've talked about the goal of expressing "function-like" objects that can be evaluated in either log space or linear space. But without a base measure, it's not really a density.

In MeasureTheory we have a Likelihood type that's a lot like this. It's not a measure, but it can "act on" one measure to yield a new one using pointwise multiplication, ⊙. So for a measure μ and Likelihood ℓ,

μ ⊙ ℓ

gives a new measure. If μ is considered as a prior, then μ ⊙ ℓ is a posterior. But everything still works if μ is Lebesgue or CountingMeasure. The effect is to take enrich ℓ, so you end up with a new measure.

[oschulz]

Thanks @cscherrer , that cleared things up a bit more for me. I think we're on to something that could fit together very naturally here.

Instead we define
logdensity(::SomeMeasure, x)
basemeasure(::SomeMeasure)

Oh, I thought the three argument methods were the primary implementation for new densities. Nice, I think that will simplify integrating MeasureTheory with DensityInterface.

We do need a good name for the three-argument function though, IMHO. If something like rd_derivative and the like seems to technical, how about reldensity (for "relative density" resp. "related density") or so? Would that be a correct term?

Question: Does MeasureTheory allow for measures that don't define basemeasure? Or is that not allowed, do they all have to? The latter would make things easier here, I think.

The three-argument methods are then derived from this. [...]

So the three-argument logdensity function does not get specialized by new measure types, right?

A ContinuousDistribution doesn't necessarily have Lebesgue as its base measure. For example, a Dirichlet

we should say that logfuncdensity returns a density or something that has a density (e.g. a measure)
We can't define a logfuncdensity like this without knowing a base measure to use for it. That's not just a MeasureTheory thing.

You mean for the case where logfuncdensity returns a density, or for the case where it returns something that has a density (like for logfuncdensity(logdensityof(d::Distribution)))?

I assume it's just a problem in the latter case (returning something that has a density)? Since logfuncdensity returning something that is a density just means turning an existing log-density function into a log/non-log-independent "density object" without adding information or changing semantics?

The function-like object passed to logfuncdensity might "know" the appropriate base measure. But it seems dangerous to not have that available as part of the interface.

Yes, I was thinking about the same thing. For the case logfuncdensity(logdensityof(d::Distribution)) == d we could say that

• logdensity_of_d = logdensityof(d::Distribution) "turns" d into a density via it's base measure (RD-derivative).

• logfuncdensity(logdensity_of_d) "turns" logdensity_of_d into a measure via an inherent base measure, it remembers where it "came from" (RD-integral).

Would it be Ok to allow basemeasure to operate on both measures and (log-)density-functions? If so, I think that would solve this in a clean fashion. In this case, we'd have basemeasure(d) == basemeasure(logdensity_of_d).

Naturally, plain black-box log-density functions that don't support basemeasure would not get turned into measures by logfuncdensity, just into "density objects" (no change in what it is, just not fixed in term of log/non-log anymore).

This approach would also mean that for any object (e.g. a measure) m that supports logdensityof and basemeasure we could easily provide:

logfuncdensity(logdensityof(m)) == m
basemeasure(logdensityof(m)) == basemeasure(m)

The first of these DensityInterface already provides, by default, via

logfuncdensity(log_f::Base.Fix1{typeof(logdensityof)}) = log_f.x

The second one we can provide just as easily via

basemeasure(log_f::Base.Fix1{typeof(logdensityof)}) = basemeasure(log_f.x)

If it's Ok to say that density functions (can) have a base measure, that would consitute a strong argument for adding basemeasure to DensityInterface instead of a secondary package. It we did, would we also want/need a function hasbasemeasure, to complement hasdensity?

A general observation... We've talked about the goal of expressing "function-like" objects that can be evaluated in either log space or linear space. But without a base measure, it's not really a density.

You mean because it's not clear what space it lives on?

Does that mean I can't do μ ⊙ ℓ if ℓ has no base measure? Or would the base measure of ℓ just be inferred from the base measure of μ?

The latter is actually what BAT.jl does right now: The user can pass the likelihood as a plain log-function (wrapped in logfuncdensity), and the space it lives on (e.g. degrees of freedom) are inferred from the prior, if not specified for the likelihood.

In MeasureTheory we have a Likelihood type that's a lot like this. It's not a measure, but it can "act on" one measure to yield a new one using pointwise multiplication, ⊙. So for a measure μ and Likelihood ℓ,

Is Likelihood required to provide base measure information?

If we go this road: I think the name of logdensityof fit it's semantics. But I do realize that logfuncdensity is a bit of an unfortunate name for something that may return something that is a measure. But maybe we can come up with a better name for "something-that-has-a-density-derived-from-this-log-density-function". It's not widely used in the ecosystem yet (just in BAT.jl so far, I think) and I'd be open to renaming it.

Just to make sure I'm understanding things correctly - which of these statements could you subscribe to, from a MeasureTheory point of view (without cringing too much)?

1. Every measure ν has exactly one (primary) base measure μ.
2. Every measure ν has exactly one "natural" density (function) f = dν/dμ, via ν's (primary) base measure μ.
3. Every (useful) density (function) f stems from a Radon–Nikodym derivative f = dν/dμ of two measures.
4. Every (useful) density (function) f can be said to have exactly one (primary) base measure. The base measure of f = dν/dμ is μ.
5. A density (function) can only be evaluated in the "context" of a measure. For a density function f = dν/dμ this can be μ or any of μ's recursive base measures.

Oh, a question in that context: What is the base measure of the Lebesgue measure? Or are there "root" measures that don't have a base measure?