diff --git a/dev/.documenter-siteinfo.json b/dev/.documenter-siteinfo.json index 4076b84..f46fc9d 100644 --- a/dev/.documenter-siteinfo.json +++ b/dev/.documenter-siteinfo.json @@ -1 +1 @@ -{"documenter":{"julia_version":"1.10.4","generation_timestamp":"2024-08-12T02:29:09","documenter_version":"1.5.0"}} \ No newline at end of file +{"documenter":{"julia_version":"1.10.4","generation_timestamp":"2024-08-12T20:54:51","documenter_version":"1.5.0"}} \ No newline at end of file diff --git a/dev/abstractmcmc_demo.svg b/dev/abstractmcmc_demo.svg index 57e9245..40649ce 100644 --- a/dev/abstractmcmc_demo.svg +++ b/dev/abstractmcmc_demo.svg @@ -1,48 +1,46 @@ - + - + - + - + - + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/dev/general/index.html b/dev/general/index.html index cf21095..89a341e 100644 --- a/dev/general/index.html +++ b/dev/general/index.html @@ -43,8 +43,8 @@ Iterations = 1:1:10000 Number of chains = 1 Samples per chain = 10000 -Wall duration = 4.19 seconds -Compute duration = 4.19 seconds +Wall duration = 4.34 seconds +Compute duration = 4.34 seconds parameters = s, m internals = lp @@ -52,14 +52,63 @@ parameters mean std mcse ess_bulk ess_tail rhat ⋯ Symbol Float64 Float64 Float64 Float64 Float64 Float64 ⋯ - s 1.5240 3.5438 0.0400 5837.6548 4743.5614 1.0002 ⋯ - m -0.0051 1.2245 0.0129 8940.8834 5833.9653 0.9999 ⋯ + s 1.5399 1.4716 0.0245 5546.7713 4499.4552 1.0001 ⋯ + m -0.0006 1.2366 0.0135 8507.9126 5682.3507 1.0006 ⋯ 1 column omitted Quantiles parameters 2.5% 25.0% 50.0% 75.0% 97.5% Symbol Float64 Float64 Float64 Float64 Float64 - s 0.4194 0.7580 1.1085 1.7324 4.8389 - m -2.5320 -0.7309 -0.0006 0.7249 2.4190 -

Drawing Samples

For drawing samples using the algorithms provided by SliceSampling, the user only needs to call:

sample([rng,] model, slice, N; initial_params)
  • slice::AbstractSliceSampling: Any slice sampling algorithm provided by SliceSampling.
  • model: A model implementing the LogDensityProblems interface.
  • N: The number of samples

The output is a SliceSampling.Transition object, which contains the following:

SliceSampling.Transition — Type
struct Transition

Struct containing the results of the transition.

Fields

  • params: Samples generated by the transition.
  • lp::Real: Log-target density of the samples.
  • info::NamedTuple: Named tuple containing information about the transition.
source

For the keyword arguments, SliceSampling allows:

  • initial_params: The intial state of the Markov chain (default: nothing).

If initial_params is nothing, the following function can be implemented to provide an initialization:

SliceSampling.initial_sample — Function
initial_sample(rng, model)

Return the initial sample for the model using the random number generator rng.

Arguments

  • rng::Random.AbstractRNG: Random number generator.
  • model: The target LogDensityProblem.
source

Performing a Single Transition

For more fined-grained control, the user can call AbstractMCMC.step. That is, the chain can be initialized by calling:

transition, state = AbstractMCMC.steps([rng,] model, slice; initial_params)

and then each MCMC transition on state can be performed by calling:

transition, state = AbstractMCMC.steps([rng,] model, slice, state)

For more details, refer to the documentation of AbstractMCMC.

+ s 0.4144 0.7681 1.1400 1.7783 5.0932 + m -2.4238 -0.7387 -0.0069 0.7363 2.4593 +

Conditional sampling in a Turing.Experimental.Gibbs sampler

SliceSampling.jl be used as a conditional sampler in Turing.Experimental.Gibbs.

using Distributions
+using Turing
+using SliceSampling
+
+@model function simple_choice(xs)
+    p ~ Beta(2, 2)
+    z ~ Bernoulli(p)
+    for i in 1:length(xs)
+        if z == 1
+            xs[i] ~ Normal(0, 1)
+        else
+            xs[i] ~ Normal(2, 1)
+        end
+    end
+end
+
+sampler = Turing.Experimental.Gibbs(
+    (
+        p = externalsampler(SliceSteppingOut(2.0)),
+        z = PG(20, :z)
+    )
+)
+
+n_samples = 1000
+model     = simple_choice([1.5, 2.0, 0.3])
+sample(model, sampler, n_samples)
Chains MCMC chain (1000×3×1 Array{Float64, 3}):
+
+Iterations        = 1:1:1000
+Number of chains  = 1
+Samples per chain = 1000
+Wall duration     = 20.19 seconds
+Compute duration  = 20.19 seconds
+parameters        = p, z
+internals         = lp
+
+Summary Statistics
+  parameters      mean       std      mcse   ess_bulk   ess_tail      rhat   e ⋯
+      Symbol   Float64   Float64   Float64    Float64    Float64   Float64     ⋯
+
+           p    0.6038    0.1975    0.0063   964.8447   619.8975    1.0005     ⋯
+           z    0.3010    0.4589    0.0166   767.5188        NaN    0.9994     ⋯
+                                                                1 column omitted
+
+Quantiles
+  parameters      2.5%     25.0%     50.0%     75.0%     97.5% 
+      Symbol   Float64   Float64   Float64   Float64   Float64 
+
+           p    0.2117    0.4566    0.6201    0.7511    0.9319
+           z    0.0000    0.0000    0.0000    1.0000    1.0000
+

Drawing Samples

For drawing samples using the algorithms provided by SliceSampling, the user only needs to call:

sample([rng,] model, slice, N; initial_params)
  • slice::AbstractSliceSampling: Any slice sampling algorithm provided by SliceSampling.
  • model: A model implementing the LogDensityProblems interface.
  • N: The number of samples

The output is a SliceSampling.Transition object, which contains the following:

SliceSampling.Transition — Type
struct Transition

Struct containing the results of the transition.

Fields

  • params: Samples generated by the transition.
  • lp::Real: Log-target density of the samples.
  • info::NamedTuple: Named tuple containing information about the transition.
source

For the keyword arguments, SliceSampling allows:

  • initial_params: The intial state of the Markov chain (default: nothing).

If initial_params is nothing, the following function can be implemented to provide an initialization:

SliceSampling.initial_sample — Function
initial_sample(rng, model)

Return the initial sample for the model using the random number generator rng.

Arguments

  • rng::Random.AbstractRNG: Random number generator.
  • model: The target LogDensityProblem.
source

Performing a Single Transition

For more fined-grained control, the user can call AbstractMCMC.step. That is, the chain can be initialized by calling:

transition, state = AbstractMCMC.steps([rng,] model, slice; initial_params)

and then each MCMC transition on state can be performed by calling:

transition, state = AbstractMCMC.steps([rng,] model, slice, state)

For more details, refer to the documentation of AbstractMCMC.

diff --git a/dev/gibbs_polar/index.html b/dev/gibbs_polar/index.html index 72e6e1d..51aff17 100644 --- a/dev/gibbs_polar/index.html +++ b/dev/gibbs_polar/index.html @@ -8,7 +8,7 @@ \theta_n &\sim \operatorname{Uniform}\left\{ \theta \in \mathbb{S}^{d-1} \mid \varrho^{(1)}\left(r_{n-1} \theta\right) > T_n \right\} \\ r_n &\sim \operatorname{Uniform}\left\{ r \in \mathbb{R}_{\geq 0} \mid \varrho^{(1)}\left(r \theta_n\right) > T_n \right\} \\ x_n &= \theta_n r_n, -\end{aligned}\]

where $T_n$ is the usual acceptance threshold auxiliary variable, while $\theta$ and $r$ are the sampler states in polar coordinates. The Gibbs steps on $\theta$ and $r$ are implemented through specialized shrinkage procedures.

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

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$.

SliceSampling.GibbsPolarSlice — Type
GibbsPolarSlice(w; max_proposals)

Gibbsian polar slice sampling algorithm by P. Schär, M. Habeck, and D. Rudolf [SHR2023].

Arguments

  • w::Real: Initial window size for the radius shrinkage procedure.

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)).
source
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 toinitial_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):

using Distributions
+\end{aligned}\]
where $T_n$ is the usual acceptance threshold auxiliary variable, while $\theta$ and $r$ are the sampler states in polar coordinates. The Gibbs steps on $\theta$ and $r$ are implemented through specialized shrinkage procedures.
The only tunable parameter of the algorithm is the size of the search interval (window) of the shrinkage sampler for the radius variable $r$.
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$.
SliceSampling.GibbsPolarSlice — Type
GibbsPolarSlice(w; max_proposals)
Gibbsian polar slice sampling algorithm by P. Schär, M. Habeck, and D. Rudolf [SHR2023].
Arguments
  • w::Real: Initial window size for the radius shrinkage procedure.
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)).
source
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 toinitial_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):
using Distributions
 using Turing
 using SliceSampling
 using LinearAlgebra
@@ -27,4 +27,4 @@
 p1 = Plots.plot(1:n_samples, latent_chain[:,1,:], ylims=[-10,10], label="LSS")
 p2 = Plots.plot(1:n_samples, polar_chain[:,1,:],  ylims=[-10,10], label="GPSS")
 plot(p1, p2, layout = l)
-savefig("student_latent_gpss.svg")
"/home/runner/work/SliceSampling.jl/SliceSampling.jl/docs/build/student_latent_gpss.svg"
Clearly, GPSS is better at exploring the deep tails compared to the latent slice sampler (LSS) despite having a similar per-iteration cost.
  • SHR2023Schär, P., Habeck, M., & Rudolf, D. (2023, July). Gibbsian polar slice sampling. In International Conference on Machine Learning.
  • RR2002Roberts, G. O., & Rosenthal, J. S. (2002). The polar slice sampler. Stochastic Models, 18(2), 257-280.

+savefig("student_latent_gpss.svg")
"/home/runner/work/SliceSampling.jl/SliceSampling.jl/docs/build/student_latent_gpss.svg"

Clearly, GPSS is better at exploring the deep tails compared to the latent slice sampler (LSS) despite having a similar per-iteration cost.

  • SHR2023Schär, P., Habeck, M., & Rudolf, D. (2023, July). Gibbsian polar slice sampling. In International Conference on Machine Learning.
  • RR2002Roberts, G. O., & Rosenthal, J. S. (2002). The polar slice sampler. Stochastic Models, 18(2), 257-280.
diff --git a/dev/index.html b/dev/index.html index a9bff8c..082eee4 100644 --- a/dev/index.html +++ b/dev/index.html @@ -1,2 +1,2 @@ -Home · SliceSampling.jl

SliceSampling

This package implements slice sampling algorithms. Slice sampling finds its roots in the Swendsen–Wang algorithm for Ising models[SW1987][ES1988]. It later came into the interest of the statistical community through Besag and Green[BG1993], and popularized by Neal [N2003]. Furthermore, Neal introduced various ways to efficiently implement slice samplers. This package provides the original slice sampling algorithms by Neal and their later extensions.

  • SW1987Swendsen, R. H., & Wang, J. S. (1987). Nonuniversal critical dynamics in Monte Carlo simulations. Physical review letters, 58(2), 86.
  • ES1988Edwards, R. G., & Sokal, A. D. (1988). Generalization of the fortuin-kasteleyn-swendsen-wang representation and monte carlo algorithm. Physical review D, 38(6), 2009.
  • BG1993Besag, J., & Green, P. J. (1993). Spatial statistics and Bayesian computation. Journal of the Royal Statistical Society Series B: Statistical Methodology, 55(1), 25-37.
  • N2003Neal, R. M. (2003). Slice sampling. The annals of statistics, 31(3), 705-767.
