Several non-exported functions that may facilitate working with structures that contain uncertain parameters (struct-of-arrays, SoA) exist. These are not to be considered part of the API and are subject to breakage at any time, but may nevertheless be of use in special situations.
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.
replace_particles(x; condition=P->P isa AbstractParticles,replacer = P->vecindex(P, 1))
This function recursively scans through the structure x, every time a field that matches condition is found, replacer is called on that field and the result is used instead of P. See function mean_object, which uses this function to replace all instances of Particles with their mean.
Exectues f on each instance of arg represented by internal particles of arg. This is useful as a last resort if all other methods to propagate particles through f fails. The function returns an array (length = num. particles) of structs rather than particles, each struct is the result of f(replace_particles(arg, p->p[i])).
Several non-exported functions that may facilitate working with structures that contain uncertain parameters (struct-of-arrays, SoA) exist. These are not to be considered part of the API and are subject to breakage at any time, but may nevertheless be of use in special situations.
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.
replace_particles(x; condition=P->P isa AbstractParticles,replacer = P->vecindex(P, 1))
This function recursively scans through the structure x, every time a field that matches condition is found, replacer is called on that field and the result is used instead of P. See function mean_object, which uses this function to replace all instances of Particles with their mean.
Exectues f on each instance of arg represented by internal particles of arg. This is useful as a last resort if all other methods to propagate particles through f fails. The function returns an array (length = num. particles) of structs rather than particles, each struct is the result of f(replace_particles(arg, p->p[i])).
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
+API · MonteCarloMeasurements Documentation
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}
@@ -26,7 +26,7 @@
julia> c = Particles(500, Poisson(3.)) # Create uncertain numbers distributed according to a given distribution
Particles{Int64,500}
- 2.882 ± 1.7
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
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.
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 ±, ⊗
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.
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.
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].
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.
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.
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 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 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.
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.
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.
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
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.
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 ±, ⊗
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.
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.
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].
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.
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.
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 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 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.
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.
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.
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.:
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
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.
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).
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
Use @bymap detailed in the documentation. Applicable if all uncertain values appears as arguments to your entry function.
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:
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
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.
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).
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
Use @bymap detailed in the documentation. Applicable if all uncertain values appears as arguments to your entry function.
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)
Helper function for uncertainty propagation through complex-valued functions of complex arguments. applies f : ℂ → ℂ to z::Complex{<:AbstractParticles}.
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)
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)
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)
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)
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.
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.
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.
Helper function for uncertainty propagation through complex-valued functions of complex arguments. applies f : ℂ → ℂ to z::Complex{<:AbstractParticles}.
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)
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)
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)
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)
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.
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.
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.
as we can see, the linear method outputs a Dirac distribution (no uncertainty) at $x=0$, while there should clearly be a lot of uncertainty in the output. The histogram displays the output density as approximated by the particles. The histogram does not go below zero, and tapers off as values increase. The problem here is that the uncertainty is large in relation to the curvature of the function. As the uncertainty decreases, the true output density becomes closer and closer to a DIrac distribution.
The next function has a discontinuity (≈ infinite curvature)
once again, linear uncertainty propagation outputs a distribution with zero uncertainty. The true output distribution has two modes since the input distribution places mass on both sides of the discontinuity. This is captured in the particle distribution, where most particles end up at the right side of the discontinuity, while a smaller proportion of the particles end up to the left. If the input density would have its mean at 0, half of the particles would end up in each of the output locations. Any function containing an if-statement where the chosen branch depends on an uncertain value falls into this category.
Once again, the uncertainty is large in relation to the curvature of the function and linear uncertainty propagation places significant mass outside the interval $[-1, 1]$ which is the range of the $\sin$ function. The particle histogram respects this range. If we increase the uncertainty in the input further, the linear approximation to the function becomes increasingly worse