Exported functions and types

This package facilitates working with probability distributions by means of Monte-Carlo methods, in a way that allows for propagation of probability distributions through functions. This is useful for, e.g., nonlinear uncertainty propagation. A variable or parameter might be associated with uncertainty if it is measured or otherwise estimated from data. We provide two core types to represent probability distributions: Particles and StaticParticles, both <: Real. (The name "Particles" comes from the particle-filtering literature.) These types all form a Monte-Carlo approximation of the distribution of a floating point number, i.e., the distribution is represented by samples/particles. Correlated quantities are handled as well, see multivariate particles below.

A number of type Particles behaves just as any other Number while partaking in calculations. Particles also behave like a distribution, so after a calculation, an approximation to the complete distribution of the output is captured and represented by the output particles. mean, std etc. can be extracted from the particles using the corresponding functions pmean and pstd. 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 asked for maximum/minimum, quantile etc. using functions with a prefix p, i.e., pmaximum. If particles are plotted with plot(p), a histogram is displayed. This requires Plots.jl. A kernel-density estimate can be obtained by density(p) is StatsPlots.jl is loaded.

Quick start

julia> using MonteCarloMeasurements, Plots

julia> a = π ± 0.1 # Construct Gaussian uncertain parameters using ± (\pm)
Particles{Float64,2000}
 3.14159 ± 0.1

julia> b = 2 ∓ 0.1 # ∓ (\mp) creates StaticParticles (with StaticArrays)
StaticParticles{Float64,100}
 2.0 ± 0.0999

julia> pstd(a)      # Ask about statistical properties
0.09999231528930486

julia> sin(a)      # Use them like any real number
Particles{Float64,2000}
 1.2168e-16 ± 0.0995

julia> plot(a)     # Plot them

julia> b = sin.(1:0.1:5) .± 0.1; # Create multivariate uncertain numbers

julia> plot(b)                   # Vectors of particles can be plotted

julia> using Distributions

julia> c = Particles(500, Poisson(3.)) # Create uncertain numbers distributed according to a given distribution
Particles{Int64,500}
 2.882 ± 1.7

For further help, see the documentation, the examples folder or the arXiv paper.

struct Particles{T, N} <: AbstractParticles{T, N}

This type represents uncertainty using a cloud of particles.

Constructors:

• Particles()
• Particles(N::Integer)
• Particles([rng::AbstractRNG,] d::Distribution)
• Particles([rng::AbstractRNG,] N::Integer, d::Distribution; permute=true, systematic=true)
• Particles(v::Vector{T} where T)
• Particles(m::Matrix{T} where T): Creates multivariate particles (Vector{Particles})

struct StaticParticles{T, N} <: AbstractParticles{T, N}

See ?Particles for help. The difference between StaticParticles and Particles is that the StaticParticles store particles in a static vecetor. This makes runtimes much shorter, but compile times longer. See the documentation for some benchmarks. Only recommended for sample sizes of ≲ 300-400

struct Workspace{T1, T2, T3, T4, T5, T6}

Arguments:

• simple_input: Input object f will be called with, does not contain any particles
• simple_result: Simple output from f without particles
• result: Complete output of f including particles
• buffersetter: Helper function to shift data between objects
• resultsetter: Helper function to shift data between objects
• f: Function to call
• N: Number of particles

Workspace(f, input)

Create a Workspace object for inputs of type typeof(input). Useful if input is a structure with fields of type <: AbstractParticles (can be deeply nested). See also with_workspace.

a .. b

Creates 2000 Particles with a Uniform distribution between a and b. See also ±, ⊗

μ ± σ

Creates 2000 Particles with mean μ and std σ. It can also be used as a unary operator, a mean of 0 is then used with std σ. If μ is a vector, the constructor MvNormal is used, and σ is thus treated as std if it's a scalar, and variances if it's a matrix or vector. See also ∓, ..

μ ∓ σ

Creates 100 StaticParticles with mean μ and std σ. It can also be used as a unary operator, a mean of 0 is then used with std σ. If μ is a vector, the constructor MvNormal is used, and σ is thus treated as std if it's a scalar, and variances if it's a matrix or vector. See also ±, ⊗

a ⊞ Distribution()

Adds 2000 Particles sampled from a specified ::Distribution to a. Shorthand for a + Particles(Distribution()), e.g., 1 ⊞ Binomial(3).

a ⊠ Distribution()

Multiplies a by 2000 Particles sampled from a specified ::Distribution. Shorthand for a * Particles(Distribution()), e.g., a ⊠ Gamma(1).

bootstrap([rng::AbstractRNG,] p::Particles, n = nparticles(p))

Return Particles resampled with replacement. n specifies the number of samples to draw. Also works for arrays of Particles, in which case a single set of indices are drawn and used to extract samples from all elements in the array.