+Home · SliceSampling.jl

SliceSampling

This package implements slice sampling algorithms. Slice sampling finds its roots in the Swendsen–Wang algorithm for Ising models[SW1987][ES1988]. It later came into the interest of the statistical community through Besag and Green[BG1993], and popularized by Neal [N2003]. Furthermore, Neal introduced various ways to efficiently implement slice samplers. This package provides the original slice sampling algorithms by Neal and their later extensions.

  • SW1987Swendsen, R. H., & Wang, J. S. (1987). Nonuniversal critical dynamics in Monte Carlo simulations. Physical review letters, 58(2), 86.
  • ES1988Edwards, R. G., & Sokal, A. D. (1988). Generalization of the fortuin-kasteleyn-swendsen-wang representation and monte carlo algorithm. Physical review D, 38(6), 2009.
  • BG1993Besag, J., & Green, P. J. (1993). Spatial statistics and Bayesian computation. Journal of the Royal Statistical Society Series B: Statistical Methodology, 55(1), 25-37.
  • N2003Neal, R. M. (2003). Slice sampling. The annals of statistics, 31(3), 705-767.
diff --git a/dev/latent_slice/index.html b/dev/latent_slice/index.html index 85da91e..247fff3 100644 --- a/dev/latent_slice/index.html +++ b/dev/latent_slice/index.html @@ -4,4 +4,4 @@ s_n &\sim p(s \mid x_{n-1}, l_{n}) \\ t_n &\sim \operatorname{Uniform}\left(0, \pi\left(x_{n-1}\right)\right) \\ x_n &\sim \operatorname{Uniform}\left\{x \mid \pi\left(x\right) > t_n\right\}, -\end{aligned}\]

When $x_n$ is updated using the usual shrinkage procedure of Neal[N2003], $s_n$ and $l_n$ are used to form the initial search window. ($s_n$ is the width of the window and $l_n$ is its center point.) Therefore, the latent slice sampler can be regarded as an automatic tuning mechanism of the "initial interval" of slice samplers.

The only tunable parameter of the algorithm is then the distribution of the width $p(s)$. For this, Li and Walker recommend

\[ p(s; \beta) = \operatorname{Gamma}(s; 2, \beta),\]

where $\beta$ is a tunable parameter. The use of the gamma distribution is somewhat important since the complete conditional $p(s \mid y, l)$ needs to be available in closed-form for efficiency. (It is a truncated exponential distribution in case of the gamma.) Therefore, we only provide control over $\beta$.

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

SliceSampling.LatentSlice — Type
LatentSlice(beta)

Latent slice sampling algorithm by Li and Walker[LW2023].

Arguments

  • 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)).
source
  • LW2023Li, Y., & Walker, S. G. (2023). A latent slice sampling algorithm. Computational Statistics & Data Analysis, 179, 107652.
  • N2003Neal, R. M. (2003). Slice sampling. The annals of statistics, 31(3), 705-767.
+\end{aligned}\]

When $x_n$ is updated using the usual shrinkage procedure of Neal[N2003], $s_n$ and $l_n$ are used to form the initial search window. ($s_n$ is the width of the window and $l_n$ is its center point.) Therefore, the latent slice sampler can be regarded as an automatic tuning mechanism of the "initial interval" of slice samplers.

The only tunable parameter of the algorithm is then the distribution of the width $p(s)$. For this, Li and Walker recommend

\[ p(s; \beta) = \operatorname{Gamma}(s; 2, \beta),\]

where $\beta$ is a tunable parameter. The use of the gamma distribution is somewhat important since the complete conditional $p(s \mid y, l)$ needs to be available in closed-form for efficiency. (It is a truncated exponential distribution in case of the gamma.) Therefore, we only provide control over $\beta$.

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

SliceSampling.LatentSlice — Type
LatentSlice(beta)

Latent slice sampling algorithm by Li and Walker[LW2023].

Arguments

  • 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)).
source
  • LW2023Li, Y., & Walker, S. G. (2023). A latent slice sampling algorithm. Computational Statistics & Data Analysis, 179, 107652.
  • N2003Neal, R. M. (2003). Slice sampling. The annals of statistics, 31(3), 705-767.
diff --git a/dev/meta_multivariate/index.html b/dev/meta_multivariate/index.html index e2a4c57..ec9bc26 100644 --- a/dev/meta_multivariate/index.html +++ b/dev/meta_multivariate/index.html @@ -1,5 +1,5 @@ -Meta Multivariate Samplers · SliceSampling.jl

Meta Multivariate Samplers

To use univariate slice sampling strategies on targets with more than on dimension, one has to embed them into a "meta" multivariate sampling scheme that relies on univariate sampling elements. The two most popular approaches for this are Gibbs sampling[GG1984] and hit-and-run[BRS1993].

Random Permutation Gibbs

Gibbs sampling[GG1984] is a strategy where we sample from the posterior one coordinate at a time, conditioned on the values of all other coordinates. In practice, one can pick the coordinates in any order they want as long as it does not depend on the state of the chain. It is generally hard to know a-prior which "scan order" is best, but randomly picking coordinates tend to work well in general. Currently, we only provide random permutation scan, which guarantees that all coordinates are updated at least once after $d$ transitions. At the same time, reversibility is maintained by randomly permuting the order we go through each coordinate:

SliceSampling.RandPermGibbs — Type
RandPermGibbs(unislice)

Random permutation coordinate-wise Gibbs sampling strategy. This applies unislice coordinate-wise in a random order.

Arguments

  • unislice::Union{<:AbstractUnivariateSliceSampling,<:AbstractVector{<:AbstractUnivariateSliceSampling}}: a single or a vector of univariate slice sampling algorithms.

When unislice is a vector of samplers, each slice sampler is applied to the corresponding coordinate of the target posterior. In that case, the length(unislice) must match the dimensionality of the posterior.

source

Each call to AbstractMCMC.step internally performs $d$ Gibbs transition so that all coordinates are updated.

For example:

RandPermGibbs(SliceSteppingOut(2.))

If one wants to use a different slice sampler configuration for each coordinate, one can mix-and-match by passing a Vector of slice samplers, one for each coordinate. For instance, for a 2-dimensional target:

RandPermGibbs([SliceSteppingOut(2.; max_proposals=32), SliceDoublingOut(2.),])

Hit-and-Run

Hit-and-run is a simple meta algorithm where we sample over a random 1-dimensional projection of the space. That is, at each iteration, we sample a random direction

\[ \theta_n \sim \operatorname{Uniform}(\mathbb{S}^{d-1}),\]

and perform a Markov transition along the 1-dimensional subspace

\[\begin{aligned} +Meta Multivariate Samplers · SliceSampling.jl

Meta Multivariate Samplers

To use univariate slice sampling strategies on targets with more than on dimension, one has to embed them into a "meta" multivariate sampling scheme that relies on univariate sampling elements. The two most popular approaches for this are Gibbs sampling[GG1984] and hit-and-run[BRS1993].

Random Permutation Gibbs

Gibbs sampling[GG1984] is a strategy where we sample from the posterior one coordinate at a time, conditioned on the values of all other coordinates. In practice, one can pick the coordinates in any order they want as long as it does not depend on the state of the chain. It is generally hard to know a-prior which "scan order" is best, but randomly picking coordinates tend to work well in general. Currently, we only provide random permutation scan, which guarantees that all coordinates are updated at least once after $d$ transitions. At the same time, reversibility is maintained by randomly permuting the order we go through each coordinate:

SliceSampling.RandPermGibbs — Type
RandPermGibbs(unislice)

Random permutation coordinate-wise Gibbs sampling strategy. This applies unislice coordinate-wise in a random order.

Arguments

  • unislice::Union{<:AbstractUnivariateSliceSampling,<:AbstractVector{<:AbstractUnivariateSliceSampling}}: a single or a vector of univariate slice sampling algorithms.

When unislice is a vector of samplers, each slice sampler is applied to the corresponding coordinate of the target posterior. In that case, the length(unislice) must match the dimensionality of the posterior.

source

Each call to AbstractMCMC.step internally performs $d$ Gibbs transition so that all coordinates are updated.

For example:

RandPermGibbs(SliceSteppingOut(2.))

If one wants to use a different slice sampler configuration for each coordinate, one can mix-and-match by passing a Vector of slice samplers, one for each coordinate. For instance, for a 2-dimensional target:

RandPermGibbs([SliceSteppingOut(2.; max_proposals=32), SliceDoublingOut(2.),])

Hit-and-Run

Hit-and-run is a simple meta algorithm where we sample over a random 1-dimensional projection of the space. That is, at each iteration, we sample a random direction

\[ \theta_n \sim \operatorname{Uniform}(\mathbb{S}^{d-1}),\]

and perform a Markov transition along the 1-dimensional subspace

\[\begin{aligned} \lambda_n &\sim p\left(\lambda \mid x_{n-1}, \theta_n \right) \propto \pi\left( x_{n-1} + \lambda \theta_n \right) \\ x_{n} &= x_{n-1} + \lambda_n \theta_n, -\end{aligned}\]

where $\pi$ is the target unnormalized density. Applying slice sampling for the 1-dimensional subproblem has been popularized by David Mackay[M2003], and is, technically, also a Gibbs sampler. (Or is that Gibbs samplers are hit-and-run samplers?) Unlike RandPermGibbs, which only makes axis-aligned moves, HitAndRun can choose arbitrary directions, which could be helpful in some cases.

SliceSampling.HitAndRun — Type
HitAndRun(unislice)

Hit-and-run sampling strategy[BRS1993]. This applies unislice along a random direction uniform sampled from the sphere.

Arguments

  • unislice::AbstractUnivariateSliceSampling: Univariate slice sampling algorithm.
source

This can be used, for example, as follows:

HitAndRun(SliceSteppingOut(2.))

Unlike RandPermGibbs, HitAndRun does not provide the option of using a unique unislice object for each coordinate. This is a natural limitation of the hit-and-run sampler: it does not operate on individual coordinates.

  • GG1984Geman, S., & Geman, D. (1984). Stochastic relaxation, Gibbs distributions, and the Bayesian restoration of images. IEEE Transactions on Pattern Analysis and Machine Intelligence, (6).
  • BRS1993Bélisle, C. J., Romeijn, H. E., & Smith, R. L. (1993). Hit-and-run algorithms for generating multivariate distributions. Mathematics of Operations Research, 18(2), 255-266.
  • M2003MacKay, D. J. (2003). Information theory, inference and learning algorithms. Cambridge university press.
+\end{aligned}\]

where $\pi$ is the target unnormalized density. Applying slice sampling for the 1-dimensional subproblem has been popularized by David Mackay[M2003], and is, technically, also a Gibbs sampler. (Or is that Gibbs samplers are hit-and-run samplers?) Unlike RandPermGibbs, which only makes axis-aligned moves, HitAndRun can choose arbitrary directions, which could be helpful in some cases.

SliceSampling.HitAndRun — Type
HitAndRun(unislice)

Hit-and-run sampling strategy[BRS1993]. This applies unislice along a random direction uniform sampled from the sphere.

Arguments

  • unislice::AbstractUnivariateSliceSampling: Univariate slice sampling algorithm.
source

This can be used, for example, as follows:

HitAndRun(SliceSteppingOut(2.))

