Future option: Express mass/normalization of densities

[phipsgabler]

I've been previously thinking about whether and how a density interface should include the possibility to require or check for normalization, so I'm just posing this for discussion here.

Options that come to my mind:

# extra argument and dispatch
logdensityof(d, x; normalized=Val{false}())

Or maybe include a whole machinery for dealing with normalization constants?

# like eltype/length traits
hasnormalizingconstant(d) = false
normalizingconstant(d::SomeDensity) = ...
isnormalized(d) = normalizingconstant(d) == 1

# probably should be its own struct
normalize(d) = let logZ = log(normalizingconstant(d))
    logfuncdensity(x -> logdensityof(d, x) - logZ)
end

[oschulz]

Hm, I wonder if something like this could go into Distributions?

[phipsgabler]

The real use case of this IMHO is the query "is something a probability density". But that should then be available to PPLs too, since there we have all kinds of densities that we know are not normalized. So if this is deemed a good idea, I'd rather have it here.

[oschulz]

The real use case of this IMHO is the query "is something a probability density".

What do you think about my proposal at the (current) end of #5 regarding adding ismeasure and isdensity, as more specific queries than hasdensity? Would that serve this use case as well?

[phipsgabler]

I wouldn't think that makes a difference; being a measure doesn't tell you anything about normalization (or being a probability measure).

A third approach, where you only have to define normalizingconstant:

hasnormalizingconstant(d) = normalizingconstant(d) isa Nothing
normalizingconstant(d) = nothing
normalizingconstant(d::SomeDensity) = ...
isnormalized(d) = normalizingconstant(d) == 1

[oschulz]

We could call it logdensitynorm, and default it to logdensitynorm(d) == missing. Not sure we'd need a haslogdensitynorm, with a default of missing.

[oschulz]

Having a logdensitynorm could be useful in applications that need to operate in non-log space (e.g. subvolume integration) but need to renormalize (shift log-space) before doing exp(logdensityof(...)) because they would be out of the dynamic range of their float-type otherwise.

[phipsgabler]

I think I like that idea. It seemed to me that it would be nice also be able to express "infinite norms" (as there is a semantic difference between "finite but unknown constant", as in posteriors, and "known to be non-normalizable"). With your approach, it would work out, right?

logdensitynorm(d) = missing
logdensitynorm(d::Distribution) = 1
densitynorm(d) = exp(logdensitynorm(d))

will propagate the missing, as do subsequent computations. And then we can also have

densitynorm(d::Lebesgue) = Inf
# default: logdensitynorm(d::Lebesgue) = missing

with somewhat consistent behaviour.

So it seems nice on the first glance, but I'm not confident enough to judge downstream effects.

[oschulz]

Yes, I think so. But maybe the name logdensitynorm would be weird when used on measures/distributions - their norm is independent from which density is derived for them using different base measures. May we should just use LinearAlgebra.norm, and add lognorm on our own?

Then, for Bayesian inference, we'd have lognorm(∫exp(log_likelihood, prior)) == log_evidence.

[oschulz]

Hm, maybe using LinearAlgebra.norm is a bad idea, though, it's not specific enough. One might expect that norm(f_density) == sqrt(∫f_density(x)^2 dx), or something else. There are so many norms ...

[phipsgabler]

