Normalizing a measure

[sethaxen]

In JuliaManifolds/ManifoldMeasures.jl#7 we were discussing that we need a way to map an (unnormalized) Hausdorff measure to the corresponding probability measure by normalizing. This seems general enough that it could be defined here.

What about something like this?

LinearAlgebra.normalize(μ::AbstractMeasure) = inv(total_mass(μ)) * μ
unit() = true  # just a function to use for dispatch.
total_mass(μ::AbstractMeasure) = ∫(unit, μ)
total_mass(μ::AbstractWeightedMeasure) = exp(μ.logweight) * total_mass(μ.base)
is_probability_measure(μ::AbstractMeasure) = isone(total_mass(μ))

Perhaps there is a more measure-theoretic term than "normalization", and we should not then overload normalize. In the case of a Hausdorff measure, total_mass would compute the volume/area of the manifold in some embedded metric space.

It might be preferable to define log_total_mass instead.

[cscherrer]

Great idea! A few comments:

First, ∫ has changed around a bit. It had originally computed the total measure (answer to your other question BTW). But I think it's more convenient as a way to define a new measure. So e.g. you can write ∫(x -> -x^2/2, Lebesgue(ℝ)) (or with do notation) to define an unnormalized Gaussian. This is integration in the sense of being dual to the Radon-Nikodym derivative 𝒹.

It might be preferable to define log_total_mass instead.

I think we can make it flexible, so you can define either. For density, we have

density(μ, ν::AbstractMeasure, x) = exp(logdensity(μ, ν, x))
density(μ, x) = exp(logdensity(μ, x))

Usually logdensity is what you define. But sometimes it can be more natural to define the density, in which case you can just overload logdensity to log the result of that. Or even have specialized methods for each. Julia is pretty great in this way :)

[sethaxen]

In that case, I think it makes sense to have an un-implemented logtotalmass, and a default totalmass(x) = Exp(logtotalmass(x)). Then, in order to normalize a measure, one would need to implement one of those two functions directly. Shall I open a PR?

[cscherrer]

Or to be more precise,

I think formally it's total_measure, but I also like the sound of total_mass... or even just mass and logmass. If this is the only way we use it, the total gives no information.

Yes, a PR would be great. Thanks for starting getting this going!

[cscherrer]

@sethaxen what do you think of having mass or logmass work like Base.sum, with an optional first argument for weighting? This was the idea with Integral:

struct Integral{F,M}
    f::F
    μ::M
end

By default it just stores the request in a struct. In general there's no clear efficient way to do integration, so the user needs to be able to select between lots of different algorithms. And maybe we have an extension package to conenct to QuadGK, etc?

[gdalle]

Following up on this, I like the idea of an interface with numerical integration. Could it be a small side package (to outsource the Quadrature dependencies) that mostly defines, say, a function evaluate(::Integral, ::AbstractDomain) ?

[cscherrer]

Maybe it could be even simpler? A measure "knows" its domain, so the package could just depend on MeasureTheory and Quadrature and only have one function Integral -> QuadratureProblem.

The only reason not to have this in MeasureTheory itself is that Quadrature.jl has lots of dependencies, and we need to keep the main package a little more focused.

[gdalle]

Yes but in that case it means that you cannot apply a measure to arbitrary sets. For instance, if you have a Poisson process on R^n but you want to simulate it on a box, you need to define a local measure restricted to that box instead of evaluating the global measure on this box subset.

[cscherrer]

Oh, I guess we're starting to split the conversation in different threads. I just made a related comment here: #130 (comment)