ThorinDistributions

MultivariateGamma(α,θ)

Construct a distribution that correspond to several comonotonous gammas with common shapes α and respectives scales θ[i]. The distribution can then be used through several methods, following the Distributions.jl standard, mainly to obtain samples.

The code is type stable and handles any <:Real types, given by the parameters.

Note that the supprot of this distribution is a half-line in the d-dimensional space, [0,+∞], with angle given by θ.

Examples

julia> dist = MultivariateGamma(2,[3,4,5]);
julia> sample = zeros(Float64,(3,10));
julia> Random.rand!(dist,sample);

MultivariateGammaConvolution(α,θ)

Constructs a distribution that corresponds to the convolutions of MutlivariateGamma(α[i],θ[i,:]) distributions. The distribution can then be used through several methods, following the Distributions.jl standard

The pdf and cdf and not yet coded. The code is type stable and handles any <:Real types, given by the parameters.

Fitting th edistribution is not yet possible.

Examples

julia> dist = MultivariateGammaConvolution([2,3,1],[3 0;4 0;0 1]);
julia> sample = zeros(Float64,(2,10));
julia> Random.rand!(dist,sample);

UnivariateGammaConvolution(α,θ)

Constructs a distribution that corresponds to the convolutions of Gamma(α[i],θ[i]) distributions. The distribution can then be used through several methods, following the Distributions.jl standard, to obtain pdf, cdf, random samples...

The pdf and cdf are handled by the Moshopoulos algorithm, and random samples by simply adding random gammas. The code is type stable and handles any <:Real types, given by the parameters.

To fit the distribution, a loglikelyhood approach could be used. A more involved approach from Furman might be coded sometimes (but requires tanh-sinh integration and bigfloats...)

Examples

julia> dist = UnivariateGammaConvolution([1,0.5, 3.7],[4,2, 10]);
julia> sample = zeros(Float64,10);
julia> Random.rand!(dist,sample);
julia> pdf.((dist,),sample);

Distributions.product_distribution(d1,d2)

Implements the special case of the product of two gamma or univariate gamma convolutions distributions, and output as the result a multivariate gamma convolution with independent margins.

E_from_g(g)

Simply computes the empirical laguerre coefficients from the exponentially shifted moments, if you computed those moments from a theoretical distribution.

L2Objective(par,emp_coefs)

A L2 distance to be minimized between laguerre coefficients of MultivariateGammaConvolutions and empirical laguerre coefficients. This distance is some kind of very high degree polynomial * exponentials, so minimizing it is very hard. We recomand using some global minimisation routine, we found the Particle Swarm optimizers from Optim.jl particularly efficient.

L2ObjectiveWithPenalty(par,emp_coefs)

A L2 distance to be minimized between laguerre coefficients of MultivariateGammaConvolutions and empirical laguerre coefficients. This distance is some kind of very high degree polynomial * exponentials, so minimizing it is very hard. This version includes a penalty to force parameters to go towards 0. but it is not yet working correctly.

MFK_Projection(g_integrals,n_gammas)

From a set of gintegrals computed by the computeg function, this function runs the algorithm from Miles, Furman & Kuznetsov to produce a univariate gamma convolution.

This algorithm works only if the distribution you computed the g_integrals from is, indeed, a generalized gamma convolution. otherwise, it might fail and return garbage.

compute_g(dist,n, integrator)

This function will compute the n first theretical (-1)-exponentially shifted moments of a distribution, using takashi-mori tanh-sinh exponentials formulas from the DoubleExponentialFormulas package.

This might take a long time. As an integrator, you should pass the result of QuadDE(BigFloat) after setting enough precision.

empirical_coefs(x,maxp)

Efficiently Compute the empirical laguerre coefficients of the density of the random vector x, given as a Matrix with n row (number of samples) and d columns (number of dimensions)

get_coefficients(α,θ,m)

α should be a vector of shapes of length n, θ should be a matrix of θ of size (n,d) for a MultivariateGammaConvolution with d marginals.

This function produce the coefficients of the multivariate (tensorised) expensions of the density. It is type stable and quite optimized. it is based on some precomputations that are done in a certain precision.

get_precomp(type,m)

Query and eventualy compute and store the precomputations for laguerre things.

laguerre_L_2x(x,p)

Compute univariate laguerre polynomials Lₚ(2x).

laguerre_density(x,coefs)

Given some laguerre coefficients, computes the correpsonding function at the point x.

laguerre_phi(x,p)

Computes the function ϕₚ(x) = √2 e⁻ˣ Lₚ(2x). This functions form an orthonormal basis of R₊

laguerre_phi_several_pts(x,max_p)

Computes the laguerrephi for each row of x and each p in CartesianIndex(maxp) This is a lot more efficient than broadcasting the laguerrephi function, but this is ugly in the sense that we re-wrote some of the code. This machanisme could be re-written more efficiently.

minimum_m(n,d)

Computes the minimum integer m such that m^d > (d+1)n