Yeah, I really wouldn't use LinearAlgebra.norm; for distributions, it's also already taken (JuliaStats/Distributions.jl#1321). And for that reason, if something like this is implemented, we should probably refrain from "norm" as the name.

[oschulz]

What could we call it so that it would also make sense for densities (as the integral of the density function over it's support)?

[oschulz]

What is the norm of MeasureTheory.Dirichlet([2,3]), and can I compute it with `MeasureTheory?

[cscherrer]

"Norm" and "normalize" are very overloaded, which is why we went with mass and logmass
cc @sethaxen

[oschulz]

I like those. I think defining those as abstract function in DensityInterface could be very handy. We often do use the term "norm of a density", we say the "PDF is normalized" and so on. So it's not just useful for MeasureBase, it's also useful for densities (density functions). We do of course mean the integral in respect to an implied base measure - but that's the same as using an implied base measure in `[log]densityof], so it would fit in quite neatly.

What do you think about this, @devmotion ?

[cscherrer]

Here's some more context:
JuliaMath/MeasureTheory.jl#128

For a trait-based approach, I can see a few situations we'd want to distinguish:

1. Measures that are normalized
2. Measures with a known mass, which may or may not be one
3. Measures known to be finite, with unknown mass
4. Measures known to have infinite mass
5. Measures where the mass is unknown (maybe finite, maybe infinite)

[oschulz]

Do we have a good package for static numbers that can represent infinity? Then logmass(measure) could deliver all of this - return either static zero (we could take that from Zeros.jl), a real number, FiniteMass() (not sure we need that one), a static Inf or missing.

All of those results could be dispatched upon, giving you that trait-based approach.

[cscherrer]

I think I'd go with Static.jl and Infinities.jl, so static(∞)

[oschulz]

Nice!

We wouldn't need these as dependencies for DensityInterface, even - we would just define

function logmass(object)
    if densitykind(object) <: IsOrHasDensity
        missing
    else
        throw(ArgumentArror("object is not and does not have a density")
    end
end

Implementations of densities and measures can then return a real number or a static something as fits the case.

[cscherrer]

When the user gives a function and says "treat this as a density", they need to be able to specify which of these cases we're in:

1. Measures that are normalized
2. Measures with a known mass, which may or may not be one
3. Measures known to be finite, with unknown mass
4. Measures known to have infinite mass
5. Measures where the mass is unknown (maybe finite, maybe infinite)

Otherwise we'll end up wasting lots of cycles on quadrature.

[oschulz]

Maybe this goes a bit for things "that are a density or a measure", though? This may be better situated in MeasureBase.

[cscherrer]

So things defined outside of MeasureBase will just be treated as completely unknown, and always use quadrature to get an approximate mass? That seems very limiting.

[oschulz]

So things defined outside of MeasureBase will just be treated as completely unknown

Well, there's isn't that much, right now: Distributions are normalized, and densities created for logfuncdensity don't have mass information in the first place, and code handling measures should use MeasureBase (at least).

Other densitiy types can have a mass, too, of course - but the best place to take care of that it is the the package that provides that density type, I would say.

Fully generic and automatic quadrature is hard to provide anyhow - there are many different algorithms and not all suit every situation. If the mass has to be estimated numerically, this won't work without involving the users in higher-dimensional cases, I would say.

[cscherrer]

Ok, so we need to be able to specify these things somewhere. It sounds like that will need to go into MeasureBase, and users of logfuncdensity will not have any way of specifying the mass. Is that right?

[oschulz]

I thought we can maybe add something like

mass(object) = missing
logmass(object) = missing

to DensityInterface (though the function name mass is probably too generic, should find a different name that can be exported), so users can specify a mass for a density if they know it (it would be the "integral in respect to an implicit base measure", of course :-) ).

The more detailed stuff I would propose to leave to MeasureBase.

[cscherrer]

How about massof and logmassof?

users can specify a mass for a density if they know it

I don't see how this would work. The type produced by logfuncdensity is owned by DensityInterface, right? So how could users add methods without type piracy?

[oschulz]

Yes, massof and logmassof could work.

users can specify a mass for a density if they know it
I don't see how this would work. The type produced by logfuncdensity

No, indeed - logfuncdensity is a quick convenient way for users to turn a log-density function into a density object, nothing more. If users want to do something more complex, they should define a type and implement (at least) isdensity and logdensityof for it - and in this case, logmassof too, for example.

[phipsgabler]

Of course @cscherrer has already put thought into this before us 😀

I like (log)massof, it's consistent with the density functions, meaningful, and distinct from norm.

Of the five possibilities Chad mentioned, all except (3) are simple to handle: values of either missing or static(something). Essentially instances of a disjoint union between reals, infinity, and Missing.

(3), however, is an outlier -- a set-like restriction of the above ("any finite real number") and harder to represent by types. I can see why it's very useful, though. I suppose an additional type like Finite, corresponding to Infinity, could solve that, but is that really a good option?

I think I'd go with Static.jl and Infinities.jl, so static(∞)

You mean StaticValues.jl?

And what is logmassof supposed to be when massof is infinite? Error-throwing? Missing?

[oschulz]

And what is logmassof supposed to be when massof is infinite?

I would say infinite as well, since log(Inf) == Inf and exp(Inf) == Inf, so that would be well-defined.

[phipsgabler]

Oh, right... for whatever reason I have assumed that log(Inf) is undefined. So that's a non-issue, sorry.

[oschulz]

No worries. :-)

[phipsgabler]

Oh, I'm using "probability measure" as "the mass is known to be 1", and would only use it as a contingent notion. Are you referring to usages such as "unnormalized probability measure"? Of course, if you understand it that way, the four properties are overdetermined.

[cscherrer]

Right, things get a lot easier if f is known to integrate to one. I thought this discussion was mostly about other cases?

[phipsgabler]

