`rand(::SuperpositionMeasure)`

[cscherrer]

@cscherrer it seems SuperpositionMeasure does not have a rand method:

julia> μ = Dirichlet(10, 1.0)
Dirichlet(α = Fill(1.0, 10),)

julia> result = multipathfinder(μ, 1_000; nruns=20)
┌ Warning: Pareto shape k = 0.8 > 0.7. Resulting importance sampling estimates are likely to be unstable.
└ @ PSIS ~/.julia/packages/PSIS/Rf0S3/src/core.jl:268
Multi-path Pathfinder result
  runs: 20
  draws: 9
  Pareto shape diagnostic: 0.8 (bad)

julia> result.fit_distribution_transformed |> typeof
Pushforward{TransformVariables.CallableTransform{TransformVariables.UnitSimplex}, SuperpositionMeasure{Vector{MvNormal{(:μ, :σ), Tuple{Vector{Float64}, Pathfinder.WoodburyPDMat{Float64, LinearAlgebra.Diagonal{Float64, Vector{Float64}}, Matrix{Float64}, Matrix{Float64}, LinearAlgebra.Diagonal{Float64, Vector{Float64}}, LinearAlgebra.QRCompactWYQ{Float64, Matrix{Float64}}, LinearAlgebra.UpperTriangular{Float64, Matrix{Float64}}}}}}}, Static.True}

julia> result.fit_distribution_transformed |> rand
ERROR: MethodError: no method matching rand(::Random._GLOBAL_RNG, ::Type{Float64}, ::SuperpositionMeasure{Vector{MvNormal{(:μ, :σ), Tuple{Vector{Float64}, Pathfinder.WoodburyPDMat{Float64, LinearAlgebra.Diagonal{Float64, Vector{Float64}}, Matrix{Float64}, Matrix{Float64}, LinearAlgebra.Diagonal{Float64, Vector{Float64}}, LinearAlgebra.QRCompactWYQ{Float64, Matrix{Float64}}, LinearAlgebra.UpperTriangular{Float64, Matrix{Float64}}}}}}})
...

julia> result.draws_transformed
1000-element Vector{Vector{Float64}}:
 [0.12646688194244224, 0.025977504820487007, 0.012354054767482447, 0.1334736170622469, 0.03246415592570656, 0.11141269659298245, 0.024034133186899103, 0.24475990778479786, 0.1342558624335393, 0.15480118548341593]
 [0.0169086648766402, 0.09911231516043292, 0.34416889221132596, 0.10180645019569193, 0.14439126397502686, 0.11713472197731346, 0.08444424836455564, 0.07675136145856229, 0.004312850832038056, 0.010969230948412829]
 [0.204562486395576, 0.10418880831568512, 0.14665258035543918, 0.03270907277182203, 0.025198320818369802, 0.08760170721317947, 0.22661983988004777, 0.0454094719713609, 0.08236705364451188, 0.04469065863400787]
 [0.01816521614098672, 0.08744121201467416, 0.08194354696600549, 0.21065542194335074, 0.05446274719421834, 0.2462246106035293, 0.08504411121447199, 0.09401263968332611, 0.05085735912825384, 0.07119313511118336]
 ⋮
 [0.05276387864922493, 0.06859177668105768, 0.12000714757911174, 0.013786343525336437, 0.0936220599997711, 0.014944597056716297, 0.15220938553755983, 0.32035991301423006, 0.07951434052610347, 0.08420055743088842]
 [0.037353097348125405, 0.1003300343889233, 0.014071365428468915, 0.19619179988878918, 0.002071999284771993, 0.15604889644757544, 0.394789935358081, 0.017367043688102463, 0.02483761912275614, 0.05693820904440616]
 [0.092177514123393, 0.03885924399407327, 0.11357006993749118, 0.040281419867287165, 0.0017297788961319048, 0.30301174382535745, 0.08705530131888517, 0.10436044162668617, 0.12720672435099528, 0.09174776205969937]
 [0.1036299013468673, 0.0237354567117644, 0.00808528404000608, 0.08094192248237413, 0.012319402714347172, 0.17985596931554573, 0.044975357793939554, 0.09246673694172126, 0.33590939793731694, 0.11808057071611726]

Originally posted by @sethaxen in mlcolab/Pathfinder.jl#46 (comment)

[cscherrer]

Thanks @sethaxen , we really need to nail down the design for this. I think the main issue comes down to normalization. Some examples;

• Sampling from a Normal() is no problem.
• Should we be allowed to sample from 3 * Normal()?
• What about 2.0 * Normal() + 3.0 * StudentT(2)?

That might start to seem more like a clear "no", but then what about 0.4 * Normal() + 0.6 * StudentT(2)? The types can't help us distinguish these, but we want to avoid dynamic dispatch.

@theogf and I were talking recently about how a lot of what we do with smart constructors is like the usual data structure business of worrying about invariants. So one possibility would be to enforce that the weights in a superposition must add to one. Then 2.0 * Normal() + 3.0 * StudentT(2) would become 5.0 * (0.4 * Normal() + 0.6 * StudentT(2)).

My main concerns with that are (1) the overhead of maintaining the invariant, and (2) that is still leaves out that instead of a distribution we might have something like Lebesgue() as a leaf.

I should maybe mention another thing I was considering for that. rebase is defined as

rebase(μ, ν) = ∫(𝒹(μ,ν), ν)

I think for sampling this, it's reasonable to use importance sampling. Or at least, that feels "primitive" in some sense. So we would sample ν and give 𝒹(μ,ν) as the log-weight - the result could represent a sample from μ.

This leads into some other things about representing samples as measures, etc, but I'll hold off on that for now.

[sethaxen]

TBH I haven't thought about this at all from a design perspective. For the moment I just need a way to represent a uniformly weighted mixture of multivariate normals. I certainly could roll my own type for that (I'm going to keep the integration code in Pathfinder as not part of the API for now, so I'm not concerned about type changes in the future) as well instead of using SuperpositionMeasure.

But this also isn't urgent because I don't think the Pathfinder/MeasureTheory integration will be used much until Soss/Tilde use the latest MeasureTheory version.

[cscherrer]

Ok, that's a nice special case. Previously I've had an EqualMix constructor - I should get that into MeaureTheory. We'll have different ways of representing superpositions, so I think it will be something like EqualMixMeasure <: AbstractSuperpositionMeasure

[theogf]

That's currently my work-around which of course involve type piracy:

weight(::MeasureBase.AbstractMeasure) = 1
weight(μ::MeasureBase.WeightedMeasure) = exp(μ.logweight)

function Base.rand(rng::AbstractRNG, μs::AbstractVector{<:MeasureBase.AbstractMeasure})
  ws = weight.(μs)
  ws /= sum(ws) # Normalization of the weights
  μs[rand(rng, Categorical(ws))]
end

function Base.rand(rng::AbstractRNG, ::Type{T}, μ::SuperpositionMeasure) where {T}
  rand(rng, T, rand(rng, μ.components))
end

function Base.rand(rng::AbstractRNG, ::Type{T}, μ::WeightedMeasure) where {T}
  rand(rng, T, μ.base)
end

EDIT: bug fix