bootstrap([rng::AbstractRNG,] p::Particles, n = nparticles(p))

Return Particles resampled with replacement. n specifies the number of samples to draw. Also works for arrays of Particles, in which case a single set of indices are drawn and used to extract samples from all elements in the array.

bymap(f, args...)

Uncertainty propagation using the map function.

Call f with particles or vectors of particles by using map. This can be utilized if registering f using register_primitive fails. See also Workspace if bymap fails.

Distributed uncertainty propagation using the pmap function. See bymap for more details.

errorbarplot(x,y,[q=0.025])

Plots a vector of particles with error bars at quantile q. If q::Tuple, then you can specify both lower and upper quantile, e.g., (0.01, 0.99).

ess(p::AbstractParticles{T,N})

Calculates the effective sample size. This is useful if particles come from MCMC sampling and are correlated in time. The ESS is a number between [0,N].

Initial source: https://github.com/tpapp/MCMCDiagnostics.jl

has_particles(P)

Determine whether or no the object P has some kind of particles inside it. This function examins fields of P recursively and looks inside arrays etc.

mcplot(x,y,[N=0])

Plots all trajectories represented by a vector of particles. N > 1 controls the number of trajectories to plot.

mean_object(x)

Returns an object similar to x, but where all internal instances of Particles are replaced with their mean. The generalization of this function is replace_particles.

nominal(p)

Return the nominal value of p (assumes that p has been endowed with a nominal value using with_nominal).

p = outer_product([rng::AbstractRNG,] dists::Vector{<:Distribution}, N=100_000)

Creates a multivariate systematic sample where each dimension is sampled according to the corresponding univariate distribution in dists. Returns p::Vector{Particles} where each Particles has a length approximately equal to N. The particles form the outer product between d systematically sampled vectors with length given by the d:th root of N, where d is the length of dists, All particles will be independent and have marginal distributions given by dists.

See also MonteCarloMeasurements.⊗

register_primitive(f, eval=eval)

Register both single and multi-argument function so that it works with particles. If you want to register functions from within a module, you must pass the modules eval function.

register_primitive_multi(ff, eval=eval)

Register a multi-argument function so that it works with particles. If you want to register functions from within a module, you must pass the modules eval function.

register_primitive_single(ff, eval=eval)

Register a single-argument function so that it works with particles. If you want to register functions from within a module, you must pass the modules eval function.

ribbonplot(x,y,[q=0.025]; N=true)

Plots a vector of particles with a ribbon covering quantiles q, 1-q. If q::Tuple, then you can specify both lower and upper quantile, e.g., (0.01, 0.99).

If a positive number N is provided, N sample trajectories will be plotted on top of the ribbon.

set_comparison_function(f)

Change the Function used to reduce particles to a number for comparison operators Toggle the use of a comparison Function without warning using the Function unsafe_comparisons.

sigmapoints(m, Σ)
sigmapoints(d::Normal)
sigmapoints(d::MvNormal)

The unscented transform uses a small number of points to propagate the first and second moments of a probability density, called sigma points. We provide a function sigmapoints(μ, Σ) that creates a Matrix of 2n+1 sigma points, where n is the dimension. This can be used to initialize any kind of AbstractParticles, e.g.:

julia> m = [1,2]

julia> Σ = [3. 1; 1 4]

julia> p = StaticParticles(sigmapoints(m,Σ))
2-element Array{StaticParticles{Float64,5},1}:
 (5 StaticParticles: 1.0 ± 1.73)
 (5 StaticParticles: 2.0 ± 2.0)

julia> cov(p) ≈ Σ
true

julia> mean(p) ≈ m
true

Make sure to pass the variance (not std) as second argument in case μ and Σ are scalars.

Caveat

If you are creating several one-dimensional uncertain values using sigmapoints independently, they will be strongly correlated. Use the multidimensional constructor! Example:

p = StaticParticles(sigmapoints(1, 0.1^2))               # Wrong!
ζ = StaticParticles(sigmapoints(0.3, 0.1^2))             # Wrong!
ω = StaticParticles(sigmapoints(1, 0.1^2))               # Wrong!

p,ζ,ω = StaticParticles(sigmapoints([1, 0.3, 1], 0.1^2)) # Correct

systematic_sample([rng::AbstractRNG,] N, d=Normal(0,1); permute=true)

returns a Vector of length N sampled systematically from the distribution d. If permute=false, this vector will be sorted.

Y = transform_moments(X::Matrix, m, Σ; preserve_latin=false)

Transforms X such that it get the specified mean and covariance.

m, Σ   = [1,2], [2 1; 1 4] # Desired mean and covariance
particles = transform_moments(X, m, Σ)
julia> cov(particles) ≈ Σ
true