Yes, but then I don't get your point about the "overdetermined system".

We can have q(x) = exp(-x^2), "implicitely with respect to Lebesgue", and pretend not knowing the integration constant. Then the real mass Z is to be determined.

Or we have p(x) = (1/Z) * exp(-x^2), "implicitely with respect to Lebesgue", and massof(p) = 1.

Are you talking about the interpretation of "q is a probability measure over the reals and we don't know Z"? Then it might indeed look like a conflict between

• q against Lebesgue(), with mass Z or
• q against (1/Z) * Lebesgue(), with mass 1.

But I don't see such a conflict arising. If an density f is a measure, then the mass Z is explicitly defined through the base measure mu. If not, then the base measure of f is assumed to be an implicit but fixed mu. There is then also a mass Z with respect to mu. In both cases, we might construct a new density, f / Z, which with respect to mu has mass 1. We just can't talk about both Z and 1 separately.

If that's not what you're talking about, please give me a counterexample :D

[phipsgabler]

I suspect (but am not sure) that our misunderstanding arises because I for now think about the density objects first, with the base measure just being additional information telling us about "support" (really just "PMF or PDF?"), while you are always starting from measures first and most importantly, and the RD-derivative resulting as a consequent therefrom. But that is because I try to put myself into the head of the possible measure-theoretically unaware end user.

[oschulz]

Oh, but maybe that's just
https://github.com/JuliaCollections/Memoize.jl

Yes, cached results could be useful to handle "properties" of things that can be expensive to calculate while still keeping things immutable (semantically).

[oschulz]

As for examples: posteriors

Ah, but posteriors are measures, not densities. So they have a mass that doesn't depend on the choice of a base measre. So they should implement the MeasureBase interface, like Distributions hopefully will, soon. Correct, @cscherrer ?

[cscherrer]

Thanks @phipsgabler . I agree with most of this, but a couple of pieces are very confusing.

If an density f is a measure

Those are very different things

then the mass Z is explicitly defined through the base measure mu.

No, it's easy to define a measure without knowing what it integrates to.

If not, then the base measure of f is assumed to be an implicit but fixed mu. There is then also a mass Z with respect to mu. In both cases, we might construct a new density, f / Z, which with respect to mu has mass 1. We just can't talk about both Z and 1 separately.

This sounds right to me.

Let's step back for a second. If a user gives a density fit's easy to imagine creating a new density, say normalized(f). Maybe that's the idea? If this value happens to be known (user-provided or computed using quadrature) then normalized(f) can actually be computed. Otherwise it would need to be a "known up to a constant" kind of thing.

[oschulz]

If a user gives a density fit's easy to imagine creating a new density, say normalized(f)

But if I can do that, I've kind of turned my density into a measure already, right? After all, to normalize it I have to "measure" it.

[oschulz]

@phipsgabler I guess the question is, are there use cases where we would want to do something like integrate a likelihood (without the prior).

[cscherrer]

But if I can do that, I've kind of turned my density into a measure already, right?

Even before that! If you (implicitly) assume Lebesgue base measure, then the density already determines a measure.

[oschulz]

Right, that's what I meant. It's definitely a useful concept get the "default" density of a measure (or we wouldn't have Distributions.logpdf :-) ). But in the other direction - I'm not so sure. In the PPL case, you're already building a measure (the user doesn't need to know that), so defining a massof makes perfect sense. But regarding mass of densities I'm not sure.

DensityInterface provides a way to get log/non-log density functions for both density objects and measure objects (the latter using an implied base measure). One typical use case is the (not uncommon, in practice) formulation of Bayes theorem that use density functions for likelihood, prior and posterior (with an unconsciously implied common base measure like Lebesgue).

DensityInterface also provides a way to construct density object from a log/non-log density function. This allows users to define (just as an example) a likelihood, a representation of a (density) function that can be passed around without having to track if it's defined via a log- or non-log-function.

The question is, I think: Should DensityInterface provide a way to construct a measure object from a log/non-log density function (using an implied base measure). Maybe no, because taking a measure (in [log]densityof) and extracting something (a log/non-log function) from it reduces information. But making a measure from a log/non-log density function means adding information, and there's no API in DensityInterface to deal with measures beyond their default density. But if there are use cases, we may want to consider it.

[oschulz]

Hm, maybe two thought experiments:

1. DensityInterface it not meant to be limited to statistics. Let's take a use case that deals with inhomogeneous objects. The are not densities, but they have a (mass) density. Since we have a specific application (we're not interested in charge densities or so), densityof(object, x) makes sense - the base measure is implied by the type of the object. The objects also have a mass. But an object's density doesn't have a mass, I would say.

