Say I have two variables a and b and I know the full 2x2 covariance between them. Is there a Julia package to conveniently propagate Gaussian uncertainy through arbitrary functions of them in the same way that Measurements.jl works for single variables? Like I want:
(a,b) = multivariate_measurement([10, 12], [1 0.5; 0.5 1]) # mean, covariance
b^2 / a # 14.4 ± 2.1
Or is there some way to do this in Measurements already?
Maybe MonteCarloMeasurements.jl?
That was fast 
julia> a, b = Particles(10000, MvNormal([10,12], [1 0.5; 0.5 1]))
2-element Vector{Particles{Float64, 10000}}:
10.0 ± 1.0
12.0 ± 1.0
julia> b^2 / a
14.4971 ± 2.11 Particles{Float64, 10000}
Beware though, it’s particle-based and not formula-based (at least I think so)
Will keep in mind, I think I’m in an ok regime though. I guess I would be interested still if there’s some AD-like formula-based thing.
If you ever see something like that, let me know. If not, I’d be interested in coding it, in case you wanna help?
Not exactly what you are asking for, but the Reachability Analysis community has some tools for “imprecise probabilities” that may be of value. They are guaranteed to bound the probabilities, however in some cases the over approximation can be too conservative to be useful.
Another thought is that if the function is invertible, computing the push-forward density is simple
Pf(\boldsymbol{x})=f(S^{-1}(\boldsymbol{x}))|\frac{d S^{-1}}{d\boldsymbol{x}}(\boldsymbol{x})|
where S is the function and f is the pdf.
Even if the function is invertible, knowing it is another issue. I don’t know if a multi-variate version of Series Reversion exists (I think I’ve seen it before), but if so, adding that to TaylorSeries.jl would make a pretty general solution, albeit approximate.
You are right, MCM.jl is particle-based, the constructor is even called Particles 
Since you are asking about AD-based solutions, I’ll note that you get exact propagation of multivariate normal distributions through linear functions using the unscented transform. These situations are also the only situations in which AD-based solutions are exact (the ones that use first-order approximations). The unscented transform selects a small number of samples carefully, in a way that allows for exact propagation of Gaussians. In any other case than linear and Gaussian, simple linearization for uncertainty propagation is fraught with a lot of pitfalls, while the unscented transform typically fails slightly less catastrophically.
You can perform the unscented transform using MCM.jl by constructing sigmapoints.
Fyi, the package Covar.jl might be of interest.