Skip to content

Wishart - #52

Draft
simonsteiger wants to merge 2 commits into
JuliaBayes:mainfrom
simonsteiger:wishart
Draft

Wishart#52
simonsteiger wants to merge 2 commits into
JuliaBayes:mainfrom
simonsteiger:wishart

Conversation

@simonsteiger

Copy link
Copy Markdown
Member

Adds Wishart and closes #28.

This work was done on top of #51.

`Gamma(α, θ)` takes a shape and a scale, so its mean is `α * θ`, and
`Gamma(α)` sets the scale to one. The density is closed form and traces;
the distribution functions are not.

`loggammap` and `loggammaq` give the regularized incomplete gamma
integrals in log space. Returning logarithms is what keeps `logcdf` and
`logccdf` finite where the probabilities underflow, and costs nothing:
each tail is already built from a logarithmic prefactor. They sum a
series below `x = a + 1` and run a continued fraction above it, both in
the type they are given, so `BigFloat` keeps its precision. `quantile`
then inverts the tail with Newton's method on `log(x)`, which is what
lets the deep lower tail return a subnormal rather than zero. All of
them loop until their terms stop changing the result, so none can run in
traced or device-side code. Relative accuracy sits at the rounding error
of the argument type for shapes up to about 1000 and falls off roughly
in proportion to the shape after that, since the prefactor's terms grow
while their sum does not.

Sampling uses Marsaglia and Tsang's rejection method, boosted by
`Gamma(α, θ) = Gamma(α + 1, θ) · U^(1/α)` below a unit shape. The accept
step reads `basevalue(α)`, a new interface function returning the plain
floating-point value inside a wrapped number, with methods in the
ForwardDiff and ReverseDiff extensions and an Enzyme inactivity rule.
The loop therefore runs on plain numbers whatever type the parameters
carry, and the accepted noise enters the draw through arithmetic on `α`
and `θ`, which leaves the draw differentiable with respect to both.

The README lists of implemented measures also gain `Poisson`, which they
had been missing.

Assisted-by: Claude Opus 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01HiwLngD494i4dSLJ2DGPBG
`Wishart(ν, L)` is the first matrix-variate measure, so it brings the
`Matrixvariate` variate form, the `ContinuousMatrixvariateMeasure` alias
it dispatches on, and the `PositiveDefiniteMatrices` support.

`L` is the lower-triangular factor of the scale matrix, the convention
`MvNormal` already uses. It keeps `log|S|` a sum over a diagonal and the
trace term a triangular solve rather than an inversion, and it makes
Bartlett's decomposition, `L A A' L'`, plain arithmetic in the
parameters. Sampling needs one chi-squared draw per dimension and so
inherits `Gamma`'s sampler along with its derivative. Draws are
symmetrized, which is what puts them exactly in the support rather than
a rounding error away from it.

`mean`, `var` and `std` take the shape of a draw. `cov` covers every
pair of entries and so is indexed the way `vec` orders them; the
conformance suite gains a `matrixsummaries` group for that shape.

`src/core/linalg.jl` collects the triangular linear algebra the two
factored measures share: `rowdot`, moved from `MvNormal`, alongside
`rowsdot`, `forwardsolve`, `logdetdiag` and `cholfactor`. Each builds
new arrays instead of writing into one, so reverse-mode backends, which
reject array mutation, can follow them. `cholfactor` takes its pivots
through the new `sqrtt`, so an indefinite argument yields a non-finite
factor and the log-density reports `NaN` instead of throwing.

Assisted-by: Claude Opus 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01HiwLngD494i4dSLJ2DGPBG
@simonsteiger

simonsteiger commented Sep 3, 2026

Copy link
Copy Markdown
Member Author

Hmm, I don't like that these are stacked. I'll separate them first!

I regret working from my fork instead of this repo directly. I'd like to rebase this PR onto gamma but can't because gamma is a branch on my fork. 🤦🏻‍♂️ Should merge #51 first for a cleaner diff.

@simonsteiger simonsteiger reopened this Sep 3, 2026
@simonsteiger
simonsteiger marked this pull request as draft September 3, 2026 05:33
@simonsteiger simonsteiger mentioned this pull request Sep 3, 2026
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

Add Wishart

1 participant