Rework Weighted Particles

[cscherrer]

I've been thinking about about the approach for WeightedParticles. When I've needed this I haven't been able to use the built-in type. As it stands, the weighting vector is per-Particles, but I've needed a shared weighting vector across several Particles.

Is the non-sharing an intentional design choice, or are you open to this sharing idea? I think it would look something like this:

struct LogWeights{T,N} <: AbstractParticles{T,N}
    particles::Vector{T}
    next 
end

struct WeightedParticles{T,N} <: AbstractParticles{T,N}
    particles::Vector{T}
    logweights::Ref{LogWeights{W,N} where W <: Real}
end

Functions like +, *, etc on WeightedParticles would require === on the logweights, to be sure that the entire system evolves as one. The next function is mostly for resampling. The idea is that when a resample is triggered, existing particles become "stale" and need to pass through this function to reindex the values.

I think there are a few ways of implementing this last bit, but I wanted to get your take on the shared weights idea first before digging into it much more.

[baggepinnen]

I am quite unhappy with how weighted particles work at the moment and was actually considering removing support for them, keeping only the weight-manipulating functions. You have probably noticed how the use of weighted particles causes quite a lot of warnings of poor support...

Certainly, the current implementation carries around a lot of redundant weights and some form of weight sharing is absolutely required. Your proposed approach above seems interesting. I can't really visualize now how the details would look, but it could work. Possibly, the particles field in

struct WeightedParticles{T,N} <: AbstractParticles{T,N}
    particles::Vector{T}

could be of type Particles{T,N} to reuse the functionality for manipulating the particles themselves, with weight manipulations added on outside.

[cscherrer]

In that case, maybe it would make sense to remove them from master and add a branch to prototype an alternative? I've got a sprint coming up for the JuliaCon proceedings deadline, but after a couple of weeks things will clear back up, and I'd love to help get something more robust in place for this

[baggepinnen]

There's now a PR that removes support for WeightedParticles until they are better supported
#33

[cscherrer]

Hi @baggepinnen , I was thinking some more a bout this today... I think we could do something like

struct Weights{T,N}
    vals::Vector{T}
    members :: MutableLinkedList{WeightedParticles{T,N}}
end


struct WeightedParticles{T,N} <: AbstractParticles{T,N}
    particles::Vector{T}
    weights::Weights{T}

    function WeightedParticles(x,w)
        x = Particles(x).particles
        N = length(x)
        w = Weights(zeros(N), MutableLinkedList{WeightedParticles{T,N}}())
        wp = WeightedParticles{eltype(x),N}(x,w)
        push!(w.members, wp)
    end
end

members links weights to all the particles using those weights, important for resampling. And weights can be === to save on space and computation. So e.g.,

julia> x = WeightedParticles(Particles().particles, ones(2000))
WeightedParticles{Float64,Array{Float64,1},2000}
 1.77636e-18 ± 1.0

julia> y = WeightedParticles(Particles(), x.weights)
WeightedParticles{Float64,Array{Float64,1},2000}
 1.24345e-17 ± 1.0

julia> x.weights[1] = 2.0
2.0

julia> y.weights[1]
2.0

All the stats need to be reworked for the weighted case, but hopefully this gets the idea across. And I guess you might have some old code for all of those.

[baggepinnen]

Hi Chad, your proposed approach looks much better than what I had before. Let's see if I have understood it properly

The members field has a list of all WeightedParticles that are part of the same "multivariate" particles, i.e., Vector{WeightedParticles}? Each element in such a vector in turn has a reference to the same Weights. An operation between two WeightedParticles creates a new WeightedParticles with new ´Weights`?

I did have some functionality in place to operate on weighted particles, like *, std, mean etc.
I did not really figure out a good way to implement +,- between WeightedParticles that did not have the same Weights, maybe it's trivial but I just couldn't easily see how to?

[cscherrer]

The members field has a list of all WeightedParticles that are part of the same "multivariate" particles, i.e., Vector{WeightedParticles}? Each element in such a vector in turn has a reference to the same Weights.

Right. I'd think the multivariate-ness of it would usually be implicit, the same way Particles are now. I assume that's what you mean, but it's good to confirm we're thinking of this in the same way.

An operation between two WeightedParticles creates a new WeightedParticles with new ´Weights`?

I'd think it would be just like Particles, that operations would be elementwise. So by default, operations would require the weights be identical. But we really need to be careful here, and maybe you've worked through the details more than I have. Do elementwise operations give the right distributional approximation when the weights are nonuniform? I'll think about this some more, may have missed something obvious

[cscherrer]

For a moment I thought this might be entirely wrong, but I'm pretty sure this is correct. In particle filtering I've done (maybe there are variations?) each particle maintains a state of the "world", together with a weight (usually as a log). With Particles, each world is an index that's used across particles. But the weights are still global, so it seems sensible to me to optimize for the case of one "global" vector of weights shared across all particles. Ideally, the weight vector could even be encoded in the type, maybe initiated as a Val(gensym()) field or something.

Maybe it would be helpful to work through an importance sampling example, like sampling from a Student's T and weighting it to act like a Gaussian, then going through some elementwise operations. Since we know the right answers, this could also be a good basis for tests.

[cscherrer]

Ok here's an example. I think all of this checks out:

using MonteCarloMeasurements
using Distributions

# Expected value of x log-weighted by ℓ
function expect(x,ℓ)
    w = exp(ℓ - maximum(ℓ))
    mean(x*w) / mean(w)
end

x = Particles(5000,TDist(3))
ℓ = logpdf(Normal(),x) - logpdf(TDist(3), x)

expect(x,ℓ)
expect(x^2,ℓ)
expect(x^3,ℓ)
expect(x^4,ℓ)

# Now sample z ~ Normal()

z = Particles(5000,Normal())

expect(z,ℓ)
expect(z^2,ℓ)
expect(z^3,ℓ)
expect(z^4,ℓ)

# And now a new weighted variable

y = Particles(5000,TDist(3))
ℓ += logpdf(Normal(),y) - logpdf(TDist(3), y)

expect(y,ℓ)
expect(y^2,ℓ)
expect(y^3,ℓ)
expect(y^4,ℓ)

expect(x+y,ℓ)
expect((x+y)^2,ℓ)

expect(x+z,ℓ)
expect((x+z)^2,ℓ)