Unlike RandPermGibbs, HitAndRun does not provide the option of using a unique unislice object for each coordinate. This is a natural limitation of the hit-and-run sampler: it does not operate on individual coordinates.

  • GG1984Geman, S., & Geman, D. (1984). Stochastic relaxation, Gibbs distributions, and the Bayesian restoration of images. IEEE Transactions on Pattern Analysis and Machine Intelligence, (6).
  • BRS1993Bélisle, C. J., Romeijn, H. E., & Smith, R. L. (1993). Hit-and-run algorithms for generating multivariate distributions. Mathematics of Operations Research, 18(2), 255-266.
  • M2003MacKay, D. J. (2003). Information theory, inference and learning algorithms. Cambridge university press.
diff --git a/dev/objects.inv b/dev/objects.inv index 038f627..af7d8cf 100644 Binary files a/dev/objects.inv and b/dev/objects.inv differ diff --git a/dev/search_index.js b/dev/search_index.js index a1fed3d..6e4e0df 100644 --- a/dev/search_index.js +++ b/dev/search_index.js @@ -1,3 +1,3 @@ var documenterSearchIndex = {"docs": -[{"location":"latent_slice/#latent","page":"Latent Slice Sampling","title":"Latent Slice Sampling","text":"","category":"section"},{"location":"latent_slice/#Introduction","page":"Latent Slice Sampling","title":"Introduction","text":"","category":"section"},{"location":"latent_slice/","page":"Latent Slice Sampling","title":"Latent Slice Sampling","text":"Latent slice sampling is a recent vector-valued slice sampling algorithm proposed by Li and Walker[LW2023]. Unlike other slice sampling algorithms, it treats the \"search intervals\" as auxiliary variables and adapts them along the samples from the log-target in a Gibbs-type scheme.","category":"page"},{"location":"latent_slice/#Description","page":"Latent Slice Sampling","title":"Description","text":"","category":"section"},{"location":"latent_slice/","page":"Latent Slice Sampling","title":"Latent Slice Sampling","text":"Specifically, the extended joint density of the latent slice sampler is as follows:","category":"page"},{"location":"latent_slice/","page":"Latent Slice Sampling","title":"Latent Slice Sampling","text":" p(x t s l) = pi(x) p(s) operatornameUniformleft(t 0 pileft(xright)right) operatornameUniformleft(l x - s2 x + s2right)","category":"page"},{"location":"latent_slice/","page":"Latent Slice Sampling","title":"Latent Slice Sampling","text":"where y is the parameters of the log-target pi, s is the width of the search interval and l is the centering of the search interval relative to y. Naturally, the sampler operates as a blocked-Gibbs sampler ","category":"page"},{"location":"latent_slice/","page":"Latent Slice Sampling","title":"Latent Slice Sampling","text":"beginaligned\nl_n sim operatornameUniformleft(l x_n-1 - s_n-12 x_n-1 + s_n-12right) \ns_n sim p(s mid x_n-1 l_n) \nt_n sim operatornameUniformleft(0 pileft(x_n-1right)right) \nx_n sim operatornameUniformleftx mid pileft(xright) t_nright\nendaligned","category":"page"},{"location":"latent_slice/","page":"Latent Slice Sampling","title":"Latent Slice Sampling","text":"When x_n is updated using the usual shrinkage procedure of Neal[N2003], s_n and l_n are used to form the initial search window. (s_n is the width of the window and l_n is its center point.) Therefore, the latent slice sampler can be regarded as an automatic tuning mechanism of the \"initial interval\" of slice samplers.","category":"page"},{"location":"latent_slice/","page":"Latent Slice Sampling","title":"Latent Slice Sampling","text":"The only tunable parameter of the algorithm is then the distribution of the width p(s). For this, Li and Walker recommend","category":"page"},{"location":"latent_slice/","page":"Latent Slice Sampling","title":"Latent Slice Sampling","text":" p(s beta) = operatornameGamma(s 2 beta)","category":"page"},{"location":"latent_slice/","page":"Latent Slice Sampling","title":"Latent Slice Sampling","text":"where beta is a tunable parameter. The use of the gamma distribution is somewhat important since the complete conditional p(s mid y l) needs to be available in closed-form for efficiency. (It is a truncated exponential distribution in case of the gamma.) Therefore, we only provide control over beta.","category":"page"},{"location":"latent_slice/","page":"Latent Slice Sampling","title":"Latent Slice Sampling","text":"info: Info\nThe 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.","category":"page"},{"location":"latent_slice/#Interface","page":"Latent Slice Sampling","title":"Interface","text":"","category":"section"},{"location":"latent_slice/","page":"Latent Slice Sampling","title":"Latent Slice Sampling","text":"LatentSlice","category":"page"},{"location":"latent_slice/#SliceSampling.LatentSlice","page":"Latent Slice Sampling","title":"SliceSampling.LatentSlice","text":"LatentSlice(beta)\n\nLatent slice sampling algorithm by Li and Walker[LW2023].\n\nArguments\n\nbeta::Real: Beta parameter of the Gamma distribution of the auxiliary variables.\n\nKeyword Arguments\n\nmax_proposals::Int: Maximum number of proposals allowed until throwing an error (default: typemax(Int)).\n\n\n\n\n\n","category":"type"},{"location":"latent_slice/","page":"Latent Slice Sampling","title":"Latent Slice Sampling","text":"[LW2023]: Li, Y., & Walker, S. G. (2023). A latent slice sampling algorithm. Computational Statistics & Data Analysis, 179, 107652.","category":"page"},{"location":"latent_slice/","page":"Latent Slice Sampling","title":"Latent Slice Sampling","text":"[N2003]: Neal, R. M. (2003). Slice sampling. The annals of statistics, 31(3), 705-767.","category":"page"},{"location":"gibbs_polar/#polar","page":"Gibbsian Polar Slice Sampling","title":"Gibbsian Polar Slice Sampling","text":"","category":"section"},{"location":"gibbs_polar/#Introduction","page":"Gibbsian Polar Slice Sampling","title":"Introduction","text":"","category":"section"},{"location":"gibbs_polar/","page":"Gibbsian Polar Slice Sampling","title":"Gibbsian Polar Slice Sampling","text":"Gibbsian polar slice sampling (GPSS) is a recent vector-valued slice sampling algorithm proposed by P. Schär, M. Habeck, and D. Rudolf[SHR2023]. It is an computationally efficient variant of polar slice sampler previously proposed by Roberts and Rosenthal[RR2002]. Unlike other slice sampling algorithms, it operates a Gibbs sampler over polar coordinates, reminiscent of the elliptical slice sampler (ESS). 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.","category":"page"},{"location":"gibbs_polar/#Description","page":"Gibbsian Polar Slice Sampling","title":"Description","text":"","category":"section"},{"location":"gibbs_polar/","page":"Gibbsian Polar Slice Sampling","title":"Gibbsian Polar Slice Sampling","text":"For a d-dimensional target distribution pi, GPSS utilizes the following augmented target distribution:","category":"page"},{"location":"gibbs_polar/","page":"Gibbsian Polar Slice Sampling","title":"Gibbsian Polar Slice Sampling","text":"beginaligned\n p(x T) = varrho_pi^(0)(x) varrho_pi^(1)(x) operatornameUniformleft(T 0 varrho^1(x)right) \n varrho_pi^(0)(x) = lVert x rVert^1 - d \n varrho_pi^(1)(x) = lVert x rVert^d-1 pileft(xright)\nendaligned","category":"page"},{"location":"gibbs_polar/","page":"Gibbsian Polar Slice Sampling","title":"Gibbsian Polar Slice Sampling","text":"As described in Appendix A of the GPSS paper, sampling from varrho^(1)(x) in polar coordinates magically targets the augmented target distribution.","category":"page"},{"location":"gibbs_polar/","page":"Gibbsian Polar Slice Sampling","title":"Gibbsian Polar Slice Sampling","text":"In a high-level view, GPSS operates a Gibbs sampler in the following fashion:","category":"page"},{"location":"gibbs_polar/","page":"Gibbsian Polar Slice Sampling","title":"Gibbsian Polar Slice Sampling","text":"beginaligned\nT_n sim operatornameUniformleft(0 varrho^(1)left(x_n-1right)right) \ntheta_n sim operatornameUniformleft theta in mathbbS^d-1 mid varrho^(1)left(r_n-1 thetaright) T_n right \nr_n sim operatornameUniformleft r in mathbbR_geq 0 mid varrho^(1)left(r theta_nright) T_n right \nx_n = theta_n r_n\nendaligned","category":"page"},{"location":"gibbs_polar/","page":"Gibbsian Polar Slice Sampling","title":"Gibbsian Polar Slice Sampling","text":"where T_n is the usual acceptance threshold auxiliary variable, while theta and r are the sampler states in polar coordinates. The Gibbs steps on theta and r are implemented through specialized shrinkage procedures.","category":"page"},{"location":"gibbs_polar/","page":"Gibbsian Polar Slice Sampling","title":"Gibbsian Polar Slice Sampling","text":"The only tunable parameter of the algorithm is the size of the search interval (window) of the shrinkage sampler for the radius variable r.","category":"page"},{"location":"gibbs_polar/","page":"Gibbsian Polar Slice Sampling","title":"Gibbsian Polar Slice Sampling","text":"info: Info\nThe 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.","category":"page"},{"location":"gibbs_polar/#Interface","page":"Gibbsian Polar Slice Sampling","title":"Interface","text":"","category":"section"},{"location":"gibbs_polar/","page":"Gibbsian Polar Slice Sampling","title":"Gibbsian Polar Slice Sampling","text":"info: Info\nBy the nature of polar coordinates, GPSS only works reliably for targets with dimension at least d geq 2.","category":"page"},{"location":"gibbs_polar/","page":"Gibbsian Polar Slice Sampling","title":"Gibbsian Polar Slice Sampling","text":"GibbsPolarSlice","category":"page"},{"location":"gibbs_polar/#SliceSampling.GibbsPolarSlice","page":"Gibbsian Polar Slice Sampling","title":"SliceSampling.GibbsPolarSlice","text":"GibbsPolarSlice(w; max_proposals)\n\nGibbsian polar slice sampling algorithm by P. Schär, M. Habeck, and D. Rudolf [SHR2023].\n\nArguments\n\nw::Real: Initial window size for the radius shrinkage procedure.\n\nKeyword Arguments\n\nw::Real: Initial window size for the radius shrinkage procedure\nmax_proposals::Int: Maximum number of proposals allowed until throwing an error (default: typemax(Int)).\n\n\n\n\n\n","category":"type"},{"location":"gibbs_polar/","page":"Gibbsian Polar Slice Sampling","title":"Gibbsian Polar Slice Sampling","text":"warning: Warning\nWhen 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. ","category":"page"},{"location":"gibbs_polar/","page":"Gibbsian Polar Slice Sampling","title":"Gibbsian Polar Slice Sampling","text":"info: Info\nFor 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 toinitial_params has a sufficiently large norm. If this is suspected, simply try a different initialization value.","category":"page"},{"location":"gibbs_polar/#Demonstration","page":"Gibbsian Polar Slice Sampling","title":"Demonstration","text":"","category":"section"},{"location":"gibbs_polar/","page":"Gibbsian Polar Slice Sampling","title":"Gibbsian Polar Slice Sampling","text":"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):","category":"page"},{"location":"gibbs_polar/","page":"Gibbsian Polar Slice Sampling","title":"Gibbsian Polar Slice Sampling","text":"using Distributions\nusing Turing\nusing SliceSampling\nusing LinearAlgebra\nusing Plots\n\n@model function demo()\n x ~ MvTDist(1, zeros(10), Matrix(I,10,10))\nend\nmodel = demo()\n\nn_samples = 1000\nlatent_chain = sample(model, externalsampler(LatentSlice(10)), n_samples; initial_params=ones(10))\npolar_chain = sample(model, externalsampler(GibbsPolarSlice(10)), n_samples; initial_params=ones(10))\n\nl = @layout [a; b]\np1 = Plots.plot(1:n_samples, latent_chain[:,1,:], ylims=[-10,10], label=\"LSS\")\np2 = Plots.plot(1:n_samples, polar_chain[:,1,:], ylims=[-10,10], label=\"GPSS\")\nplot(p1, p2, layout = l)\nsavefig(\"student_latent_gpss.svg\")","category":"page"},{"location":"gibbs_polar/","page":"Gibbsian Polar Slice Sampling","title":"Gibbsian Polar Slice Sampling","text":"(Image: )","category":"page"},{"location":"gibbs_polar/","page":"Gibbsian Polar Slice Sampling","title":"Gibbsian Polar Slice Sampling","text":"Clearly, GPSS is better at exploring the deep tails compared to the latent slice sampler (LSS) despite having a similar per-iteration cost.","category":"page"},{"location":"gibbs_polar/","page":"Gibbsian Polar Slice Sampling","title":"Gibbsian Polar Slice Sampling","text":"[SHR2023]: Schär, P., Habeck, M., & Rudolf, D. (2023, July). Gibbsian polar slice sampling. In International Conference on Machine Learning.","category":"page"},{"location":"gibbs_polar/","page":"Gibbsian Polar Slice Sampling","title":"Gibbsian Polar Slice Sampling","text":"[RR2002]: Roberts, G. O., & Rosenthal, J. S. (2002). The polar slice sampler. Stochastic Models, 18(2), 257-280.","category":"page"},{"location":"meta_multivariate/#meta","page":"Meta Multivariate Samplers","title":"Meta Multivariate Samplers","text":"","category":"section"},{"location":"meta_multivariate/","page":"Meta Multivariate Samplers","title":"Meta Multivariate Samplers","text":"To use univariate slice sampling strategies on targets with more than on dimension, one has to embed them into a \"meta\" multivariate sampling scheme that relies on univariate sampling elements. The two most popular approaches for this are Gibbs sampling[GG1984] and hit-and-run[BRS1993].","category":"page"},{"location":"meta_multivariate/#Random-Permutation-Gibbs","page":"Meta Multivariate Samplers","title":"Random Permutation Gibbs","text":"","category":"section"},{"location":"meta_multivariate/","page":"Meta Multivariate Samplers","title":"Meta Multivariate Samplers","text":"Gibbs sampling[GG1984] is a strategy where we sample from the posterior one coordinate at a time, conditioned on the values of all other coordinates. In practice, one can pick the coordinates in any order they want as long as it does not depend on the state of the chain. It is generally hard to know a-prior which \"scan order\" is best, but randomly picking coordinates tend to work well in general. Currently, we only provide random permutation scan, which guarantees that all coordinates are updated at least once after d transitions. At the same time, reversibility is maintained by randomly permuting the order we go through each coordinate:","category":"page"},{"location":"meta_multivariate/","page":"Meta Multivariate Samplers","title":"Meta Multivariate Samplers","text":"RandPermGibbs","category":"page"},{"location":"meta_multivariate/#SliceSampling.RandPermGibbs","page":"Meta Multivariate Samplers","title":"SliceSampling.RandPermGibbs","text":"RandPermGibbs(unislice)\n\nRandom permutation coordinate-wise Gibbs sampling strategy. This applies unislice coordinate-wise in a random order.\n\nArguments\n\nunislice::Union{<:AbstractUnivariateSliceSampling,<:AbstractVector{<:AbstractUnivariateSliceSampling}}: a single or a vector of univariate slice sampling algorithms.\n\nWhen unislice is a vector of samplers, each slice sampler is applied to the corresponding coordinate of the target posterior. In that case, the length(unislice) must match the dimensionality of the posterior.\n\n\n\n\n\n","category":"type"},{"location":"meta_multivariate/","page":"Meta Multivariate Samplers","title":"Meta Multivariate Samplers","text":"Each call to AbstractMCMC.step internally performs d Gibbs transition so that all coordinates are updated.","category":"page"},{"location":"meta_multivariate/","page":"Meta Multivariate Samplers","title":"Meta Multivariate Samplers","text":"For example:","category":"page"},{"location":"meta_multivariate/","page":"Meta Multivariate Samplers","title":"Meta Multivariate Samplers","text":"RandPermGibbs(SliceSteppingOut(2.))","category":"page"},{"location":"meta_multivariate/","page":"Meta Multivariate Samplers","title":"Meta Multivariate Samplers","text":"If one wants to use a different slice sampler configuration for each coordinate, one can mix-and-match by passing a Vector of slice samplers, one for each coordinate. For instance, for a 2-dimensional target:","category":"page"},{"location":"meta_multivariate/","page":"Meta Multivariate Samplers","title":"Meta Multivariate Samplers","text":"RandPermGibbs([SliceSteppingOut(2.; max_proposals=32), SliceDoublingOut(2.),])","category":"page"},{"location":"meta_multivariate/#Hit-and-Run","page":"Meta Multivariate Samplers","title":"Hit-and-Run","text":"","category":"section"},{"location":"meta_multivariate/","page":"Meta Multivariate Samplers","title":"Meta Multivariate Samplers","text":"Hit-and-run is a simple meta algorithm where we sample over a random 1-dimensional projection of the space. That is, at each iteration, we sample a random direction","category":"page"},{"location":"meta_multivariate/","page":"Meta Multivariate Samplers","title":"Meta Multivariate Samplers","text":" theta_n sim operatornameUniform(mathbbS^d-1)","category":"page"},{"location":"meta_multivariate/","page":"Meta Multivariate Samplers","title":"Meta Multivariate Samplers","text":"and perform a Markov transition along the 1-dimensional subspace","category":"page"},{"location":"meta_multivariate/","page":"Meta Multivariate Samplers","title":"Meta Multivariate Samplers","text":"beginaligned\n lambda_n sim pleft(lambda mid x_n-1 theta_n right) propto pileft( x_n-1 + lambda theta_n right) \n x_n = x_n-1 + lambda_n theta_n\nendaligned","category":"page"},{"location":"meta_multivariate/","page":"Meta Multivariate Samplers","title":"Meta Multivariate Samplers","text":"where pi is the target unnormalized density. Applying slice sampling for the 1-dimensional subproblem has been popularized by David Mackay[M2003], and is, technically, also a Gibbs sampler. (Or is that Gibbs samplers are hit-and-run samplers?) Unlike RandPermGibbs, which only makes axis-aligned moves, HitAndRun can choose arbitrary directions, which could be helpful in some cases.","category":"page"},{"location":"meta_multivariate/","page":"Meta Multivariate Samplers","title":"Meta Multivariate Samplers","text":"HitAndRun","category":"page"},{"location":"meta_multivariate/#SliceSampling.HitAndRun","page":"Meta Multivariate Samplers","title":"SliceSampling.HitAndRun","text":"HitAndRun(unislice)\n\nHit-and-run sampling strategy[BRS1993]. This applies unislice along a random direction uniform sampled from the sphere.\n\nArguments\n\nunislice::AbstractUnivariateSliceSampling: Univariate slice sampling algorithm.\n\n\n\n\n\n","category":"type"},{"location":"meta_multivariate/","page":"Meta Multivariate Samplers","title":"Meta Multivariate Samplers","text":"This can be used, for example, as follows:","category":"page"},{"location":"meta_multivariate/","page":"Meta Multivariate Samplers","title":"Meta Multivariate Samplers","text":"HitAndRun(SliceSteppingOut(2.))","category":"page"},{"location":"meta_multivariate/","page":"Meta Multivariate Samplers","title":"Meta Multivariate Samplers","text":"Unlike RandPermGibbs, HitAndRun does not provide the option of using a unique unislice object for each coordinate. This is a natural limitation of the hit-and-run sampler: it does not operate on individual coordinates.","category":"page"},{"location":"meta_multivariate/","page":"Meta Multivariate Samplers","title":"Meta Multivariate Samplers","text":"[GG1984]: Geman, S., & Geman, D. (1984). Stochastic relaxation, Gibbs distributions, and the Bayesian restoration of images. IEEE Transactions on Pattern Analysis and Machine Intelligence, (6).","category":"page"},{"location":"meta_multivariate/","page":"Meta Multivariate Samplers","title":"Meta Multivariate Samplers","text":"[BRS1993]: Bélisle, C. J., Romeijn, H. E., & Smith, R. L. (1993). Hit-and-run algorithms for generating multivariate distributions. Mathematics of Operations Research, 18(2), 255-266.","category":"page"},{"location":"meta_multivariate/","page":"Meta Multivariate Samplers","title":"Meta Multivariate Samplers","text":"[M2003]: MacKay, D. J. (2003). Information theory, inference and learning algorithms. Cambridge university press.","category":"page"},{"location":"univariate_slice/#Univariate-Slice-Sampling-Algorithms","page":"Univariate Slice Sampling","title":"Univariate Slice Sampling Algorithms","text":"","category":"section"},{"location":"univariate_slice/#Introduction","page":"Univariate Slice Sampling","title":"Introduction","text":"","category":"section"},{"location":"univariate_slice/","page":"Univariate Slice Sampling","title":"Univariate Slice Sampling","text":"These algorithms are the \"single-variable\" slice sampling algorithms originally described by Neal[N2003]. Since these algorithms are univariate, one has to incorporate them into a \"meta\" multivariate sampler, which are discussed in this section.","category":"page"},{"location":"univariate_slice/#Fixed-Initial-Interval-Slice-Sampling","page":"Univariate Slice Sampling","title":"Fixed Initial Interval Slice Sampling","text":"","category":"section"},{"location":"univariate_slice/","page":"Univariate Slice Sampling","title":"Univariate Slice Sampling","text":"This is the most basic form of univariate slice sampling, where the proposals are generated within a fixed interval formed by the window.","category":"page"},{"location":"univariate_slice/","page":"Univariate Slice Sampling","title":"Univariate Slice Sampling","text":"Slice","category":"page"},{"location":"univariate_slice/#SliceSampling.Slice","page":"Univariate Slice Sampling","title":"SliceSampling.Slice","text":"Slice(window; max_proposals)\n\nUnivariate slice sampling with a fixed initial interval (Scheme 2 by Neal[N2003])\n\nArguments\n\nwindow::Real: Proposal window.\n\nKeyword Arguments\n\nmax_proposals::Int: Maximum number of proposals allowed until throwing an error (default: typemax(Int)).\n\n\n\n\n\n","category":"type"},{"location":"univariate_slice/#Adaptive-Initial-Interval-Slice-Sampling","page":"Univariate Slice Sampling","title":"Adaptive Initial Interval Slice Sampling","text":"","category":"section"},{"location":"univariate_slice/","page":"Univariate Slice Sampling","title":"Univariate Slice Sampling","text":"These algorithms try to adaptively set the initial interval through a simple search procedure. The \"stepping-out\" procedure grows the initial window on a linear scale, while the \"doubling-out\" procedure grows it geometrically. window controls the scale of the increase.","category":"page"},{"location":"univariate_slice/#What-Should-I-Use?","page":"Univariate Slice Sampling","title":"What Should I Use?","text":"","category":"section"},{"location":"univariate_slice/","page":"Univariate Slice Sampling","title":"Univariate Slice Sampling","text":"This highly depends on the problem at hand. In general, the doubling-out procedure tends to be more expensive as it requires additional log-target evaluations to decide whether to accept a proposal. However, if the scale of the posterior varies drastically, doubling out might work better. In general, it is recommended to use the stepping-out procedure.","category":"page"},{"location":"univariate_slice/","page":"Univariate Slice Sampling","title":"Univariate Slice Sampling","text":"SliceSteppingOut\nSliceDoublingOut","category":"page"},{"location":"univariate_slice/#SliceSampling.SliceSteppingOut","page":"Univariate Slice Sampling","title":"SliceSampling.SliceSteppingOut","text":"SliceSteppingOut(window; max_stepping_out, max_proposals)\n\nUnivariate slice sampling by automatically adapting the initial interval through the \"stepping-out\" procedure (Scheme 3 by Neal[N2003])\n\nArguments\n\nwindow::Real: Proposal window.\n\nKeyword Arguments\n\nmax_stepping_out::Int: Maximum number of \"stepping outs\" (default: 32).\nmax_proposals::Int: Maximum number of proposals allowed until throwing an error (default: typemax(Int)).\n\n\n\n\n\n","category":"type"},{"location":"univariate_slice/#SliceSampling.SliceDoublingOut","page":"Univariate Slice Sampling","title":"SliceSampling.SliceDoublingOut","text":"SliceDoublingOut(window; max_doubling_out, max_proposals)\n\nUnivariate slice sampling by automatically adapting the initial interval through the \"doubling-out\" procedure (Scheme 4 by Neal[N2003])\n\nArguments\n\nwindow::Real: Proposal window.\n\nKeyword Arguments\n\nmax_doubling_out: Maximum number of \"doubling outs\" (default: 8).\nmax_proposals::Int: Maximum number of proposals allowed until throwing an error (default: typemax(Int)).\n\n\n\n\n\n","category":"type"},{"location":"univariate_slice/","page":"Univariate Slice Sampling","title":"Univariate Slice Sampling","text":"[N2003]: Neal, R. M. (2003). Slice sampling. The Annals of Statistics, 31(3), 705-767.","category":"page"},{"location":"general/#General-Usage","page":"General Usage","title":"General Usage","text":"","category":"section"},{"location":"general/","page":"General Usage","title":"General Usage","text":"This package implements the AbstractMCMC interface. AbstractMCMC provides a unifying interface for MCMC algorithms applied to LogDensityProblems.","category":"page"},{"location":"general/#Examples","page":"General Usage","title":"Examples","text":"","category":"section"},{"location":"general/#Drawing-Samples-From-a-LogDensityProblems-Through-AbstractMCMC","page":"General Usage","title":"Drawing Samples From a LogDensityProblems Through AbstractMCMC","text":"","category":"section"},{"location":"general/","page":"General Usage","title":"General Usage","text":"SliceSampling.jl implements the AbstractMCMC interface through LogDensityProblems. That is, one simply needs to define a LogDensityProblems and pass it to AbstractMCMC:","category":"page"},{"location":"general/","page":"General Usage","title":"General Usage","text":"using AbstractMCMC\nusing Distributions\nusing LinearAlgebra\nusing LogDensityProblems\nusing Plots\n\nusing SliceSampling\n\nstruct Target{D}\n\tdist::D\nend\n\nLogDensityProblems.logdensity(target::Target, x) = logpdf(target.dist, x)\n\nLogDensityProblems.dimension(target::Target) = length(target.distx)\n\nLogDensityProblems.capabilities(::Type{<:Target}) = LogDensityProblems.LogDensityOrder{0}()\n\nsampler = GibbsPolarSlice(2.0)\nn_samples = 10000\nmodel = Target(MvTDist(5, zeros(10), Matrix(I, 10, 10)))\nlogdensitymodel = AbstractMCMC.LogDensityModel(model)\n\nchain = sample(logdensitymodel, sampler, n_samples; initial_params=randn(10))\nsamples = hcat([transition.params for transition in chain]...)\n\nplot(samples[1,:], xlabel=\"Iteration\", ylabel=\"Trace\")\nsavefig(\"abstractmcmc_demo.svg\")","category":"page"},{"location":"general/","page":"General Usage","title":"General Usage","text":"(Image: )","category":"page"},{"location":"general/#Drawing-Samples-From-Turing-Models","page":"General Usage","title":"Drawing Samples From Turing Models","text":"","category":"section"},{"location":"general/","page":"General Usage","title":"General Usage","text":"SliceSampling.jl can also be used to sample from Turing models through Turing's externalsampler interface:","category":"page"},{"location":"general/","page":"General Usage","title":"General Usage","text":"using Distributions\nusing Turing\nusing SliceSampling\n\n@model function demo()\n s ~ InverseGamma(3, 3)\n m ~ Normal(0, sqrt(s))\nend\n\nsampler = RandPermGibbs(SliceSteppingOut(2.))\nn_samples = 10000\nmodel = demo()\nsample(model, externalsampler(sampler), n_samples)","category":"page"},{"location":"general/#Drawing-Samples","page":"General Usage","title":"Drawing Samples","text":"","category":"section"},{"location":"general/","page":"General Usage","title":"General Usage","text":"For drawing samples using the algorithms provided by SliceSampling, the user only needs to call:","category":"page"},{"location":"general/","page":"General Usage","title":"General Usage","text":"sample([rng,] model, slice, N; initial_params)","category":"page"},{"location":"general/","page":"General Usage","title":"General Usage","text":"slice::AbstractSliceSampling: Any slice sampling algorithm provided by SliceSampling.\nmodel: A model implementing the LogDensityProblems interface.\nN: The number of samples","category":"page"},{"location":"general/","page":"General Usage","title":"General Usage","text":"The output is a SliceSampling.Transition object, which contains the following:","category":"page"},{"location":"general/","page":"General Usage","title":"General Usage","text":"SliceSampling.Transition","category":"page"},{"location":"general/#SliceSampling.Transition","page":"General Usage","title":"SliceSampling.Transition","text":"struct Transition\n\nStruct containing the results of the transition.\n\nFields\n\nparams: Samples generated by the transition.\nlp::Real: Log-target density of the samples.\ninfo::NamedTuple: Named tuple containing information about the transition. \n\n\n\n\n\n","category":"type"},{"location":"general/","page":"General Usage","title":"General Usage","text":"For the keyword arguments, SliceSampling allows:","category":"page"},{"location":"general/","page":"General Usage","title":"General Usage","text":"initial_params: The intial state of the Markov chain (default: nothing).","category":"page"},{"location":"general/","page":"General Usage","title":"General Usage","text":"If initial_params is nothing, the following function can be implemented to provide an initialization:","category":"page"},{"location":"general/","page":"General Usage","title":"General Usage","text":"SliceSampling.initial_sample","category":"page"},{"location":"general/#SliceSampling.initial_sample","page":"General Usage","title":"SliceSampling.initial_sample","text":"initial_sample(rng, model)\n\nReturn the initial sample for the model using the random number generator rng.\n\nArguments\n\nrng::Random.AbstractRNG: Random number generator.\nmodel: The target LogDensityProblem.\n\n\n\n\n\n","category":"function"},{"location":"general/#Performing-a-Single-Transition","page":"General Usage","title":"Performing a Single Transition","text":"","category":"section"},{"location":"general/","page":"General Usage","title":"General Usage","text":"For more fined-grained control, the user can call AbstractMCMC.step. That is, the chain can be initialized by calling:","category":"page"},{"location":"general/","page":"General Usage","title":"General Usage","text":"transition, state = AbstractMCMC.steps([rng,] model, slice; initial_params)","category":"page"},{"location":"general/","page":"General Usage","title":"General Usage","text":"and then each MCMC transition on state can be performed by calling:","category":"page"},{"location":"general/","page":"General Usage","title":"General Usage","text":"transition, state = AbstractMCMC.steps([rng,] model, slice, state)","category":"page"},{"location":"general/","page":"General Usage","title":"General Usage","text":"For more details, refer to the documentation of AbstractMCMC.","category":"page"},{"location":"","page":"Home","title":"Home","text":"CurrentModule = SliceSampling","category":"page"},{"location":"#SliceSampling","page":"Home","title":"SliceSampling","text":"","category":"section"},{"location":"","page":"Home","title":"Home","text":"This package implements slice sampling algorithms. Slice sampling finds its roots in the Swendsen–Wang algorithm for Ising models[SW1987][ES1988]. It later came into the interest of the statistical community through Besag and Green[BG1993], and popularized by Neal [N2003]. Furthermore, Neal introduced various ways to efficiently implement slice samplers. This package provides the original slice sampling algorithms by Neal and their later extensions.","category":"page"},{"location":"","page":"Home","title":"Home","text":"[SW1987]: Swendsen, R. H., & Wang, J. S. (1987). Nonuniversal critical dynamics in Monte Carlo simulations. Physical review letters, 58(2), 86.","category":"page"},{"location":"","page":"Home","title":"Home","text":"[ES1988]: Edwards, R. G., & Sokal, A. D. (1988). Generalization of the fortuin-kasteleyn-swendsen-wang representation and monte carlo algorithm. Physical review D, 38(6), 2009.","category":"page"},{"location":"","page":"Home","title":"Home","text":"[BG1993]: Besag, J., & Green, P. J. (1993). Spatial statistics and Bayesian computation. Journal of the Royal Statistical Society Series B: Statistical Methodology, 55(1), 25-37.","category":"page"},{"location":"","page":"Home","title":"Home","text":"[N2003]: Neal, R. M. (2003). Slice sampling. The annals of statistics, 31(3), 705-767.","category":"page"},{"location":"","page":"Home","title":"Home","text":"","category":"page"}] +[{"location":"latent_slice/#latent","page":"Latent Slice Sampling","title":"Latent Slice Sampling","text":"","category":"section"},{"location":"latent_slice/#Introduction","page":"Latent Slice Sampling","title":"Introduction","text":"","category":"section"},{"location":"latent_slice/","page":"Latent Slice Sampling","title":"Latent Slice Sampling","text":"Latent slice sampling is a recent vector-valued slice sampling algorithm proposed by Li and Walker[LW2023]. Unlike other slice sampling algorithms, it treats the \"search intervals\" as auxiliary variables and adapts them along the samples from the log-target in a Gibbs-type scheme.","category":"page"},{"location":"latent_slice/#Description","page":"Latent Slice Sampling","title":"Description","text":"","category":"section"},{"location":"latent_slice/","page":"Latent Slice Sampling","title":"Latent Slice Sampling","text":"Specifically, the extended joint density of the latent slice sampler is as follows:","category":"page"},{"location":"latent_slice/","page":"Latent Slice Sampling","title":"Latent Slice Sampling","text":" p(x t s l) = pi(x) p(s) operatornameUniformleft(t 0 pileft(xright)right) operatornameUniformleft(l x - s2 x + s2right)","category":"page"},{"location":"latent_slice/","page":"Latent Slice Sampling","title":"Latent Slice Sampling","text":"where y is the parameters of the log-target pi, s is the width of the search interval and l is the centering of the search interval relative to y. Naturally, the sampler operates as a blocked-Gibbs sampler ","category":"page"},{"location":"latent_slice/","page":"Latent Slice Sampling","title":"Latent Slice Sampling","text":"beginaligned\nl_n sim operatornameUniformleft(l x_n-1 - s_n-12 x_n-1 + s_n-12right) \ns_n sim p(s mid x_n-1 l_n) \nt_n sim operatornameUniformleft(0 pileft(x_n-1right)right) \nx_n sim operatornameUniformleftx mid pileft(xright) t_nright\nendaligned","category":"page"},{"location":"latent_slice/","page":"Latent Slice Sampling","title":"Latent Slice Sampling","text":"When x_n is updated using the usual shrinkage procedure of Neal[N2003], s_n and l_n are used to form the initial search window. (s_n is the width of the window and l_n is its center point.) Therefore, the latent slice sampler can be regarded as an automatic tuning mechanism of the \"initial interval\" of slice samplers.","category":"page"},{"location":"latent_slice/","page":"Latent Slice Sampling","title":"Latent Slice Sampling","text":"The only tunable parameter of the algorithm is then the distribution of the width p(s). For this, Li and Walker recommend","category":"page"},{"location":"latent_slice/","page":"Latent Slice Sampling","title":"Latent Slice Sampling","text":" p(s beta) = operatornameGamma(s 2 beta)","category":"page"},{"location":"latent_slice/","page":"Latent Slice Sampling","title":"Latent Slice Sampling","text":"where beta is a tunable parameter. The use of the gamma distribution is somewhat important since the complete conditional p(s mid y l) needs to be available in closed-form for efficiency. (It is a truncated exponential distribution in case of the gamma.) Therefore, we only provide control over beta.","category":"page"},{"location":"latent_slice/","page":"Latent Slice Sampling","title":"Latent Slice Sampling","text":"info: Info\nThe 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.","category":"page"},{"location":"latent_slice/#Interface","page":"Latent Slice Sampling","title":"Interface","text":"","category":"section"},{"location":"latent_slice/","page":"Latent Slice Sampling","title":"Latent Slice Sampling","text":"LatentSlice","category":"page"},{"location":"latent_slice/#SliceSampling.LatentSlice","page":"Latent Slice Sampling","title":"SliceSampling.LatentSlice","text":"LatentSlice(beta)\n\nLatent slice sampling algorithm by Li and Walker[LW2023].\n\nArguments\n\nbeta::Real: Beta parameter of the Gamma distribution of the auxiliary variables.\n\nKeyword Arguments\n\nmax_proposals::Int: Maximum number of proposals allowed until throwing an error (default: typemax(Int)).\n\n\n\n\n\n","category":"type"},{"location":"latent_slice/","page":"Latent Slice Sampling","title":"Latent Slice Sampling","text":"[LW2023]: Li, Y., & Walker, S. G. (2023). A latent slice sampling algorithm. Computational Statistics & Data Analysis, 179, 107652.","category":"page"},{"location":"latent_slice/","page":"Latent Slice Sampling","title":"Latent Slice Sampling","text":"[N2003]: Neal, R. M. (2003). Slice sampling. The annals of statistics, 31(3), 705-767.","category":"page"},{"location":"gibbs_polar/#polar","page":"Gibbsian Polar Slice Sampling","title":"Gibbsian Polar Slice Sampling","text":"","category":"section"},{"location":"gibbs_polar/#Introduction","page":"Gibbsian Polar Slice Sampling","title":"Introduction","text":"","category":"section"},{"location":"gibbs_polar/","page":"Gibbsian Polar Slice Sampling","title":"Gibbsian Polar Slice Sampling","text":"Gibbsian polar slice sampling (GPSS) is a recent vector-valued slice sampling algorithm proposed by P. Schär, M. Habeck, and D. Rudolf[SHR2023]. It is an computationally efficient variant of polar slice sampler previously proposed by Roberts and Rosenthal[RR2002]. Unlike other slice sampling algorithms, it operates a Gibbs sampler over polar coordinates, reminiscent of the elliptical slice sampler (ESS). 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.","category":"page"},{"location":"gibbs_polar/#Description","page":"Gibbsian Polar Slice Sampling","title":"Description","text":"","category":"section"},{"location":"gibbs_polar/","page":"Gibbsian Polar Slice Sampling","title":"Gibbsian Polar Slice Sampling","text":"For a d-dimensional target distribution pi, GPSS utilizes the following augmented target distribution:","category":"page"},{"location":"gibbs_polar/","page":"Gibbsian Polar Slice Sampling","title":"Gibbsian Polar Slice Sampling","text":"beginaligned\n p(x T) = varrho_pi^(0)(x) varrho_pi^(1)(x) operatornameUniformleft(T 0 varrho^1(x)right) \n varrho_pi^(0)(x) = lVert x rVert^1 - d \n varrho_pi^(1)(x) = lVert x rVert^d-1 pileft(xright)\nendaligned","category":"page"},{"location":"gibbs_polar/","page":"Gibbsian Polar Slice Sampling","title":"Gibbsian Polar Slice Sampling","text":"As described in Appendix A of the GPSS paper, sampling from varrho^(1)(x) in polar coordinates magically targets the augmented target distribution.","category":"page"},{"location":"gibbs_polar/","page":"Gibbsian Polar Slice Sampling","title":"Gibbsian Polar Slice Sampling","text":"In a high-level view, GPSS operates a Gibbs sampler in the following fashion:","category":"page"},{"location":"gibbs_polar/","page":"Gibbsian Polar Slice Sampling","title":"Gibbsian Polar Slice Sampling","text":"beginaligned\nT_n sim operatornameUniformleft(0 varrho^(1)left(x_n-1right)right) \ntheta_n sim operatornameUniformleft theta in mathbbS^d-1 mid varrho^(1)left(r_n-1 thetaright) T_n right \nr_n sim operatornameUniformleft r in mathbbR_geq 0 mid varrho^(1)left(r theta_nright) T_n right \nx_n = theta_n r_n\nendaligned","category":"page"},{"location":"gibbs_polar/","page":"Gibbsian Polar Slice Sampling","title":"Gibbsian Polar Slice Sampling","text":"where T_n is the usual acceptance threshold auxiliary variable, while theta and r are the sampler states in polar coordinates. The Gibbs steps on theta and r are implemented through specialized shrinkage procedures.","category":"page"},{"location":"gibbs_polar/","page":"Gibbsian Polar Slice Sampling","title":"Gibbsian Polar Slice Sampling","text":"The only tunable parameter of the algorithm is the size of the search interval (window) of the shrinkage sampler for the radius variable r.","category":"page"},{"location":"gibbs_polar/","page":"Gibbsian Polar Slice Sampling","title":"Gibbsian Polar Slice Sampling","text":"info: Info\nThe 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.","category":"page"},{"location":"gibbs_polar/#Interface","page":"Gibbsian Polar Slice Sampling","title":"Interface","text":"","category":"section"},{"location":"gibbs_polar/","page":"Gibbsian Polar Slice Sampling","title":"Gibbsian Polar Slice Sampling","text":"info: Info\nBy the nature of polar coordinates, GPSS only works reliably for targets with dimension at least d geq 2.","category":"page"},{"location":"gibbs_polar/","page":"Gibbsian Polar Slice Sampling","title":"Gibbsian Polar Slice Sampling","text":"GibbsPolarSlice","category":"page"},{"location":"gibbs_polar/#SliceSampling.GibbsPolarSlice","page":"Gibbsian Polar Slice Sampling","title":"SliceSampling.GibbsPolarSlice","text":"GibbsPolarSlice(w; max_proposals)\n\nGibbsian polar slice sampling algorithm by P. Schär, M. Habeck, and D. Rudolf [SHR2023].\n\nArguments\n\nw::Real: Initial window size for the radius shrinkage procedure.\n\nKeyword Arguments\n\nw::Real: Initial window size for the radius shrinkage procedure\nmax_proposals::Int: Maximum number of proposals allowed until throwing an error (default: typemax(Int)).\n\n\n\n\n\n","category":"type"},{"location":"gibbs_polar/","page":"Gibbsian Polar Slice Sampling","title":"Gibbsian Polar Slice Sampling","text":"warning: Warning\nWhen 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. ","category":"page"},{"location":"gibbs_polar/","page":"Gibbsian Polar Slice Sampling","title":"Gibbsian Polar Slice Sampling","text":"info: Info\nFor 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 toinitial_params has a sufficiently large norm. If this is suspected, simply try a different initialization value.","category":"page"},{"location":"gibbs_polar/#Demonstration","page":"Gibbsian Polar Slice Sampling","title":"Demonstration","text":"","category":"section"},{"location":"gibbs_polar/","page":"Gibbsian Polar Slice Sampling","title":"Gibbsian Polar Slice Sampling","text":"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):","category":"page"},{"location":"gibbs_polar/","page":"Gibbsian Polar Slice Sampling","title":"Gibbsian Polar Slice Sampling","text":"using Distributions\nusing Turing\nusing SliceSampling\nusing LinearAlgebra\nusing Plots\n\n@model function demo()\n x ~ MvTDist(1, zeros(10), Matrix(I,10,10))\nend\nmodel = demo()\n\nn_samples = 1000\nlatent_chain = sample(model, externalsampler(LatentSlice(10)), n_samples; initial_params=ones(10))\npolar_chain = sample(model, externalsampler(GibbsPolarSlice(10)), n_samples; initial_params=ones(10))\n\nl = @layout [a; b]\np1 = Plots.plot(1:n_samples, latent_chain[:,1,:], ylims=[-10,10], label=\"LSS\")\np2 = Plots.plot(1:n_samples, polar_chain[:,1,:], ylims=[-10,10], label=\"GPSS\")\nplot(p1, p2, layout = l)\nsavefig(\"student_latent_gpss.svg\")","category":"page"},{"location":"gibbs_polar/","page":"Gibbsian Polar Slice Sampling","title":"Gibbsian Polar Slice Sampling","text":"(Image: )","category":"page"},{"location":"gibbs_polar/","page":"Gibbsian Polar Slice Sampling","title":"Gibbsian Polar Slice Sampling","text":"Clearly, GPSS is better at exploring the deep tails compared to the latent slice sampler (LSS) despite having a similar per-iteration cost.","category":"page"},{"location":"gibbs_polar/","page":"Gibbsian Polar Slice Sampling","title":"Gibbsian Polar Slice Sampling","text":"[SHR2023]: Schär, P., Habeck, M., & Rudolf, D. (2023, July). Gibbsian polar slice sampling. In International Conference on Machine Learning.","category":"page"},{"location":"gibbs_polar/","page":"Gibbsian Polar Slice Sampling","title":"Gibbsian Polar Slice Sampling","text":"[RR2002]: Roberts, G. O., & Rosenthal, J. S. (2002). The polar slice sampler. Stochastic Models, 18(2), 257-280.","category":"page"},{"location":"meta_multivariate/#meta","page":"Meta Multivariate Samplers","title":"Meta Multivariate Samplers","text":"","category":"section"},{"location":"meta_multivariate/","page":"Meta Multivariate Samplers","title":"Meta Multivariate Samplers","text":"To use univariate slice sampling strategies on targets with more than on dimension, one has to embed them into a \"meta\" multivariate sampling scheme that relies on univariate sampling elements. The two most popular approaches for this are Gibbs sampling[GG1984] and hit-and-run[BRS1993].","category":"page"},{"location":"meta_multivariate/#Random-Permutation-Gibbs","page":"Meta Multivariate Samplers","title":"Random Permutation Gibbs","text":"","category":"section"},{"location":"meta_multivariate/","page":"Meta Multivariate Samplers","title":"Meta Multivariate Samplers","text":"Gibbs sampling[GG1984] is a strategy where we sample from the posterior one coordinate at a time, conditioned on the values of all other coordinates. In practice, one can pick the coordinates in any order they want as long as it does not depend on the state of the chain. It is generally hard to know a-prior which \"scan order\" is best, but randomly picking coordinates tend to work well in general. Currently, we only provide random permutation scan, which guarantees that all coordinates are updated at least once after d transitions. At the same time, reversibility is maintained by randomly permuting the order we go through each coordinate:","category":"page"},{"location":"meta_multivariate/","page":"Meta Multivariate Samplers","title":"Meta Multivariate Samplers","text":"RandPermGibbs","category":"page"},{"location":"meta_multivariate/#SliceSampling.RandPermGibbs","page":"Meta Multivariate Samplers","title":"SliceSampling.RandPermGibbs","text":"RandPermGibbs(unislice)\n\nRandom permutation coordinate-wise Gibbs sampling strategy. This applies unislice coordinate-wise in a random order.\n\nArguments\n\nunislice::Union{<:AbstractUnivariateSliceSampling,<:AbstractVector{<:AbstractUnivariateSliceSampling}}: a single or a vector of univariate slice sampling algorithms.\n\nWhen unislice is a vector of samplers, each slice sampler is applied to the corresponding coordinate of the target posterior. In that case, the length(unislice) must match the dimensionality of the posterior.\n\n\n\n\n\n","category":"type"},{"location":"meta_multivariate/","page":"Meta Multivariate Samplers","title":"Meta Multivariate Samplers","text":"Each call to AbstractMCMC.step internally performs d Gibbs transition so that all coordinates are updated.","category":"page"},{"location":"meta_multivariate/","page":"Meta Multivariate Samplers","title":"Meta Multivariate Samplers","text":"For example:","category":"page"},{"location":"meta_multivariate/","page":"Meta Multivariate Samplers","title":"Meta Multivariate Samplers","text":"RandPermGibbs(SliceSteppingOut(2.))","category":"page"},{"location":"meta_multivariate/","page":"Meta Multivariate Samplers","title":"Meta Multivariate Samplers","text":"If one wants to use a different slice sampler configuration for each coordinate, one can mix-and-match by passing a Vector of slice samplers, one for each coordinate. For instance, for a 2-dimensional target:","category":"page"},{"location":"meta_multivariate/","page":"Meta Multivariate Samplers","title":"Meta Multivariate Samplers","text":"RandPermGibbs([SliceSteppingOut(2.; max_proposals=32), SliceDoublingOut(2.),])","category":"page"},{"location":"meta_multivariate/#Hit-and-Run","page":"Meta Multivariate Samplers","title":"Hit-and-Run","text":"","category":"section"},{"location":"meta_multivariate/","page":"Meta Multivariate Samplers","title":"Meta Multivariate Samplers","text":"Hit-and-run is a simple meta algorithm where we sample over a random 1-dimensional projection of the space. That is, at each iteration, we sample a random direction","category":"page"},{"location":"meta_multivariate/","page":"Meta Multivariate Samplers","title":"Meta Multivariate Samplers","text":" theta_n sim operatornameUniform(mathbbS^d-1)","category":"page"},{"location":"meta_multivariate/","page":"Meta Multivariate Samplers","title":"Meta Multivariate Samplers","text":"and perform a Markov transition along the 1-dimensional subspace","category":"page"},{"location":"meta_multivariate/","page":"Meta Multivariate Samplers","title":"Meta Multivariate Samplers","text":"beginaligned\n lambda_n sim pleft(lambda mid x_n-1 theta_n right) propto pileft( x_n-1 + lambda theta_n right) \n x_n = x_n-1 + lambda_n theta_n\nendaligned","category":"page"},{"location":"meta_multivariate/","page":"Meta Multivariate Samplers","title":"Meta Multivariate Samplers","text":"where pi is the target unnormalized density. Applying slice sampling for the 1-dimensional subproblem has been popularized by David Mackay[M2003], and is, technically, also a Gibbs sampler. (Or is that Gibbs samplers are hit-and-run samplers?) Unlike RandPermGibbs, which only makes axis-aligned moves, HitAndRun can choose arbitrary directions, which could be helpful in some cases.","category":"page"},{"location":"meta_multivariate/","page":"Meta Multivariate Samplers","title":"Meta Multivariate Samplers","text":"HitAndRun","category":"page"},{"location":"meta_multivariate/#SliceSampling.HitAndRun","page":"Meta Multivariate Samplers","title":"SliceSampling.HitAndRun","text":"HitAndRun(unislice)\n\nHit-and-run sampling strategy[BRS1993]. This applies unislice along a random direction uniform sampled from the sphere.\n\nArguments\n\nunislice::AbstractUnivariateSliceSampling: Univariate slice sampling algorithm.\n\n\n\n\n\n","category":"type"},{"location":"meta_multivariate/","page":"Meta Multivariate Samplers","title":"Meta Multivariate Samplers","text":"This can be used, for example, as follows:","category":"page"},{"location":"meta_multivariate/","page":"Meta Multivariate Samplers","title":"Meta Multivariate Samplers","text":"HitAndRun(SliceSteppingOut(2.))","category":"page"},{"location":"meta_multivariate/","page":"Meta Multivariate Samplers","title":"Meta Multivariate Samplers","text":"Unlike RandPermGibbs, HitAndRun does not provide the option of using a unique unislice object for each coordinate. This is a natural limitation of the hit-and-run sampler: it does not operate on individual coordinates.","category":"page"},{"location":"meta_multivariate/","page":"Meta Multivariate Samplers","title":"Meta Multivariate Samplers","text":"[GG1984]: Geman, S., & Geman, D. (1984). Stochastic relaxation, Gibbs distributions, and the Bayesian restoration of images. IEEE Transactions on Pattern Analysis and Machine Intelligence, (6).","category":"page"},{"location":"meta_multivariate/","page":"Meta Multivariate Samplers","title":"Meta Multivariate Samplers","text":"[BRS1993]: Bélisle, C. J., Romeijn, H. E., & Smith, R. L. (1993). Hit-and-run algorithms for generating multivariate distributions. Mathematics of Operations Research, 18(2), 255-266.","category":"page"},{"location":"meta_multivariate/","page":"Meta Multivariate Samplers","title":"Meta Multivariate Samplers","text":"[M2003]: MacKay, D. J. (2003). Information theory, inference and learning algorithms. Cambridge university press.","category":"page"},{"location":"univariate_slice/#Univariate-Slice-Sampling-Algorithms","page":"Univariate Slice Sampling","title":"Univariate Slice Sampling Algorithms","text":"","category":"section"},{"location":"univariate_slice/#Introduction","page":"Univariate Slice Sampling","title":"Introduction","text":"","category":"section"},{"location":"univariate_slice/","page":"Univariate Slice Sampling","title":"Univariate Slice Sampling","text":"These algorithms are the \"single-variable\" slice sampling algorithms originally described by Neal[N2003]. Since these algorithms are univariate, one has to incorporate them into a \"meta\" multivariate sampler, which are discussed in this section.","category":"page"},{"location":"univariate_slice/#Fixed-Initial-Interval-Slice-Sampling","page":"Univariate Slice Sampling","title":"Fixed Initial Interval Slice Sampling","text":"","category":"section"},{"location":"univariate_slice/","page":"Univariate Slice Sampling","title":"Univariate Slice Sampling","text":"This is the most basic form of univariate slice sampling, where the proposals are generated within a fixed interval formed by the window.","category":"page"},{"location":"univariate_slice/","page":"Univariate Slice Sampling","title":"Univariate Slice Sampling","text":"Slice","category":"page"},{"location":"univariate_slice/#SliceSampling.Slice","page":"Univariate Slice Sampling","title":"SliceSampling.Slice","text":"Slice(window; max_proposals)\n\nUnivariate slice sampling with a fixed initial interval (Scheme 2 by Neal[N2003])\n\nArguments\n\nwindow::Real: Proposal window.\n\nKeyword Arguments\n\nmax_proposals::Int: Maximum number of proposals allowed until throwing an error (default: typemax(Int)).\n\n\n\n\n\n","category":"type"},{"location":"univariate_slice/#Adaptive-Initial-Interval-Slice-Sampling","page":"Univariate Slice Sampling","title":"Adaptive Initial Interval Slice Sampling","text":"","category":"section"},{"location":"univariate_slice/","page":"Univariate Slice Sampling","title":"Univariate Slice Sampling","text":"These algorithms try to adaptively set the initial interval through a simple search procedure. The \"stepping-out\" procedure grows the initial window on a linear scale, while the \"doubling-out\" procedure grows it geometrically. window controls the scale of the increase.","category":"page"},{"location":"univariate_slice/#What-Should-I-Use?","page":"Univariate Slice Sampling","title":"What Should I Use?","text":"","category":"section"},{"location":"univariate_slice/","page":"Univariate Slice Sampling","title":"Univariate Slice Sampling","text":"This highly depends on the problem at hand. In general, the doubling-out procedure tends to be more expensive as it requires additional log-target evaluations to decide whether to accept a proposal. However, if the scale of the posterior varies drastically, doubling out might work better. In general, it is recommended to use the stepping-out procedure.","category":"page"},{"location":"univariate_slice/","page":"Univariate Slice Sampling","title":"Univariate Slice Sampling","text":"SliceSteppingOut\nSliceDoublingOut","category":"page"},{"location":"univariate_slice/#SliceSampling.SliceSteppingOut","page":"Univariate Slice Sampling","title":"SliceSampling.SliceSteppingOut","text":"SliceSteppingOut(window; max_stepping_out, max_proposals)\n\nUnivariate slice sampling by automatically adapting the initial interval through the \"stepping-out\" procedure (Scheme 3 by Neal[N2003])\n\nArguments\n\nwindow::Real: Proposal window.\n\nKeyword Arguments\n\nmax_stepping_out::Int: Maximum number of \"stepping outs\" (default: 32).\nmax_proposals::Int: Maximum number of proposals allowed until throwing an error (default: typemax(Int)).\n\n\n\n\n\n","category":"type"},{"location":"univariate_slice/#SliceSampling.SliceDoublingOut","page":"Univariate Slice Sampling","title":"SliceSampling.SliceDoublingOut","text":"SliceDoublingOut(window; max_doubling_out, max_proposals)\n\nUnivariate slice sampling by automatically adapting the initial interval through the \"doubling-out\" procedure (Scheme 4 by Neal[N2003])\n\nArguments\n\nwindow::Real: Proposal window.\n\nKeyword Arguments\n\nmax_doubling_out: Maximum number of \"doubling outs\" (default: 8).\nmax_proposals::Int: Maximum number of proposals allowed until throwing an error (default: typemax(Int)).\n\n\n\n\n\n","category":"type"},{"location":"univariate_slice/","page":"Univariate Slice Sampling","title":"Univariate Slice Sampling","text":"[N2003]: Neal, R. M. (2003). Slice sampling. The Annals of Statistics, 31(3), 705-767.","category":"page"},{"location":"general/#General-Usage","page":"General Usage","title":"General Usage","text":"","category":"section"},{"location":"general/","page":"General Usage","title":"General Usage","text":"This package implements the AbstractMCMC interface. AbstractMCMC provides a unifying interface for MCMC algorithms applied to LogDensityProblems.","category":"page"},{"location":"general/#Examples","page":"General Usage","title":"Examples","text":"","category":"section"},{"location":"general/#Drawing-Samples-From-a-LogDensityProblems-Through-AbstractMCMC","page":"General Usage","title":"Drawing Samples From a LogDensityProblems Through AbstractMCMC","text":"","category":"section"},{"location":"general/","page":"General Usage","title":"General Usage","text":"SliceSampling.jl implements the AbstractMCMC interface through LogDensityProblems. That is, one simply needs to define a LogDensityProblems and pass it to AbstractMCMC:","category":"page"},{"location":"general/","page":"General Usage","title":"General Usage","text":"using AbstractMCMC\nusing Distributions\nusing LinearAlgebra\nusing LogDensityProblems\nusing Plots\n\nusing SliceSampling\n\nstruct Target{D}\n\tdist::D\nend\n\nLogDensityProblems.logdensity(target::Target, x) = logpdf(target.dist, x)\n\nLogDensityProblems.dimension(target::Target) = length(target.distx)\n\nLogDensityProblems.capabilities(::Type{<:Target}) = LogDensityProblems.LogDensityOrder{0}()\n\nsampler = GibbsPolarSlice(2.0)\nn_samples = 10000\nmodel = Target(MvTDist(5, zeros(10), Matrix(I, 10, 10)))\nlogdensitymodel = AbstractMCMC.LogDensityModel(model)\n\nchain = sample(logdensitymodel, sampler, n_samples; initial_params=randn(10))\nsamples = hcat([transition.params for transition in chain]...)\n\nplot(samples[1,:], xlabel=\"Iteration\", ylabel=\"Trace\")\nsavefig(\"abstractmcmc_demo.svg\")","category":"page"},{"location":"general/","page":"General Usage","title":"General Usage","text":"(Image: )","category":"page"},{"location":"general/#Drawing-Samples-From-Turing-Models","page":"General Usage","title":"Drawing Samples From Turing Models","text":"","category":"section"},{"location":"general/","page":"General Usage","title":"General Usage","text":"SliceSampling.jl can also be used to sample from Turing models through Turing's externalsampler interface:","category":"page"},{"location":"general/","page":"General Usage","title":"General Usage","text":"using Distributions\nusing Turing\nusing SliceSampling\n\n@model function demo()\n s ~ InverseGamma(3, 3)\n m ~ Normal(0, sqrt(s))\nend\n\nsampler = RandPermGibbs(SliceSteppingOut(2.))\nn_samples = 10000\nmodel = demo()\nsample(model, externalsampler(sampler), n_samples)","category":"page"},{"location":"general/#Conditional-sampling-in-a-Turing.Experimental.Gibbs-sampler","page":"General Usage","title":"Conditional sampling in a Turing.Experimental.Gibbs sampler","text":"","category":"section"},{"location":"general/","page":"General Usage","title":"General Usage","text":"SliceSampling.jl be used as a conditional sampler in Turing.Experimental.Gibbs.","category":"page"},{"location":"general/","page":"General Usage","title":"General Usage","text":"using Distributions\nusing Turing\nusing SliceSampling\n\n@model function simple_choice(xs)\n p ~ Beta(2, 2)\n z ~ Bernoulli(p)\n for i in 1:length(xs)\n if z == 1\n xs[i] ~ Normal(0, 1)\n else\n xs[i] ~ Normal(2, 1)\n end\n end\nend\n\nsampler = Turing.Experimental.Gibbs(\n (\n p = externalsampler(SliceSteppingOut(2.0)),\n z = PG(20, :z)\n )\n)\n\nn_samples = 1000\nmodel = simple_choice([1.5, 2.0, 0.3])\nsample(model, sampler, n_samples)","category":"page"},{"location":"general/#Drawing-Samples","page":"General Usage","title":"Drawing Samples","text":"","category":"section"},{"location":"general/","page":"General Usage","title":"General Usage","text":"For drawing samples using the algorithms provided by SliceSampling, the user only needs to call:","category":"page"},{"location":"general/","page":"General Usage","title":"General Usage","text":"sample([rng,] model, slice, N; initial_params)","category":"page"},{"location":"general/","page":"General Usage","title":"General Usage","text":"slice::AbstractSliceSampling: Any slice sampling algorithm provided by SliceSampling.\nmodel: A model implementing the LogDensityProblems interface.\nN: The number of samples","category":"page"},{"location":"general/","page":"General Usage","title":"General Usage","text":"The output is a SliceSampling.Transition object, which contains the following:","category":"page"},{"location":"general/","page":"General Usage","title":"General Usage","text":"SliceSampling.Transition","category":"page"},{"location":"general/#SliceSampling.Transition","page":"General Usage","title":"SliceSampling.Transition","text":"struct Transition\n\nStruct containing the results of the transition.\n\nFields\n\nparams: Samples generated by the transition.\nlp::Real: Log-target density of the samples.\ninfo::NamedTuple: Named tuple containing information about the transition. \n\n\n\n\n\n","category":"type"},{"location":"general/","page":"General Usage","title":"General Usage","text":"For the keyword arguments, SliceSampling allows:","category":"page"},{"location":"general/","page":"General Usage","title":"General Usage","text":"initial_params: The intial state of the Markov chain (default: nothing).","category":"page"},{"location":"general/","page":"General Usage","title":"General Usage","text":"If initial_params is nothing, the following function can be implemented to provide an initialization:","category":"page"},{"location":"general/","page":"General Usage","title":"General Usage","text":"SliceSampling.initial_sample","category":"page"},{"location":"general/#SliceSampling.initial_sample","page":"General Usage","title":"SliceSampling.initial_sample","text":"initial_sample(rng, model)\n\nReturn the initial sample for the model using the random number generator rng.\n\nArguments\n\nrng::Random.AbstractRNG: Random number generator.\nmodel: The target LogDensityProblem.\n\n\n\n\n\n","category":"function"},{"location":"general/#Performing-a-Single-Transition","page":"General Usage","title":"Performing a Single Transition","text":"","category":"section"},{"location":"general/","page":"General Usage","title":"General Usage","text":"For more fined-grained control, the user can call AbstractMCMC.step. That is, the chain can be initialized by calling:","category":"page"},{"location":"general/","page":"General Usage","title":"General Usage","text":"transition, state = AbstractMCMC.steps([rng,] model, slice; initial_params)","category":"page"},{"location":"general/","page":"General Usage","title":"General Usage","text":"and then each MCMC transition on state can be performed by calling:","category":"page"},{"location":"general/","page":"General Usage","title":"General Usage","text":"transition, state = AbstractMCMC.steps([rng,] model, slice, state)","category":"page"},{"location":"general/","page":"General Usage","title":"General Usage","text":"For more details, refer to the documentation of AbstractMCMC.","category":"page"},{"location":"","page":"Home","title":"Home","text":"CurrentModule = SliceSampling","category":"page"},{"location":"#SliceSampling","page":"Home","title":"SliceSampling","text":"","category":"section"},{"location":"","page":"Home","title":"Home","text":"This package implements slice sampling algorithms. Slice sampling finds its roots in the Swendsen–Wang algorithm for Ising models[SW1987][ES1988]. It later came into the interest of the statistical community through Besag and Green[BG1993], and popularized by Neal [N2003]. Furthermore, Neal introduced various ways to efficiently implement slice samplers. This package provides the original slice sampling algorithms by Neal and their later extensions.","category":"page"},{"location":"","page":"Home","title":"Home","text":"[SW1987]: Swendsen, R. H., & Wang, J. S. (1987). Nonuniversal critical dynamics in Monte Carlo simulations. Physical review letters, 58(2), 86.","category":"page"},{"location":"","page":"Home","title":"Home","text":"[ES1988]: Edwards, R. G., & Sokal, A. D. (1988). Generalization of the fortuin-kasteleyn-swendsen-wang representation and monte carlo algorithm. Physical review D, 38(6), 2009.","category":"page"},{"location":"","page":"Home","title":"Home","text":"[BG1993]: Besag, J., & Green, P. J. (1993). Spatial statistics and Bayesian computation. Journal of the Royal Statistical Society Series B: Statistical Methodology, 55(1), 25-37.","category":"page"},{"location":"","page":"Home","title":"Home","text":"[N2003]: Neal, R. M. (2003). Slice sampling. The annals of statistics, 31(3), 705-767.","category":"page"},{"location":"","page":"Home","title":"Home","text":"","category":"page"}] } diff --git a/dev/student_latent_gpss.svg b/dev/student_latent_gpss.svg index 456734d..dd18dfc 100644 --- a/dev/student_latent_gpss.svg +++ b/dev/student_latent_gpss.svg @@ -1,78 +1,78 @@ - + - + - + - + - + - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + - + - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/dev/univariate_slice/index.html b/dev/univariate_slice/index.html index 93b50c2..fe81f87 100644 --- a/dev/univariate_slice/index.html +++ b/dev/univariate_slice/index.html @@ -1,2 +1,2 @@ -Univariate Slice Sampling · SliceSampling.jl

