Implement gradlogpartition for Exponential Family Distributions

[Nimrais]

This PR aims to close issue #130: it implements getgradlogpartition for exponential family distribution types from Distributions.jl

Derivations to do this feature we already prepared.

To test new functionality a new type of test for the exponential family interface was added: test_gradlogpartion_against_expectation. Based on the fact that expectation of the sufficient statistics is equal to the gradient of the logpartition:

$$\mathbf{E}[T(x)] = \nabla_{\eta} A(\eta)$$

However the straightforward check would be to long to run:

$$\mathbf{E}[T(x)]$$

can be generically approximated with a monte-carlo estimator

$$\frac{1}{N}\sum T(z_{i})$$

but this estimator is too slow to converge for the testing purposes.

So another with a faster convergence rate (the linear convergence rate)

$$\mathbf{E}[T(x)]^{T} F(\eta)^{-1} \mathbf{E}[T(x)]= \nabla_{\eta} A(\eta)^{T} F(\eta)^{-1} \nabla_{\eta} A(\eta)$$

is used and a sanity check that dimensionality of the gradient and the natural parameters are the same.

Tasks to do:

• Bernoulli @HoangMHNguyen
• Beta (Add gradient of beta #154) @bartvanerp (Thanks @ismailsenoz)
• Univariate Normal (Add gradient calculation for log partition #157) @bartvanerp
• Multivariate Normal @bartvanerp @Nimrais
• Wishart @Nimrais @ismailsenoz
• InverseWishart @Nimrais @ismailsenoz
• MvNormalWishart @ismailsenoz
• Contingency
• Dirichlet @wouterwln
• MatrixDirichlet @ismailsenoz
• Inverse Gamma @ismailsenoz
• Gamma @wouterwln
• Pareto @Sepideh-Adamiat
• Poisson @HoangMHNguyen
• Laplace @abpolym @ismailsenoz
• Exponential @HoangMHNguyen
• Log Normal @Sepideh-Adamiat @ismailsenoz
• Rayleigh @ismailsenoz
• Negative binomial @ismailsenoz
• Weibull @abpolym @ismailsenoz
• Von Mises @Nimrais
• Von Mises-Fisher @Nimrais @ismailsenoz
• Continuous Bernoulli @Nimrais
• Chi Squared @Nimrais
• Normal-Gamma @abpolym @ismailsenoz
• Geometric @wouterwln
• Erlang
• Multinomial
• Binomial @Sepideh-Adamiat @ismailsenoz
• Categorical @ismailsenoz

So please open a PR with getgradlogpartition for your type implemented into this branch.

[Nimrais]

@bvdmitri, I believe this PR is now ready for review.

Please note that the Contingency, Multinomial, and Continuous Bernoulli distributions will not be implemented in this PR. This is because we currently lack comprehensive functionality for them, and they haven’t been added to the package yet.

I’ve suppressed the test for MvNormalWishart. This is because, first, we need to implement its variate type, I believe. Subsequently, we should refactor this PR with its gradlogpartition: #173 accordingly.

The only issue I foresee is with the MatrixDirichlet. Sampling for it seems to be numerically problematic. Perhaps we could choose a non-random point for evaluation as a workaround?

[Nimrais]

@bvdmitri @ismailsenoz The problem is completely numerical as I see for MatrixDirichlet, it manifests it with even Beta sometimes samples are 0 exactly and log sufficient statistics function is producing NaNs. I tested scipy for this scenario as well and it also can produce samples that are exactly 0.

[bvdmitri]

@Nimrais can you provide a MWE? My very simple thus naive test does not produce 0 samples

julia> any(iszero, rand(Beta(2, 7), 1000000))
false

julia> any(sample -> sample ≈ 0, rand(Beta(2, 7), 1000000))
false

[Nimrais]

@Nimrais can you provide a MWE? My very simple thus naive test does not produce 0 samples

julia> any(iszero, rand(Beta(2, 7), 1000000))
false

julia> any(sample -> sample ≈ 0, rand(Beta(2, 7), 1000000))
false

julia> rand(Beta(0.001, 10), 1000)
1000-element Vector{Float64}:
 2.501697090499602e-61
 1.5363386828054128e-117
 0.0
 2.2934573763809164e-19
 0.0
 5.025249983628666e-188
 3.5214471095207316e-46
 ⋮
 1.2648082681768658e-96
 0.0
 8.892674549717065e-308
 4.820971913271643e-100
 0.0
 5.115976241095818e-165
 0.0

[Nimrais]

You can check our Direchlet test distributions, some marginal of them are pretty close to ill-behaved Betas, one parameter is super small and another is big.

[bvdmitri]

Got it! Lets address this in a separate issue

[codecov]

Codecov Report

Attention: 47 lines in your changes are missing coverage. Please review.

Comparison is base (65b1623) 79.77% compared to head (6f31e89) 80.21%.
Report is 23 commits behind head on main.

Files	Patch %	Lines
src/distributions/wishart.jl	46.15%	7 Missing ⚠️
src/distributions/lognormal.jl	50.00%	5 Missing ⚠️
src/distributions/matrix_dirichlet.jl	50.00%	5 Missing ⚠️
src/distributions/normal_family/normal_family.jl	76.47%	4 Missing ⚠️
src/distributions/weibull.jl	42.85%	4 Missing ⚠️
src/distributions/dirichlet.jl	50.00%	3 Missing ⚠️
src/distributions/gamma_family/gamma_family.jl	50.00%	3 Missing ⚠️
src/distributions/geometric.jl	50.00%	3 Missing ⚠️
src/distributions/negative_binomial.jl	50.00%	3 Missing ⚠️
src/distributions/pareto.jl	50.00%	3 Missing ⚠️
... and 5 more

Additional details and impacted files

@@            Coverage Diff             @@
##             main     #149      +/-   ##
==========================================
+ Coverage   79.77%   80.21%   +0.44%     
==========================================
  Files          39       39              
  Lines        2887     3094     +207     
==========================================
+ Hits         2303     2482     +179     
- Misses        584      612      +28     

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

[bvdmitri]

I trust the tests, didn't check the actual math. Well done guys!!!