# \[ANN\] MonteCarloMeasurements.jl

**URL:** <https://discourse.julialang.org/t/ann-montecarlomeasurements-jl/23021>\
**Category:** Package Announcements\
**Tags:** announcement, monte-carlo, measurements\
**Created:** [April 11, 2019, 7:21am UTC](https://discourse.julialang.org/t/ann-montecarlomeasurements-jl/23021 "2019-04-11T07:21:14Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [April 11, 2019, 7:21am UTC](https://discourse.julialang.org/t/ann-montecarlomeasurements-jl/23021/1 "2019-04-11T07:21:14Z")

</div>

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](https://github.com/baggepinnen/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](https://github.com/JuliaPhysics/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](https://github.com/JuliaPhysics/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](https://github.com/baggepinnen/MonteCarloMeasurements.jl) contains some examples, benchmarks and [pretty plots](https://github.com/baggepinnen/MonteCarloMeasurements.jl#plotting).

## 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.

---

<div class="post-metadata">

**Author:** ![ajkeith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ajkeith/32/970_2.png) [@ajkeith](https://discourse.julialang.org/u/ajkeith)\
**Post date:** [April 11, 2019, 10:49am UTC](https://discourse.julialang.org/t/ann-montecarlomeasurements-jl/23021/2 "2019-04-11T10:49:01Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [April 11, 2019, 10:58am UTC](https://discourse.julialang.org/t/ann-montecarlomeasurements-jl/23021/3 "2019-04-11T10:58:25Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![giordano](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/giordano/32/2166_2.png) [@giordano](https://discourse.julialang.org/u/giordano)\
**Post date:** [April 11, 2019, 11:40am UTC](https://discourse.julialang.org/t/ann-montecarlomeasurements-jl/23021/4 "2019-04-11T11:40:43Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![mschauer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mschauer/32/13946_2.png) [@mschauer](https://discourse.julialang.org/u/mschauer)\
**Post date:** [April 11, 2019, 12:11pm UTC](https://discourse.julialang.org/t/ann-montecarlomeasurements-jl/23021/5 "2019-04-11T12:11:08Z")

</div>

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?

---

<div class="post-metadata">

**Author:** ![carstenbauer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/carstenbauer/32/4981_2.png) [@carstenbauer](https://discourse.julialang.org/u/carstenbauer)\
**Post date:** [April 11, 2019, 12:13pm UTC](https://discourse.julialang.org/t/ann-montecarlomeasurements-jl/23021/6 "2019-04-11T12:13:42Z")

</div>

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](https://github.com/crstnbr/BinningAnalysis.jl). Just a thought on my side.

---

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [April 11, 2019, 12:22pm UTC](https://discourse.julialang.org/t/ann-montecarlomeasurements-jl/23021/7 "2019-04-11T12:22:53Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [April 11, 2019, 12:24pm UTC](https://discourse.julialang.org/t/ann-montecarlomeasurements-jl/23021/8 "2019-04-11T12:24:09Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![mschauer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mschauer/32/13946_2.png) [@mschauer](https://discourse.julialang.org/u/mschauer)\
**Post date:** [April 11, 2019, 12:26pm UTC](https://discourse.julialang.org/t/ann-montecarlomeasurements-jl/23021/9 "2019-04-11T12:26:52Z")

</div>

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!

---

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [April 11, 2019, 12:32pm UTC](https://discourse.julialang.org/t/ann-montecarlomeasurements-jl/23021/10 "2019-04-11T12:32:01Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [April 11, 2019, 12:33pm UTC](https://discourse.julialang.org/t/ann-montecarlomeasurements-jl/23021/11 "2019-04-11T12:33:48Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![mschauer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mschauer/32/13946_2.png) [@mschauer](https://discourse.julialang.org/u/mschauer)\
**Post date:** [April 11, 2019, 12:36pm UTC](https://discourse.julialang.org/t/ann-montecarlomeasurements-jl/23021/12 "2019-04-11T12:36:51Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [April 11, 2019, 12:47pm UTC](https://discourse.julialang.org/t/ann-montecarlomeasurements-jl/23021/13 "2019-04-11T12:47:28Z")

</div>

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](https://github.com/baggepinnen/LowLevelParticleFilters.jl/blob/master/src/resample.jl) for resampling as well as normalization, logsumexp etc. That would be a project for another time, I’m afraid, but a cool idea nontheless.

---

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [April 12, 2019, 1:01pm UTC](https://discourse.julialang.org/t/ann-montecarlomeasurements-jl/23021/14 "2019-04-12T13:01:18Z")

</div>

Thinking a bit more along those lines… This type could be used for black-box optimization with particle methods, e.g., [particle-swarm optimization](https://en.wikipedia.org/wiki/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 😛

---

<div class="post-metadata">

**Author:** ![mschauer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mschauer/32/13946_2.png) [@mschauer](https://discourse.julialang.org/u/mschauer)\
**Post date:** [April 12, 2019, 5:08pm UTC](https://discourse.julialang.org/t/ann-montecarlomeasurements-jl/23021/15 "2019-04-12T17:08:51Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![JeffreySarnoff](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jeffreysarnoff/32/1980_2.png) [@JeffreySarnoff](https://discourse.julialang.org/u/JeffreySarnoff)\
**Post date:** [April 12, 2019, 5:10pm UTC](https://discourse.julialang.org/t/ann-montecarlomeasurements-jl/23021/16 "2019-04-12T17:10:58Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![ajkeith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ajkeith/32/970_2.png) [@ajkeith](https://discourse.julialang.org/u/ajkeith)\
**Post date:** [April 14, 2019, 3:52pm UTC](https://discourse.julialang.org/t/ann-montecarlomeasurements-jl/23021/17 "2019-04-14T15:52:12Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [April 16, 2019, 3:24pm UTC](https://discourse.julialang.org/t/ann-montecarlomeasurements-jl/23021/18 "2019-04-16T15:24:10Z")

</div>

The package is now registrated and can be added with

```julia
] 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](https://github.com/baggepinnen/MonteCarloMeasurements.jl/blob/master/examples/controlsystems.jl) 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](http://juliadiffeq.org/DiffEqTutorials.jl/html/type_handling/uncertainties.html) 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](https://en.wikipedia.org/wiki/CMA-ES) optimization algorithm as a proof of concept. It’s available [here](https://github.com/baggepinnen/MonteCarloMeasurements.jl/blob/master/src/optimize.jl) and seems to work okay. I do not intend to take this further though, as there is already the well developed [BlackBoxOptim.jl](https://github.com/robertfeldt/BlackBoxOptim.jl)

---

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [April 24, 2019, 12:17pm UTC](https://discourse.julialang.org/t/ann-montecarlomeasurements-jl/23021/19 "2019-04-24T12:17:05Z")

</div>

# Updates in v0.1.4

I implemented the weighted particles with [resampling functionality](https://github.com/baggepinnen/MonteCarloMeasurements.jl/blob/master/src/resampling.jl) (using [numerically stable logsumexp](https://discourse.julialang.org/t/fast-logsumexp/22827/7)).  
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.

```julia
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:

```julia
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.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [April 25, 2019, 4:07am UTC](https://discourse.julialang.org/t/ann-montecarlomeasurements-jl/23021/20 "2019-04-25T04:07:28Z")

</div>

> [@baggepinnen](#):
>
> 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

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

[Next page](https://discourse.julialang.org/t/ann-montecarlomeasurements-jl/23021.md?page=2)