 # Uncertainty propagation

Hey,

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)
end
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
end
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)])
end
return r
end
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: GitHub - AnderGray/MomentArithmetic.jl: Rigorous moment propagation with partial information about moments and dependencies in Julia

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