Univariate Slice Sampling Algorithms

Introduction

These algorithms are the "single-variable" slice sampling algorithms originally described by Neal[N2003]. Since these algorithms are univariate, one has to incorporate them into a "meta" multivariate sampler, which are discussed in this section.

Fixed Initial Interval Slice Sampling

This is the most basic form of univariate slice sampling, where the proposals are generated within a fixed interval formed by the window.

SliceSampling.Slice — Type
Slice(window; max_proposals)

Univariate slice sampling with a fixed initial interval (Scheme 2 by Neal[N2003])

Arguments

  • window::Real: Proposal window.

Keyword Arguments

  • max_proposals::Int: Maximum number of proposals allowed until throwing an error (default: typemax(Int)).
source

Adaptive Initial Interval Slice Sampling

These algorithms try to adaptively set the initial interval through a simple search procedure. The "stepping-out" procedure grows the initial window on a linear scale, while the "doubling-out" procedure grows it geometrically. window controls the scale of the increase.

What Should I Use?

This highly depends on the problem at hand. In general, the doubling-out procedure tends to be more expensive as it requires additional log-target evaluations to decide whether to accept a proposal. However, if the scale of the posterior varies drastically, doubling out might work better. In general, it is recommended to use the stepping-out procedure.

SliceSampling.SliceSteppingOut — Type
SliceSteppingOut(window; max_stepping_out, max_proposals)

