We expose a problem of estimation of multivariate convolutions of gamma random variables, which has very bad numerical properties. This bad numerical behavior literally forced us to use Julia. We describe why Python, R or C++ were not capable of solving our problem and argue that the multiple dispatch paradigm in Julia was the reason we were able to reuse existing code.
The estimation of multivariate generalized Gamma convolutions via their projections into the Laguerre basis is easy to deal with mathematically. However, this deconvolution problem gives a loss function that requires the optimization routine to be:
In both Python, R and C++, these three conditions were incompatible and forces us to re-code a full library. However, the ease of code reuse in the Julia ecosystem allowed us to solve our problem using existing libraries.