Skip to content

Inverse gamma - #53

Draft
simonsteiger wants to merge 3 commits into
JuliaBayes:mainfrom
simonsteiger:inverse-gamma
Draft

Inverse gamma#53
simonsteiger wants to merge 3 commits into
JuliaBayes:mainfrom
simonsteiger:inverse-gamma

Conversation

@simonsteiger

Copy link
Copy Markdown
Member

Closes #29.

Same situation as for #52. I'd have liked to point this onto the wishart branch but can't do so because I was working from my fork. Merge #51 and #52 first, then this one.

`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
`InverseGamma(α, θ)` is the measure of `1/X` for `X` distributed as
`Gamma(α, 1/θ)`. Its mean is infinite at or below a unit shape and its
variance at or below a shape of two, both of which it reports rather
than returning a negative number.

Its distribution functions read the upper incomplete gamma integral at
`θ/x`, which puts its lower tail on the integral's upper one. Inverting
that through the lower tail would mean forming `1 - p`, and a
probability of `1e-300` has no complement in floating point. So
`gammaquantile` now takes which tail its probability measures, and
`quantile` inverts the upper one directly, keeping small probabilities
down to the point where the quantile itself underflows.

`valuetype(d, x)` moves to `core/interface.jl`, where both gamma
measures reach it, and `masstype` becomes its discrete case rather than
a second copy of the same promotion.

Assisted-by: Claude Opus 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01HiwLngD494i4dSLJ2DGPBG
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 Inverse gamma

1 participant