Uncertainty propagation


Consider a random dataset:

using Distributions
using Random
using Statistics
using StatsBase
X = rand(LogNormal(0,1),10000)

Consider a first function that ‘summarises’ the dataset into a simple vector of k statistics (here k=16 for this simple example):

function moms(x)
    return mean(x .^ (0:15)' .* exp.(-x),dims=1)
m = moms(X)

and a second function that does a lot of heavy computations on these statistics, which append to also output a k-sized vector :

function heavy_computations(m)
    # do a lot a stuff with m 
    result = log.(1 .+abs.(m)) .+ 1
    return result 
r = heavy_computations(m)

Now i want to asses the empirical covariance of the final k statistics. I could bootstrap:

function boot(M,data,f)
    n_obs = length(data)
    r = zeros((M,length(f(data))))
    for i in 1:M
        r[i,:] = f(data[sample(1:n_obs,n_obs,replace=true)])
    return r
M = 1000
cov_r = cov(boot(M,X,heavy_computations ∘ moms))

but it computes the heavy_computations function M times ! It’s on the other hand way faster to only compute the covariance of m:

cov_m = cov(boot(1000,X,moms))

Could I leverage the fact that the heavy_computations function is written in Julia, a little like Measurement.jl but multivariate ?

I saw this package the other day for moment propagation, it may be related to your goals: https://github.com/AnderGray/MomentArithmetic.jl

Unfortunately, this is not exactly what I want, but thanks for the catch, really interesting stuff.

What I want is a mechanisme to propagate a second order information through a piece of code. I’m not even sure it’s possible…

GitHub - baggepinnen/MonteCarloMeasurements.jl: Propagation of distributions by Monte-Carlo sampling: Real number types with uncertainty represented by samples. does this, it uses samples, but it’s not exactly the same as running heavy_comoutation M times due to the way its implemented. There are some details on this in the paper linked in the Readme.

1 Like