Hey,
The loss function that i am optimizing was previously posted there, but is now incorporated in this package, with a type-agnostic interface.
This simple example should be enough to see the problem:
import Random, ThorinDistributions
using DoubleFloats
using ArbNumerics
setprecision(256)
setprecision(ArbFloat,256)
n = 8
d = 2
m = (40,40);
Random.seed!(123)
alpha = rand(n)
scales = reshape(rand(n*d),(n,d))
Errors = Dict(Float64 => (0.0,1.0,0.0), Double64 => (0.0,1.0,0.0), BigFloat => (0.0,1.0,0.0), ArbFloat => (0.0, 1.0, 0.0))
for T in (Float64, Double64, BigFloat, ArbFloat)
alpha = T.(alpha);
scales = T.(scales);
d = ThorinDistributions.MultivariateGammaConvolution(alpha,scales);
samples = zeros(T,(2,10000));
Random.rand!(d,samples);
E = ThorinDistributions.empirical_coefs(samples,m); # empirical coefs
coefs = ThorinDistributions.get_coefficients(alpha, scales,m); # theoreticla coefs
time = Base.@elapsed(ThorinDistributions.get_coefficients(alpha,scales,m))
aloc = @allocated ThorinDistributions.get_coefficients(alpha,scales,m)
# All these should have the same magnitude, very small.
Errors[T] = (Float64(sum((E - coefs)^2)),time,aloc)
end
# Type => (Error, time, allocations)
display(Errors)
julia> # Type => (Error, time, allocations)
julia> display(Errors)
Dict{Type,Tuple{Float64,Float64,Float64}} with 4 entries:
BigFloat => (6.09923e-5, 1.01388, 2.9704e8)
Float64 => (1.01502e20, 0.320554, 1.37131e8)
ArbFloat => (0.000182493, 0.893221, 2.22239e8)
Double64 => (0.000136315, 0.397989, 2.7385e8)
julia>
E and coefs should be very close. Overflow occurs in both the empirical and theoretical coefs, and incrising m, n or d just makes it worse. Note that i did not manage to save the time returned by @btime, so timing might be wrong, but i still ran the thing several time to be shure everything was compiled properly (the package does some pre-computations on ech new type it encounters).
It is weird that arbfloat get a very poor result. I was expecting something more like bigfloat since the precision should be the same.