Future option: Express mass/normalization of densities

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

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

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.

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?

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

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

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.

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.

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.

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 ...

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.

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)?

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

@sethaxen

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

@devmotion

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 ?

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:

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

@phipsgabler

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.

@cscherrer

@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?

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.

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).

@cscherrer

@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)

?

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.

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.

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.

@phipsgabler

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 ?

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

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

Navigation Menu

Repository navigation

Future option: Express mass/normalization of densities #8

39 remaining items

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

Future option: Express mass/normalization of densities #8

Description

Activity

oschulz commented on Nov 10, 2021

phipsgabler commented on Nov 10, 2021

oschulz commented on Nov 10, 2021

phipsgabler commented on Nov 10, 2021

oschulz commented on Nov 10, 2021

oschulz commented on Nov 10, 2021

phipsgabler commented on Nov 13, 2021

oschulz commented on Nov 13, 2021

oschulz commented on Nov 13, 2021

phipsgabler commented on Nov 13, 2021

oschulz commented on Nov 13, 2021

oschulz commented on Nov 13, 2021

cscherrer commented on Nov 13, 2021

oschulz commented on Nov 13, 2021

cscherrer commented on Nov 13, 2021

39 remaining items

cscherrer commented on Nov 16, 2021

oschulz commented on Nov 16, 2021

cscherrer commented on Nov 16, 2021

phipsgabler commented on Nov 16, 2021

phipsgabler commented on Nov 16, 2021

cscherrer commented on Nov 16, 2021

oschulz commented on Nov 16, 2021

cscherrer commented on Nov 23, 2021

oschulz commented on Nov 24, 2021

phipsgabler commented on Nov 24, 2021

oschulz commented on Nov 24, 2021

Metadata

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

Issue actions