2. A probability measure has a "probability mass" - for a distribution its value is one. It also has a "default" probability density. We often do say that a PDF is normalized - so we imply that a PDF does have a mass (we really mean the probability measure, but people often to talk about "normalizing a PDF").

But the thing is - given just a plain density function, I can assume a sensible base measure and integrate it. But - have I created a measure or a physical object? Given just a function, there's no context to build it, and DensityInterface doesn't have the type(s) to represent it.

[phipsgabler]

Ah, but posteriors are measures, not densities. So they have a mass that doesn't depend on the choice of a base measre. So they should implement the MeasureBase interface, like Distributions hopefully will, soon.

Yeah, but they always have densities, right? Because that's how we use them -- evaluating their densities at points. The measure always exists in a theoretical sense, but often, we can't automatically construct it, for technical reasons (it requires static analysis of a probabilistic program etc.). But you could, e.g., have the user explicitely add the base measures they intend (which will be Lebesgue and counting for contininuous and discrete parameters, right?).

If an density f is a measure

Those are very different things

Ah, yeah, sure -- sorry for being sloppy again... what I mean is of course "in the context of MT, where f is he RD derivative of a known measure". I think then the rest should make sense, too.

Let's step back for a second. If a user gives a density fit's easy to imagine creating a new density, say normalized(f). Maybe that's the idea? If this value happens to be known (user-provided or computed using quadrature) then normalized(f) can actually be computed. Otherwise it would need to be a "known up to a constant" kind of thing.

That's exactly what I'm thinking of, yes! Except that if the massof is missing/Unknown(), I would let normalize fail immediately.

---

Now, for the following I fear we stray into fundamental discussions again, but let me answer one point...

But the thing is - given just a plain density function, I can assume a sensible base measure and integrate it. But - have I created a measure or a physical object?

The question is, I think: Should DensityInterface provide a way to construct a measure object from a log/non-log density function (using an implied base measure).

I see what you mean there, but I think that's the line we should not cross in DensityInterface. IMHO it is good to set up things so that the "implicitly existing base measure" view is always sane and can be justified, but not to actually do (and thereby "force" unto users) the measure theoretic notions. When you want to get the "real measure" from the density object, do it using MeasureTheory.

Also, by just saying that we talk about masses doesn't force us to commit to a full interface for integration. Also that should, IMHO, be left to either MeasureTheory proper, or the packages doing the actual integration (which may then depend only on the density object and their own extra information, like base measure or support.)

---

That being said -- my naive idea was to just define the massof interface with a default of Unknown/missing. Interface implementors then may opt in to provide more information, if known -- with respect to the "implicit base measures" view (1 for Distributions, other finite numbers for measures, Finite or so for probabilistic programs). If the client happens to do proper measure theory, this should of course be consistent wrt the explicit base measure.

normalize is then only a helper function that can scale a given "density thing" (either a density object, or a measure-like object), when a known mass exists. The result should then be consistent again: a density object with the same implicit base measure, or a new scaled measure, etc. (Distributions can just return themselves, for example). Really just "OK, you told me my mass is Z, so obviously my normalized density is my own devided by Z".

If that is not a reasonable idea -- what, for example, is the current design of MeasureTheory's mass functions?

[cscherrer]

Thanks @phipsgabler . Normalization will be really important for MeasureTheory, so I want to be sure the interface is MT-friendly. We haven't implemented mass yet, but @sethaxen and I have talked about it. I'll give my view of the design ideas.

When a user builds a new measure, they can either specify the mass or leave it to be determined. In some cases the code should be able to prove some things about this. For example, if you start with a finite measure, restricting it will result in another finite measure. That kind of thing. This is the idea behind the different "classes" I proposed, and why I think a trait representation could make sense.

If the user doesn't specify a mass, but asks one, we need to be able to represent that. But we don't want to implement quadrature, or impose a specific method, or depend on every package we know about. So the result here would probably be a new Integral type (or Integration, or DefiniteIntegral, something like that). What we really need is a rigorous and computationally efficient + C from intro calculus.

The user should also be able to ask for a normalization of a given finite measure. Here again, we'd probably have a first-class representation like a Normalized type. Evaluating the log-density of this measure must either include the + C indeterminacy, or be over (more likely) be over a base measure that includes this indeterminacy.

An important aspect of this is that individual evaluations of log-densities include the indeterminacy, but taking differences of these must allow it to be canceled out. In this way it feels a little bit like the provenance tracking done in Measurements, that allows making sure that x - x == 0 for uncertain x, rather than incorrectly treating them as independent.