Univariate slice sampling by automatically adapting the initial interval through the "stepping-out" procedure (Scheme 3 by Neal[N2003])

Arguments

  • window::Real: Proposal window.

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)).
source
SliceSampling.SliceDoublingOut — Type
SliceDoublingOut(window; max_doubling_out, max_proposals)

Univariate slice sampling by automatically adapting the initial interval through the "doubling-out" procedure (Scheme 4 by Neal[N2003])

Arguments

  • window::Real: Proposal window.

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)).
source
  • N2003Neal, R. M. (2003). Slice sampling. The Annals of Statistics, 31(3), 705-767.
+Univariate Slice Sampling · SliceSampling.jl

Univariate Slice Sampling Algorithms

Introduction

These algorithms are the "single-variable" slice sampling algorithms originally described by Neal[N2003]. Since these algorithms are univariate, one has to incorporate them into a "meta" multivariate sampler, which are discussed in this section.

Fixed Initial Interval Slice Sampling

This is the most basic form of univariate slice sampling, where the proposals are generated within a fixed interval formed by the window.

SliceSampling.Slice — Type
Slice(window; max_proposals)

Univariate slice sampling with a fixed initial interval (Scheme 2 by Neal[N2003])

Arguments

  • window::Real: Proposal window.

Keyword Arguments

  • max_proposals::Int: Maximum number of proposals allowed until throwing an error (default: typemax(Int)).
source

Adaptive Initial Interval Slice Sampling

These algorithms try to adaptively set the initial interval through a simple search procedure. The "stepping-out" procedure grows the initial window on a linear scale, while the "doubling-out" procedure grows it geometrically. window controls the scale of the increase.

What Should I Use?

This highly depends on the problem at hand. In general, the doubling-out procedure tends to be more expensive as it requires additional log-target evaluations to decide whether to accept a proposal. However, if the scale of the posterior varies drastically, doubling out might work better. In general, it is recommended to use the stepping-out procedure.

SliceSampling.SliceSteppingOut — Type
SliceSteppingOut(window; max_stepping_out, max_proposals)

Univariate slice sampling by automatically adapting the initial interval through the "stepping-out" procedure (Scheme 3 by Neal[N2003])

Arguments

  • window::Real: Proposal window.

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)).
source
SliceSampling.SliceDoublingOut — Type
SliceDoublingOut(window; max_doubling_out, max_proposals)

Univariate slice sampling by automatically adapting the initial interval through the "doubling-out" procedure (Scheme 4 by Neal[N2003])

Arguments

  • window::Real: Proposal window.

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)).
source
  • N2003Neal, R. M. (2003). Slice sampling. The Annals of Statistics, 31(3), 705-767.