# [ANN] MonteCarloMeasurements.jl

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 `Gamma`distribution. 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

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

1 Like

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

1 Like

Wow, looks great! Let me know if youâ€™d like `MonteCarloMeasurements.jl` to be hosted under `JuliaPhysics`

2 Likes

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?

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

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

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.

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

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

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.

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

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

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

2 Likes

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

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

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

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

# 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

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