There's a lot to this, which is part of why we haven't gotten to it yet.

[oschulz]

@cscherrer , so you have something in mind like measure = ∫(density, μ, mass = static(1)) ?

That may be very user friendly. And if it that would accept a DensityInterface-compatible object as density, there would be no need to state log/non-log explicitly.

But: If there can be mass(density) that would make things awkward, because that mass might not be in relation to μ. And then we'd have mass(measure) != mass(density). We could do something about that if we knew the base measure of the density, as well as it's mass. But a density together with a known base measure is, in effect, already a measure.

@phipsgabler , maybe we should explore the concept of [log]massof function in MeasureBase first, and then - maybe - "promote it upwards" to DensityInterface if we see a bit clearer, and we need densities with a mass that should not be measures?

[cscherrer]

so you have something in mind like measure = ∫(density, μ, mass = static(1)) ?

Interesting, I hadn't thought of that. So here mass would be another field of DensityMeasure? I think a separate function is still better, because then it could be added after the fact.

But: If there can be mass(density) that would make things awkward, because that mass might not be in relation to μ.

Mass is a characteristic of a measure. Once you have a measure, the (base measure, density) pairs are just different ways of expressing it.

Or to put it another way... Given a density, the base measure is only relevant for determining the measure. So "the mass of a function" doesn't really make sense, you need a base measure. But the reason you need it is so the measure itself is determined, which in turn makes the mass well-defined.

[phipsgabler]

maybe we should explore the concept of [log]massof function in MeasureBase first, and then - maybe - "promote it upwards" to DensityInterface if we see a bit clearer, and we need densities with a mass that should not be measures?

I think that's a reasonable approach, if the ultimate goal of generalizability to DensityInterface and "implicit base measure densities" is always kept in mind (to avoid a situation like the retrofittting of logdensityof).

[phipsgabler]

@cscherrer about the possible return values -- would that then be

- Value known
  > x or static(x) for finite x
  > static(∞)
- Value unknown
  > DefiniteIntegral(μ, basemeasure(μ)) (unknown finite)
  > Unknown/missing (unknown)

?

[cscherrer]

This sounds at least pretty close. The completely unknown case could still be represented as an integral to be passed to outside quadrature routines. Then we still need to distinguish the "I only know this is finite" case.

There are lots of different types here. Should everything be under a hierarchy? There are some tradeoffs, but it could be useful for dispatch. Maybe something like

abstract type AbstractMass end

struct UnknownMassOf{T}
    measure :: T
end

struct UnknownFiniteMassOf{T}
    measure :: T
end

struct KnownFiniteMass{T}
    mass :: T
end

This is a rough sketch, but in the unknown cases we would need to refer to a fully-determined measure (Maybe density function over Lebesgue, but "density function over unknown base measure" isn't enough).

If the mass is known, I think there's no need to keep the reference to its origin.

[oschulz]

so you have something in mind like measure = ∫(density, μ, mass = static(1)) ?
Interesting, I hadn't thought of that. So here mass would be another field of DensityMeasure? I think a separate function is still better, because then it could be added after the fact.

Yes, I maybe a "decorator" approach is more composable/flexible here. So measure = withmass(∫(density, μ), static(1))), resulting in something like a MeasureWithMass{RadonNikodymIntegral{...},...}.

The completely unknown case could still be represented as an integral to be passed to outside quadrature routines.

I guess a lazy value could be returned in that case (so that integration is only run if needed)? I think we have packages for lazy values, but I haven't looked into that area for a while.

There's still the question of algorithm choice. Maybe it would be better to return a trait that says "it has a norm, but the value is unknown, you need to integrate". Then the user can integrate with a suitable algorithm, and create an new measure object that includes the mass information, e.g. via withmass(measure, mass_value) or so (see above). An advantage would be that we wouldn't need to depend on integration packages. In actual non-trivial use cases, default integrator choices would likely fail anyway.

[cscherrer]

Any updates here? I think we should commit to this being in MeasureBase or MeasureInterface. We've already complicated things there to conform with DensityInterface, and it will quickly get out of hand and be hard to reason about if we push that farther.

[oschulz]

Yes, having it in the measure theory ecosystem will be easiest for now. If we see a need and a path to move it to DensityInterface later on (one the dust around Distributions and MeasureBase has settled) we can still reopen this. Ok with you, @phipsgabler ?

[phipsgabler]

I agree for now. But in the long term, I'd absolutely want to get this working for measure-oblivious applications.

[oschulz]

Let's keep this open as a reminder for the future then.