Skip to content

Commit

Permalink
Tweak default max_proposals limit (#16)
Browse files Browse the repository at this point in the history
* tweak default max proposal limit and improve max proposal error msg
* fix docs move infos and warnings for `GibbsPolarSlice` to docstring
* bump setup-julia action version
  • Loading branch information
Red-Portal authored Nov 30, 2024
1 parent 48059bc commit 7b261ae
Show file tree
Hide file tree
Showing 9 changed files with 42 additions and 31 deletions.
4 changes: 2 additions & 2 deletions .github/workflows/CI.yml
Original file line number Diff line number Diff line change
Expand Up @@ -20,7 +20,7 @@ jobs:
matrix:
version:
- '1'
- '1.10'
- 'lts'
os:
- ubuntu-latest
- macOS-latest
Expand All @@ -29,7 +29,7 @@ jobs:
- x64
steps:
- uses: actions/checkout@v4
- uses: julia-actions/setup-julia@v1
- uses: julia-actions/setup-julia@v2
with:
version: ${{ matrix.version }}
arch: ${{ matrix.arch }}
Expand Down
14 changes: 3 additions & 11 deletions docs/src/gibbs_polar.md
Original file line number Diff line number Diff line change
Expand Up @@ -8,7 +8,6 @@ Unlike other slice sampling algorithms, it operates a Gibbs sampler over polar c
Due to the involvement of polar coordinates, GPSS only works reliably on more than one dimension.
However, unlike ESS, GPSS is applicable to any target distribution.


## Description
For a $$d$$-dimensional target distribution $$\pi$$, GPSS utilizes the following augmented target distribution:
```math
Expand All @@ -34,26 +33,19 @@ The Gibbs steps on $$\theta$$ and $$r$$ are implemented through specialized shri

The only tunable parameter of the algorithm is the size of the search interval (window) of the shrinkage sampler for the radius variable $$r$$.

!!! warning
A limitation of the current implementation of GPSS is that the acceptance rate exhibits a heavy tail. That is, occasionally, a single transition might take an excessive amount of time.

!!! info
The kernel corresponding to this sampler is defined on an **augmented state space** and cannot directly perform a transition on $$x$$.
This also means that the corresponding kernel is not reversible with respect to $$x$$.

## Interface

!!! info
By the nature of polar coordinates, GPSS only works reliably for targets with dimension at least $$d \geq 2$$.

```@docs
GibbsPolarSlice
```

!!! warning
When initializing the chain (*e.g.* the `initial_params` keyword arguments in `AbstractMCMC.sample`), it is necessary to inialize from a point $$x_0$$ that has a sensible norm $$\lVert x_0 \rVert > 0$$, otherwise, the chain will start from a pathologic point in polar coordinates. This might even result in the sampler getting stuck in an infinite loop. (This can be prevented by setting `max_proposals`.) If $$\lVert x_0 \rVert \leq 10^{-5}$$, the current implementation will display a warning.

!!! info
For Turing users: `Turing` might change `initial_params` to match the support of the posterior. This might lead to $$\lVert x_0 \rVert$$ being small, even though the vector you passed to`initial_params` has a sufficiently large norm. If this is suspected, simply try a different initialization value.


## Demonstration
As illustrated in the original paper, GPSS shows good performance on heavy-tailed targets despite being a multivariate slice sampler.
Consider a 10-dimensional Student-$$t$$ target with 1-degree of freedom (this corresponds to a multivariate Cauchy):
Expand Down
19 changes: 13 additions & 6 deletions src/SliceSampling.jl
Original file line number Diff line number Diff line change
Expand Up @@ -20,6 +20,8 @@ export sample, MCMCThreads, MCMCDistributed, MCMCSerial
# Interfaces
abstract type AbstractSliceSampling <: AbstractMCMC.AbstractSampler end

const DEFAULT_MAX_PROPOSALS = 10_000

"""
struct Transition
Expand Down Expand Up @@ -52,10 +54,9 @@ Return the initial sample for the `model` using the random number generator `rng
"""
function initial_sample(::Random.AbstractRNG, ::Any)
error(
"`initial_sample` is not implemented but an initialization wasn't provided. " *
"`initial_sample` is not implemented but an initialization wasn't provided. ",
"Consider supplying an initialization to `initial_params`."
)
println("fuck!!!")
end

# If target is from `LogDensityProblemsAD`, unwrap target before calling `initial_sample`.
Expand All @@ -66,10 +67,16 @@ initial_sample(
) = initial_sample(rng, parent(wrap))

function exceeded_max_prop(max_prop::Int)
error("Exceeded maximum number of proposal $(max_prop).\n",
"Here are possible causes:\n",
"- The model might be broken or pathologic.\n",
"- There might be a bug in the sampler.")
error("Exceeded maximum number of proposal $(max_prop), ",
"which indicates an acceptance rate less than $(1/max_prop*100)%. ",
"A quick fix is to increase `max_prop`, ",
"but an acceptance rate that is too low often indicates that there is a problem. ",
"Here are some possible causes:\n",
"- The model might be broken or degenerate (most likely cause).\n",
"- The tunable parameters of the sampler are suboptimal.\n",
"- The initialization is pathologic. (try supplying a (different) `initial_params`)\n",
"- There might be a bug in the sampler. (if this is suspected, file an issue to `SliceSampling`)\n"
)
end

## Univariate Slice Sampling Algorithms
Expand Down
16 changes: 14 additions & 2 deletions src/multivariate/gibbspolar.jl
Original file line number Diff line number Diff line change
Expand Up @@ -9,14 +9,26 @@ Gibbsian polar slice sampling algorithm by P. Schär, M. Habeck, and D. Rudolf [
# Keyword Arguments
- `w::Real`: Initial window size for the radius shrinkage procedure
- `max_proposals::Int`: Maximum number of proposals allowed until throwing an error (default: `typemax(Int)`).
- `max_proposals::Int`: Maximum number of proposals allowed until throwing an error (default: `$(DEFAULT_MAX_PROPOSALS)`).
!!! info
By the nature of polar coordinates, GPSS only works reliably for targets with dimension at least \$\$d \\geq 2\$\$.
!!! info
The initial window size `w` must be set at least an order of magnitude larger than what is sensible for other slice samplers. Otherwise, a large number of rejections might be experienced.
!!! warning
When initializing the chain (*e.g.* the `initial_params` keyword arguments in `AbstractMCMC.sample`), it is necessary to inialize from a point \$\$x_0\$\$ that has a sensible norm \$\$\\lVert x_0 \\rVert > 0\$\$, otherwise, the chain will start from a pathologic point in polar coordinates. This might even result in the sampler getting stuck in an infinite loop. (This can be prevented by setting `max_proposals`.) If \$\$\\lVert x_0 \\rVert \\leq 10^{-5}\$\$, the current implementation will display a warning.
!!! info
For Turing users: `Turing` might change `initial_params` to match the support of the posterior. This might lead to \$\$\\lVert x_0 \\rVert\$\$ being small, even though the vector you passed to`initial_params` has a sufficiently large norm. If this is suspected, simply try a different initialization value.
"""
struct GibbsPolarSlice{W <: Real} <: AbstractMultivariateSliceSampling
w::W
max_proposals::Int
end

GibbsPolarSlice(w::Real; max_proposals::Int = typemax(Int)) = GibbsPolarSlice(w, max_proposals)
GibbsPolarSlice(w::Real; max_proposals::Int = DEFAULT_MAX_PROPOSALS) = GibbsPolarSlice(w, max_proposals)

struct GibbsPolarSliceState{T <: Transition, R <: Real, D <: AbstractVector}
"Current [`Transition`](@ref)."
Expand Down
4 changes: 2 additions & 2 deletions src/multivariate/latent.jl
Original file line number Diff line number Diff line change
Expand Up @@ -8,14 +8,14 @@ Latent slice sampling algorithm by Li and Walker[^LW2023].
- `beta::Real`: Beta parameter of the Gamma distribution of the auxiliary variables.
# Keyword Arguments
- `max_proposals::Int`: Maximum number of proposals allowed until throwing an error (default: `typemax(Int)`).
- `max_proposals::Int`: Maximum number of proposals allowed until throwing an error (default: `$(DEFAULT_MAX_PROPOSALS)`).
"""
struct LatentSlice{B <: Real} <: AbstractMultivariateSliceSampling
beta ::B
max_proposals::Int
end

function LatentSlice(beta::Real; max_proposals::Int = typemax(Int))
function LatentSlice(beta::Real; max_proposals::Int = DEFAULT_MAX_PROPOSALS)
@assert beta > 0 "Beta must be strictly positive"
LatentSlice(beta, max_proposals)
end
Expand Down
4 changes: 2 additions & 2 deletions src/univariate/doublingout.jl
Original file line number Diff line number Diff line change
Expand Up @@ -9,7 +9,7 @@ Univariate slice sampling by automatically adapting the initial interval through
# Keyword Arguments
- `max_doubling_out`: Maximum number of "doubling outs" (default: 8).
- `max_proposals::Int`: Maximum number of proposals allowed until throwing an error (default: `typemax(Int)`).
- `max_proposals::Int`: Maximum number of proposals allowed until throwing an error (default: `$(DEFAULT_MAX_PROPOSALS)`).
"""
struct SliceDoublingOut{W <: Real} <: AbstractUnivariateSliceSampling
window ::W
Expand All @@ -20,7 +20,7 @@ end
function SliceDoublingOut(
window ::Real;
max_doubling_out::Int = 8,
max_proposals ::Int = typemax(Int),
max_proposals ::Int = DEFAULT_MAX_PROPOSALS,
)
@assert window > 0
SliceDoublingOut(window, max_doubling_out, max_proposals)
Expand Down
4 changes: 2 additions & 2 deletions src/univariate/fixedinterval.jl
Original file line number Diff line number Diff line change
Expand Up @@ -8,7 +8,7 @@ Univariate slice sampling with a fixed initial interval (Scheme 2 by Neal[^N2003
- `window::Real`: Proposal window.
# Keyword Arguments
- `max_proposals::Int`: Maximum number of proposals allowed until throwing an error (default: `typemax(Int)`).
- `max_proposals::Int`: Maximum number of proposals allowed until throwing an error (default: `$(DEFAULT_MAX_PROPOSALS)`).
"""
struct Slice{W <: Real} <: AbstractUnivariateSliceSampling
window ::W
Expand All @@ -17,7 +17,7 @@ end

function Slice(
window ::Real;
max_proposals::Int = typemax(Int),
max_proposals::Int = DEFAULT_MAX_PROPOSALS,
)
@assert window > 0
Slice(window, max_proposals)
Expand Down
4 changes: 2 additions & 2 deletions src/univariate/steppingout.jl
Original file line number Diff line number Diff line change
Expand Up @@ -9,7 +9,7 @@ Univariate slice sampling by automatically adapting the initial interval through
# Keyword Arguments
- `max_stepping_out::Int`: Maximum number of "stepping outs" (default: 32).
- `max_proposals::Int`: Maximum number of proposals allowed until throwing an error (default: `typemax(Int)`).
- `max_proposals::Int`: Maximum number of proposals allowed until throwing an error (default: `$(DEFAULT_MAX_PROPOSALS)`).
"""
struct SliceSteppingOut{W <: Real} <: AbstractUnivariateSliceSampling
window ::W
Expand All @@ -20,7 +20,7 @@ end
function SliceSteppingOut(
window ::Real;
max_stepping_out::Int = 32,
max_proposals ::Int = typemax(Int),
max_proposals ::Int = DEFAULT_MAX_PROPOSALS,
)
@assert window > 0
SliceSteppingOut(window, max_stepping_out, max_proposals)
Expand Down
4 changes: 2 additions & 2 deletions test/multivariate.jl
Original file line number Diff line number Diff line change
Expand Up @@ -78,7 +78,7 @@ end
LatentSlice(5),

# Gibbsian polar slice sampling
GibbsPolarSlice(10),
GibbsPolarSlice(100),
]
@testset "initial_params" begin
model = MultiModel(1.0, 1.0, [0.0])
Expand All @@ -87,7 +87,7 @@ end

θ0 = [1.0, 0.1]
chain = sample(
model,
model,
sampler,
10;
initial_params=θ0,
Expand Down

0 comments on commit 7b261ae

Please sign in to comment.