ThorinDistributions

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.

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.

resempsthorinmoments(M,D,t_star,m)

M is an integer D is the data, as a (d,n)-shaped matrix, where d is the numebr of dimensions and n the number of samples. t_star is usually taking to be -1 m is an integer.

This resamples the bare thorin moments M times, more efficienty than a simple bootstrap of the thorin_moments(D,t_star,m) function.

thorin_moments(D,t_star,m)

D should be a vector of samples from a distribution, t_star is usually taking to be -1, and m is an integer.

This compute the bare thorin moments, defined as a rescaling of the sample biaised t_star-shifted cumulants.