[ANN] MonteCarloMeasurements.jl

Package Announcements
, ,
1

As a spinoff from a homework we gave in our julia course at our department, I created a small package to do nonlinear uncertainty propagation, MonteCarloMeasurements.jl.

This package provides two types Particles <: Real and StaticParticles <: Real that represents a distribution of a floating point number, kind of like the type Measurement from Measurements.jl. The difference compared to a Measurement is that Particles represent the distribution using a vector of unweighted particles, and can thus represent arbitrary distributions and handles nonlinear uncertainty propagation well. Functions like f(x) = x² or f(x) = sign(x) at x=0, are not handled well using linear uncertainty propagation ala Measurements.jl.

The goal is to have a number of this type behave just as any other number while partaking in calculations. After a calculation, the mean, std etc. can be extracted from the number using the corresponding functions. Particles also interact with Distributions.jl, so that you can call, e.g., Normal(p) and get back a Normal type from distributions or fit(Gamma, p) to get a Gammadistribution. Particles can also be iterated, asked for maximum/minimum, quantile etc.

The README contains some examples, benchmarks and pretty plots.

Why

Convenience. Also, the benefit of using this number type instead of manually calling a function f with perturbed inputs is that, at least in theory, each intermediate operation on a Particles can exploit SIMD, since it’s performed over a vector. If the function f is called several times, however, the compiler might not be smart enough to SIMD the entire thing. Further, any dynamic dispatch is only paid for once, whereas it would be paid for N times if doing things manually. One could perhaps also make an argument for cache locality being favorable for the Particles type, but I’m not sure this holds for all examples. Primitive benchmarks in the readme show 15x-80x speedup for the considered examples.

The package will be registered once the new registration system is robust, for now
] add https://github.com/baggepinnen/MonteCarloMeasurements.jl.git
is the way to go.

23 Likes
2

This looks like a great package! This gets rid of a lot of ad hoc Monte Carlo boilerplate for me

1 Like
3

Thanks! I hope it works well for you, if not, feel free to open issues as they arise!

1 Like
4

Wow, looks great! Let me know if you’d like MonteCarloMeasurements.jl to be hosted under JuliaPhysics :wink:

2 Likes
5

Interesting, do I understand correctly that one could work out a weighted particle-cloud type which essentially acts as particle filter which incorporates a resampling step if operations involve several particle clouds?

6

Great package! If I understand correctly, you are assuming different data points to be statistically independent (since you simply forward Statictics.std etc.). Maybe it’d make sense to use/integrate some features of BinningAnalysis.jl. Just a thought on my side.

1 Like
7

You could certainly do that and I was contemplating the same. This only becomes interesting if there is something to weight them by, such as a repeated measurement etc. I do however already have a library for particle filters and didn’t find the need to experiment with this now ^^

8

Different Particles are independent if you create them independently. If you create them with the constructor Particles(N, d::MultivariateDistribution), they will be correlated as one would expect based on the distribution.

9

I would have the hope that generalised Particles could provide uncertainty quantification of generic code in a similar way as Dual numbers provides derivatives… ping me if you ever make an attempt!

1 Like
10

This package actually does exactly that, but not in they way you envision. The unweighted particle cloud encodes a distribution by varying its density instead of weighting the particles. A particle filter does both, the dynamics noise alter the density of the particle cloud, whereas the measurement density function provides weights. In this setting, there are no measurements to weight by so the weighting is never carried out.

1 Like
11

The alternative approach would be to sample particles uniformly and then weight them by the distribution, this would be suboptimal as equally many particles would be used for low-density areas as for high-density areas.

12

Yes, right, I was somehow jumping ahead with the measurements, which would not turn up in generic code. But it would be something one could want to include anyway, as there will be already resampling steps because of operations having several independent particle clouds as input (that is not yet included in your code, right?)

13

Certainly, it would be kind of cool to have a particle filter implemented using a particle type! I have no facilities for resampling of the particles at the moment since the weighting required is not a thing for uncertainty propagation, as it is for bayesian filtering (particle filtering). Having said that, the extension required would be a new type WeigtedParticles <: AbstractParticles with some additional functions for resampling as well as normalization, logsumexp etc. That would be a project for another time, I’m afraid, but a cool idea nontheless.

3 Likes
14

Thinking a bit more along those lines… This type could be used for black-box optimization with particle methods, e.g., particle-swarm optimization. The weighting would then be provided by the (inverse) cost function, and resampling would focus on promising areas of the parameter space. Maybe it would be a good idea to extract the particle type into a separate package to facilitate all these crazy ideas :stuck_out_tongue:

2 Likes
15

Not the most urgent thing to decide but Swarm.jl or ParticleSwarm.jl sounds cool.

16

ParticleSwarms.jl is clearer and reflects use in the literature (edited to pick up next comment). There are some repos named ParticleSwarmOptimization.jl

17

At the risk of exacerbating the way too early naming…if I understand the goal of the proposal correctly, I think ParticleSwarms.jl would align with the plural convention for things like Distributions.jl and StaticArrays.jl (if each subtype is a type of ParticleSwarm)

2 Likes
18

The package is now registrated and can be added with

] add MonteCarloMeasurements

(make sure to ] update first)

The readme is extended with further documentation, more plotting functions are implemented and I also added an example file where I visualize the Bode and Nyquist diagrams of an uncertain transfer-function model, as well as perform a time-domain simulation of the system and look at the uncertainties in the response.

I verified that the package works with DifferentialEquations.jl. The tutorial for solving differential equations using Measurement works for Particles as well. In the second pendulum example, I had to call solve(..., dt=0.01, adaptive=false) , otherwise the solver tried to estimate a suitable dt using the uncertain state, which created some problems.

On a side note, I implemented a very primitive CMA optimization algorithm as a proof of concept. It’s available here and seems to work okay. I do not intend to take this further though, as there is already the well developed BlackBoxOptim.jl

1 Like
19

Updates in v0.1.4

I implemented the weighted particles with resampling functionality (using numerically stable logsumexp).
I also added two macros, @bymap, @bypmap, detailed below.

Weighted particles

The type WeightedParticles contains an additional field logweights. You may modify this field as you see fit, e.g.

reweight(p,y) = (p.logweights .+= logpdf.(Normal(0,1), y .- p.particles))

where y would be some measurement. After this you can resample the particles using resample!(p). This performs a systematic resample with replacement, where each particle is sampled proportionally to exp.(logweights).

Monte-Carlo simulation by map/pmap

Some functions will not work when the input arguments are of type Particles. For this kind of function, we provide a fallback onto a traditional map(f,p.particles). The only thing you need to do is to decorate the function call with the macro @bymap like so:

f(x) = 3x^2
p = 1 ± 0.1
r = @bymap f(p)

We further provide the macro @bypmap which does exactly the same thing, but with a pmap (parallel map) instead, allowing you to run several invocations of f in a distributed fashion.

These macros will map the function f over each element of p::Particles{T,N}, such that f is only called with arguments of type T, e.g., Float64. This handles arguments that are multivaiate particles <: Vector{<:AbstractParticles} as well.

These macros will typically be slower than calling f(p), but allow for Monte-Carlo simulation of things like calls to Optim.optimize etc., which fail if called like optimize(f,p). If f is very expensive, @bypmap might prove prove faster than calling f with p, it’s worth a try. The usual caveats for distributed computing applies, all code must be loaded on all workers etc.

4 Likes
20

Add an ODE_DEFAULT_NORM overload in the init file of DiffEqBase so that the error norm returns a float if time is a float. That will make it work more nicely with adaptive time stepping. Otherwise time itself would need to be a swarm