Note, if X is a latin hypercube and Σ is non-diagonal, then the latin property is destroyed for all dimensions but the first. We provide a method preserve_latin=true) which absolutely preserves the latin property in all dimensions, but if you use this, the covariance of the sample will be slightly wrong

unsafe_comparisons(onoff=true; verbose=true)

Toggle the use of a comparison function without warning. By default mean is used to reduce particles to a floating point number for comparisons. This function can be changed, example: set_comparison_function(median)

unsafe_comparisons(mode=:reduction; verbose=true)

One can also specify a comparison mode, mode can take the values :safe, :montecarlo, :reduction. :safe is the same as calling unsafe_comparisons(false) and :reduction corresponds to true. If

wasserstein(p1::AbstractParticles,p2::AbstractParticles,p)

Returns the Wasserstein distance (Earth-movers distance) of order p, to the pth power, between p1 and p2. I.e., for p=2, this returns W₂²

pn = with_nominal(p, val)

Endow particles p with a nominal value val. The particle closest to val will be replaced with val, and moved to index 1. This operation introduces a slight bias in the statistics of pn, but the operation is asymptotically unbiased for large sample sizes. To obtain the nominal value of pn, call nominal(pn).

with_workspace(f,P)

In some cases, defining a primitive function which particles are to be propagate through is not possible but allowing unsafe comparisons are not acceptable. One such case is functions that internally calculate eigenvalues of uncertain matrices. The eigenvalue calculation makes use of comparison operators. If the uncertainty is large, eigenvalues might change place in the sorted list of returned eigenvalues, completely ruining downstream computations. For this we recommend, in order of preference

1. Use @bymap detailed in the documentation. Applicable if all uncertain values appears as arguments to your entry function.
2. Create a Workspace object and call it using your entry function. Applicable if uncertain parameters appear nested in an object that is an argument to your entry function:

# desired computation: y = f(obj), obj contains uncertain parameters inside
y = with_workspace(f, obj)
# or equivalently
w = Workspace(f, obj)
use_invokelatest = true # Set this to false to gain 0.1-1 ms, at the expense of world-age problems if w is created and used in the same function.
w(obj, use_invokelatest)

ℂ2ℂ_function(f::Function, z::Complex{<:AbstractParticles})

Helper function for uncertainty propagation through complex-valued functions of complex arguments. applies f : ℂ → ℂ to z::Complex{<:AbstractParticles}.

ℝⁿ2ℂⁿ_function(f::Function, p::AbstractArray{T})

Helper function for performing uncertainty propagation through complex-valued functions with vector inputs. Applies f : ℝⁿ → Cⁿ to an array of particles. E.g., LinearAlgebra.eigvals(p::Matrix{<:AbstractParticles}) = ℝⁿ2ℂⁿ_function(eigvals,p)

ℝⁿ2ℂⁿ_function(f::Function, p::AbstractArray{T})

Helper function for performing uncertainty propagation through complex-valued functions with vector inputs. Applies f : ℝⁿ → Cⁿ to an array of particles. E.g., LinearAlgebra.eigvals(p::Matrix{<:AbstractParticles}) = ℝⁿ2ℂⁿ_function(eigvals,p)

ℝⁿ2ℝⁿ_function(f::Function, p::AbstractArray{T})

Helper function for performing uncertainty propagation through vector-valued functions with vector inputs. Applies f : ℝⁿ → ℝⁿ to an array of particles. E.g., Base.log(p::Matrix{<:AbstractParticles}) = ℝⁿ2ℝⁿ_function(log,p)

ℝⁿ2ℝⁿ_function(f::Function, p::AbstractArray{T})

Helper function for performing uncertainty propagation through vector-valued functions with vector inputs. Applies f : ℝⁿ → ℝⁿ to an array of particles. E.g., Base.log(p::Matrix{<:AbstractParticles}) = ℝⁿ2ℝⁿ_function(log,p)

@bymap f(p, args...)

Call f with particles or vectors of particles by using map. This can be utilized if registering f using register_primitive fails. See also Workspace if bymap fails.

@bypmap f(p, args...)

Call f with particles or vectors of particles by using parallel pmap. This can be utilized if registering f using register_primitive fails. See also Workspace if bymap fails.

@prob a < b

Calculate the probability that an event on any of the forms a < b, a > b, a <= b, a >= b occurs, where a and/or b are of type AbstractParticles.

@unsafe expression

Activates unsafe comparisons for the provided expression only. The expression is surrounded by a try/catch block to robustly restore unsafe comparisons in case of exception.

p1 ≈ p2

Determine if two particles are not significantly different

p1 ≉ p2

Determine if two particles are significantly different