---
title: "`vws`: Vertical Weighted Strips in R using C++"
author: Andrew M. Raim
format:
  pdf:
    mathspec: true
    fontsize: 11pt
    indent: false
    toc: true
    number-sections: true
    colorlinks: true
    link-citations: true
    prompt: false
    fig-pos: H
    include-in-header:
      text: |
        \usepackage{common}
    cite-method: natbib
    biblio-style: plainnat
    natbiboptions:
      - round
pdf-engine: pdflatex
vignette: >
  %\VignetteIndexEntry{vws}
  %\VignetteEncoding{UTF-8}
  %\VignetteEngine{quarto::pdf}
bibliography: references.bib
editor_options: 
  chunk_output_type: console
classoption: abstract
abstract: |
  The `vws` package facilitates the vertical weighted strips (VWS) method in R.
  VWS is an approach to construct proposals for rejection sampling which is
  explored in @vws-2026. The `vws` package provides an API to program the
  necessary components for the proposal and carry out rejection sampling. The
  API of the `vws` package is in C++ and makes use of templates and
  object-oriented programming; a working understanding of these techniques is
  assumed. The primary intended usage is that sampler logic is implemented in
  C++ as a variate generation function; this is exposed in R for use in
  applications. This vignette presents the API and several complete examples;
  multiple variations of VWS samplers are presented with each example. The
  `vws` package and this vignette have been written primarily for univariate
  targets; however, it is intended that the framework can be used to handle
  multivariate problems as well.
thanks: |
  **For correspondence**`:` <andrew.raim@gmail.com>.
  Center for Statistical Research & Methodology, U.S. Census Bureau,
  Washington, DC, 20233, U.S.A.
  Document was compiled `{r} format(Sys.time(), "%Y-%m-%d %H:%M:%S %Z")`.
geometry:
  - left=0.75in
  - right=0.75in
  - top=0.75in
  - bottom=1.00in
execute:
  eval: true
code-block-bg: "#FAFAFA"
code-block-border-left: "#B2DEED"
filters:
  - include-code-files
callout-icon: false
---

```{r}
#| include: false
library(vws)
library(tidyverse)

set.seed(1234)
```

# Disclaimer and Acknowledgments {.unnumbered}

This document is released to inform interested parties of ongoing research and to encourage discussion of work in progress. Any views expressed are those of the author and not of the U.S. Census Bureau.

The author is grateful to Drs. James Livsey and Kyle Irimata for discussions that motivated development of this package. Thanks to Drs. Kyle Irimata and Tommy Wright for their reviews of this manuscript.

Although there are no guarantees of correctness of the `vws` package, reasonable efforts will be made to address shortcomings. Comments, questions, corrections, and possible improvements can be communicated through the project’s Github repository (<https://github.com/andrewraim/vws>).

# Introduction

The `vws` package facilitates vertical weighted strips (VWS) sampling [@vws-2026]. VWS is an approach to construct proposal distributions for rejection sampling. Here, the target density is regarded as a weighted density - the product of a nonnegative weight function and a base density. The weight function is majorized by a given function (which bounds it above) to yield an envelope which is recombined with the base density to form a proposal distribution. If it is convenient to generate draws from the proposal, and the distance between the original and majorized weight function is not too large, such a proposal can be effectively used in rejection sampling. When the weight function is majorized in a piecewise manner over the support, the resulting proposal can be regarded as a finite mixture. Such a mixture can be refined to reduce the rejection probability to an acceptable tolerance. The term "vertical weighted strips" emphasizes that it is an extension of a method referred to as "vertical strips" in Section 3.6.1 of @Martino-Luengo-Miguez-2018, and was previously described in Chapter VIII of @Devroye-1986. Partitioning in VWS is done to facilitate majorization of the weight function rather than to majorize the entire target density as in vertical strips.

The `vws` package provides tools to construct and utilize VWS proposal distributions. Variate generation functions are constructed in C++ and may be exposed for use in R via the Rcpp package [@Eddelbuettel-2103]. C++ is taken as the primary programming language because of its efficiency as a compiled language and its formal approach to object-orientation and template programming. The `fntl` package [@fntl] is integral to `vws` usage; it provides a simplified C++ interface to useful routines in the R application programming interface (API) which is described in @R-ext, as well as other useful numerical routines, where functions are supplied as lambdas. The `vws` package has been designed primarily to generate from univariate distributions (both continuous and discrete). It is also possible to construct proposals for multivariate distributions using the VWS approach; the `vws` package is intended to support code development for such settings as well.

This vignette presents an overview of the package, its API, and several examples of working samplers. An `R>` prompt is shown in some code displays to emphasize interaction via the console. Examples make liberal use of the `tidyverse` [@Tidyverse-2019], especially `ggplot2` to plot and `dplyr` to manipulate tables. Some of the lengthier codes for the examples are not presented in this document; however, complete codes are provided with the package. They can be accessed in the `doc/examples` folder of the installed `vws` package. This path can be accessed in R with the following command.

```{r}
#| eval: false
#| prompt: true
file.path(path.package("vws"), "doc", "examples")
```

The remainder of the vignette proceeds as follows. @sec-vws briefly reviews VWS sampling. @sec-overview gives an overview of the package, its major components, and important preliminaries for programming. @sec-api provides a detailed specification of the C++ API. @sec-vmf develops several VWS samplers to generate from the Von Mises Fisher (VMF) distribution following the application in @vws-2026. @sec-privacy develops samplers in a disclosure avoidance setting relevant to the U.S. Census Bureau, which is considered by @Direct-Sampling-DAS-2021, @DP-Simulation-2022, and @DAS-Methods-2025. @sec-bessel develops VWS samplers for the Bessel count distribution [@Devroye-2002]. In some applications, it may be desirable to have a persistent VWS proposal object in R that can gradually be refined and reused to draw samples; @sec-persistent presents several mechanisms to achieve this. @sec-privacy-within-gibbs extends the example from @sec-privacy to a Gibbs sampler and demonstrates VWS-within-Gibbs. The first variant uses persistent proposals, as in @sec-persistent, with a Gibbs sampler written in R and VWS implemented in C++. The second variant implements both in C++.

::: {.callout-tip title="Quick Start"}
Readers looking for a quick start can go directly to the first example in @sec-vmf-const-default whose contents are described with the highest amount of detail.
:::

# A Brief Review of Vertical Weighted Strips {#sec-vws}

The objective of VWS is to sample from a weighted target density
<!-- -->
\begin{align}
f(x) = f_0(x) / \psi, \quad
f_0(x) = w(x) g(x), \quad
\psi = \int_\Omega f_0(x) d\nu(x),
\label{eqn:weighted-target}
\end{align}
<!-- -->
where $\Omega$ is the support, $\nu$ is a dominating measure, $g$ is assumed to be a normalized density, $w(x)$ is a nonnegative weight function, and $\psi$ is a normalizing constant. We will construct a proposal of the form
<!-- -->
\begin{align*}
h(x) = h_0(x) / \psi_N, \quad
h_0(x) = \overline{w}(x) g(x), \quad
\psi_N = \int_\Omega h_0(x) d\nu(x).
\end{align*}
<!-- -->
We will say that a function $\phi$ *majorizes* $w$ on $A \subseteq \Omega$ if $\phi(x) \geq w(x)$ for $x \in A$. Let $\ind\{ x \in A \}$ denote the indicator function for $x \in A$. Suppose $\Omega$ is partitioned into regions $\mathscr{D}_1, \ldots, \mathscr{D}_N$ and there are corresponding functions $\overline{w}_1, \ldots, \overline{w}_N$ where $\overline{w}_j$ majorizes $w$ on $\mathscr{D}_j$. Then $\overline{w}(x) = \sum_{j=1}^N \overline{w}_j(x) \ind\{x \in \mathscr{D}_j\}$ majorizes $w$ on $\Omega$ and the unnormalized proposal becomes
<!-- -->
\begin{align*}
h_0(x) = g(x) \sum_{j=1}^N \overline{w}_j(x) \ind\{x \in \mathscr{D}_j\}.
\end{align*}
<!-- -->
With this construction, $f_0(x) \leq h_0(x)$ for all $x \in \Omega$. Therefore, classical rejection sampling can be carried out by drawing $u$ from $\text{Uniform}(0,1)$, $x$ from $h$, and accepting $x$ as a draw from $f$ if $u \leq f_0(x) / h_0(x)$. A benefit of this construction is that the functional form of the density can help to guide proposal selection. The normalized $h$ can be obtained by defining
<!-- -->
\begin{align*}
\overline{\xi}_j = \int_{\mathscr{D}_j} \overline{w}_j(x) g(x) d\nu(x)
\end{align*}
<!-- -->
and
<!-- -->
\begin{align*}
\psi_N = \sum_{j=1}^N \overline{\xi}_j,
\end{align*}
<!-- -->
giving the finite mixture
<!-- -->
\begin{align}
h(x)
= h_0(x) / \psi_N
= \sum_{j=1}^N \pi_j g_j(x),
\label{eqn:fmm-proposal}
\end{align}
<!-- -->
a finite mixture with mixing weights and component densities,
<!-- -->
\begin{align*}
\pi_j = \frac{\overline{\xi}_j}{\sum_{\ell=1}^N \overline{\xi}_\ell},
\quad \text{and} \quad
g_j(x) = 
\overline{w}_j(x) g(x) \ind\{x \in \mathscr{D}_j\} / \overline{\xi}_j,
\end{align*}
<!-- -->
respectively. The $g_j$ are truncated and reweighted versions of base distribution $g$. In addition to the majorizer, suppose that $\underline{w}_j$ is a *minorizer* of $w$ so that $0 \leq \underline{w}_j(x) \leq w(x)$ for all $x \in \mathscr{D}_j$, and let 
<!-- -->
\begin{math}
\underline{\xi}_j = \int_{\mathscr{D}_j} \underline{w}_j(x) g(x) d\nu(x).
\end{math}
<!-- -->
Note that the majorizer and minorizer need not have the same functional form. When $h$ is used as a proposal in rejection sampling, an upper bound for the probability of rejection is
<!-- -->
\begin{align}
\rho_+ = \sum_{j=1}^N \rho_j,
\quad \text{where} \quad
\rho_j = \frac{
\overline{\xi}_j - \underline{\xi}_j
}{
\sum_{j=1}^N \overline{\xi}_\ell
}.
\label{eqn:bound}
\end{align}
<!-- -->
The quantities $\rho_1, \ldots, \rho_N$ are the contributions of each region to $\rho_+$. The bound $\rho_+$ can be used to determine whether a VWS proposal will be viable for rejection sampling; if it is not much less than 1, rejections may be too frequent for practical use. The proposal may be refined by altering the partitioning or considering a different majorizer. Several specific choices of majorizer are considered by @vws-2026 and will be reviewed in the remainder of this section. @sec-vws-constant discusses the use of a constant majorizing function. A linear majorizing function is discussed in @sec-vws-linear. @sec-vws-knots reviews several knot selection methods.

## Constant Majorizer {#sec-vws-constant}

Suppose $\Omega = (\alpha_0, \alpha_N]$ is an interval whose endpoints may or may not be finite and $w(x) < \infty$ on $\Omega$. We can majorize $w$ using a constant $\overline{w}_j \geq \sup_{x \in \mathscr{D}_j} w(x)$ on each $\mathscr{D}_j$ so that the majorizer is $\overline{w}(x) = \sum_{j=1}^N \overline{w}_j \ind\{ x \in \mathscr{D}_j \}$. Here we obtain component densities of finite mixture \eqref{eqn:fmm-proposal} as
<!-- -->
\begin{align*}
g_j(x) = g(x) \ind\{x \in \mathscr{D}_j\} / \Prob(T \in \mathscr{D}_j),
\quad \text{with} \quad
\overline{\xi}_j = \overline{w}_j \Prob(T \in \mathscr{D}_j).
\end{align*}
<!-- -->
Furthermore, if we assume a constant minorizer $\underline{w}_j \leq \inf_{x \in \mathscr{D}_j} w(x)$, then $\underline{\xi}_j = \underline{w}_j \Prob(T \in \mathscr{D}_j)$.

## Linear Majorizer {#sec-vws-linear}

Suppose $\Omega = (\alpha_0, \alpha_N]$ is an interval whose endpoints may or may not be finite and $0 < w(x) < \infty$ on $\Omega$. Furthermore, suppose $\Omega$ is partitioned into regions $\mathscr{D}_j = (\alpha_{j-1}, \alpha_j]$ where $w$ is entirely either log-convex or log-concave for each $j = 1\, \ldots, N$. An exponentiated linear function $\overline{w}_j(x) = \exp\{ \overline{\beta}_{0j} + \overline{\beta}_{1j} x \}$ can be used to majorize $w$ on $\mathscr{D}_j$ with an appropriate choice of coefficients $\overline{\beta}_{0j}$ and $\overline{\beta}_{1j}$. Therefore, the piecewise linear function $\overline{w}(x) = \sum_{j=1}^N \overline{w}_j(x) \ind\{ x \in \mathscr{D}_j \}$ is a majorizer on $\Omega$. Here the component densities of finite mixture \eqref{eqn:fmm-proposal} are
<!-- -->
\begin{align*}
g_j(x) = \exp\{ \overline{\beta}_{0j} + \overline{\beta}_{1j} x \} g(x) \ind\{ x \in \mathscr{D}_j \} / \overline{\xi}_j,
%
\quad \text{with} \quad
%
\overline{\xi}_j = e^{ \overline{\beta}_{0j} } \int_{\alpha_{j-1}}^{\alpha_j}
e^{ \overline{\beta}_{1j} \cdot x } g(x) d\nu(x).
\end{align*}
<!-- -->
In certain cases, such as where $g$ is in an exponential family, the integral in $\overline{\xi}_j$ can be simplified and $g_j$ is seen to be a familiar distribution. We may also assume a linear minorizer $\overline{w}(x) = \sum_{j=1}^N \underline{w}_j(x) \ind\{ x \in \mathscr{D}_j \}$ with $\overline{w}_j(x) = \exp\{ \underline{\beta}_{0j} + \underline{\beta}_{1j} x \}$ for appropriate coefficients $\underline{\beta}_{0j}$ and $\underline{\beta}_{1j}$.

To identify coefficients $\overline{\beta}_{0j}$, $\overline{\beta}_{1j}$, $\underline{\beta}_{0j}$, and $\underline{\beta}_{1j}$, consider the following remarks corresponding to the log-concave and log-convex cases, respectively.

::: {#rem-linear-concave}
Suppose $w$ is log-concave on $\mathscr{D}_j$. For any $c \in \mathscr{D}_j$,
<!-- -->
\begin{align}
\log w(x)
&\leq \log w(c) + (x - c) \nabla(c) \nonumber \\
&= \overline{\beta}_{0j} + \overline{\beta}_{1j} \cdot x,
\label{eqn:linear-concave-majorizer}
\end{align}
<!-- -->
so that
\begin{math}
\overline{w}_j(x) = \exp\{ \overline{\beta}_{0j} + \overline{\beta}_{1j} \cdot x \}
\end{math}
<!-- -->
with
<!-- -->
\begin{math}
\overline{\beta}_{0j} = \log w(c) - c \cdot \nabla(c),
\end{math}
<!-- -->
\begin{math}
\overline{\beta}_{1j} = \nabla(c),
\end{math}
<!-- -->
and
<!-- -->
\begin{math}
\nabla(x) = \frac{d}{dx} \log w(x).
\end{math}
<!-- -->
For an accompanying minorizer, let $\lambda \in [0,1]$ for a given $x \in \mathscr{D}_j$ so that $x = (1 - \lambda) \alpha_{j-1} + \lambda \alpha_j$. Concavity of $\log w(x)$ gives
<!-- -->
\begin{align}
\log w(x) &\geq (1-\lambda) \log w(\alpha_{j-1}) + \lambda \log w(\alpha_j) \nonumber \\
&= \log w(\alpha_{j-1}) + \frac{x - \alpha_{j-1}}{\alpha_j - \alpha_{j-1}} [ \log w(\alpha_j) - \log w(\alpha_{j-1})] \nonumber \\
&= \underline{\beta}_{0j} + \underline{\beta}_{1j} \cdot x,
\label{eqn:linear-concave-minorizer}
\end{align}
<!-- -->
so that
<!-- -->
\begin{math}
\underline{w}_j(x) = \exp\{ \underline{\beta}_{0j} + \underline{\beta}_{1j} \cdot x \}
\end{math}
<!-- -->
with
<!-- -->
\begin{math}
\underline{\beta}_{0j} = \log w(\alpha_{j-1}) - \alpha_{j-1} \underline{\beta}_{1j}
\end{math}
<!-- -->
and
<!-- -->
\begin{math}
\underline{\beta}_{1j} = \{\log w(\alpha_j) - \log w(\alpha_{j-1}) \} / \{ \alpha_j - \alpha_{j-1} \}.
\end{math}
\remarkqed
:::

::: {#rem-linear-convex}
Suppose $w$ is log-convex on $\mathscr{D}_j$. A majorizer and minorizer are obtained from \eqref{eqn:linear-concave-minorizer} and \eqref{eqn:linear-concave-majorizer}, respectively. Note that they have switched roles from @rem-linear-concave.
\remarkqed
:::

The following remark describes one possible choice for the expansion point $c$ which will be utilized in the examples in Sections [-@sec-vmf-linear], [-@sec-privacy-linear], and [-@sec-bessel-linear].

::: {#rem-expansion-point}
For a majorizer in the log-concave case, we choose $c$ to minimize the L1 distance between unnormalized densities,
<!-- -->
\begin{align}
c^*
&= \argmin_{c \in \mathscr{D}_j} \int_{\mathscr{D}_j} [\overline{w}_j(x) - w(x) ] g(x) d\nu(x) \nonumber \\
%
&= \argmin_{c \in \mathscr{D}_j} \Big\{ w(c) \exp\{ -c \nabla(c) \} \int_{\mathscr{D}_j} \exp\{ x \nabla(c) \} g(x) d\nu(x) \Big\} \nonumber \\
%
&= \argmin_{c \in \mathscr{D}_j} \Big\{ \log w(c) - c \nabla(c) + \log M_j(\nabla(c)) \Big\},
\label{eqn:concave-majorize-choice}
\end{align}
<!-- -->
where 
<!-- -->
\begin{math}
\overline{w}_j(x) = \exp\{ \overline{\beta}_{0j} + \overline{\beta}_{1j} \cdot x \},
\end{math}
<!-- -->
\begin{math}
\overline{\beta}_{0j} = \log w(c) - c \cdot \nabla(c),
\end{math}
<!-- -->
\begin{math}
\overline{\beta}_{1j} = \nabla(c),
\end{math}
<!-- -->
\begin{math}
\nabla(x) = \frac{d}{dx} \log w(x),
\end{math}
<!-- -->
and
<!-- -->
\begin{math}
M_j(s) = \int_{\alpha_{j-1}}^{\alpha_j} e^{xs} g(x) d\nu(x).
\end{math}
<!-- -->
A similar choice of $c$ for a minorizer in the log-convex case is
<!-- -->
\begin{align}
c^* &= \argmin_{c \in \mathscr{D}_j} \int_{\mathscr{D}_j} [w(x) - \underline{w}_j(x) ] g(x) d\nu(x) \nonumber \\
%
&= \argmax_{c \in \mathscr{D}_j} \int_{\mathscr{D}_j} \underline{w}_j(x) g(x) d\nu(x) \nonumber \\
%
&= \argmax_{c \in \mathscr{D}_j} \Big\{ \log w(c) - c \nabla(c) + \log M_j(\nabla(c)) \Big\}.
\label{eqn:convex-minorize-choice}
\end{align}
<!-- -->
where 
<!-- -->
\begin{math}
\overline{w}_j(x) = \exp\{ \underline{\beta}_{0j} + \underline{\beta}_{1j} \cdot x \},
\end{math}
<!-- -->
\begin{math}
\underline{\beta}_{0j} = \log w(c) - c \cdot \nabla(c),
\end{math}
<!-- -->
\begin{math}
\underline{\beta}_{1j} = \nabla(c),
\end{math}
<!-- -->
\begin{math}
\nabla(x) = \frac{d}{dx} \log w(x).
\end{math}
<!-- -->
\remarkqed
:::

## Knot Selection {#sec-vws-knots}

After selecting a decomposition of the target into weighted form \eqref{eqn:weighted-target}, the choice of knots is important to obtain an efficient proposal distribution. One option is to specify knots at known locations. For example, we can consider using evenly spaced points on a given interval within the support; however, this may not take into account important features of the target and add wasteful components to the finite mixture which do not contribute much to the quality of the approximation. An automated option that often achieves much better results is the rule of thumb proposed by @vws-2026, shown here as Algorithm \ref{alg:refine}. This method sequentially refines an initial partition of the support consisting of one or more regions. At each step a region is selected and bifurcated at its midpoint. The selection is probabilistic, with probabilities proportional to $\rho_1, \ldots, \rho_N$. A deterministic ("greedy") variation can be considered by instead always selecting a region $\ell$ whose value $\rho_\ell$ is largest.

Bifurcation at a midpoint is mostly straightforward in univariate settings where the support is a subset of the real line. However, there may be multiple ways to characterize this operation when the support is multivariate, perhaps with more structure like a multinomial sample space. The necessary operations can be coded for such settings within the `vws` framework; however, this document does not consider such cases in detail.

\begin{algorithm}
\caption{Probabilistic rule of thumb for sequential knot selection.}
\label{alg:refine}
\hspace*{\algorithmicindent}\textbf{Input}: maximum number of knots to add $N$. \\
\hspace*{\algorithmicindent}\textbf{Input}: initial vector of internal knots $\alpha_1, \ldots, \alpha_{N_0-1}$; may be empty with $N_0 = 0$. \\
\hspace*{\algorithmicindent}\textbf{Input}: tolerance $\epsilon > 0$.
\begin{algorithmic}[1]
\State $j \leftarrow 0$ 
\While{$j \leq N$}
\State Let $\mathscr{D}_\ell$ with $\rho_\ell$ for $\ell \in \{ 1, \ldots, N_0 + j \}$ be current regions.
\State If $\sum_{\ell=1}^{N_0 + j} \rho_\ell < \epsilon$, break from the loop.
\State Draw $\ell \in \{ 1, \ldots N_0 + j \}$ from $\text{Discrete}(\rho_1, \ldots, \rho_{N_0 + j})$. \label{line:selection}
\State Let $\alpha^*$ be midpoint of $\mathscr{D}_\ell$; add $\alpha^*$ to vector of knots.
\State Let $j \leftarrow j + 1$. \label{line:midpoint}
\EndWhile
\State \Return $(\alpha_0, \ldots, \alpha_{N_0 + j})$.
\end{algorithmic}
\end{algorithm}

## Self Tuning {#sec-self-tune}

VWS may be used to generate from conditional distributions in a Gibbs sampler; however, the amount of computation can become arduous if new proposals are constructed in each step of the Gibbs sampler. In this setting, @saevws-2025 consider a persistent VWS proposal that is gradually tuned to track a particular conditional distribution in the Gibbs sampler as the MCMC chain evolves. If the target does not change much as the chain approaches the posterior distribution, the amount of tuning needed to maintain a rejection rate which is sufficiently below 1 should also diminish.

Two types of updates are considered when a proposed draw $x$ is rejected. If bound $\rho_+$ is larger than a given threshold $\epsilon_1 \in (0,1)$, $x$ is used to refine the current set of regions. The region $\mathscr{D}_j$ containing $x$ is selected for bifurcation, but the specific way that $x$ is used to bifurcate depends on the nature of the regions. For example, there are a number of possibilities to bifurcate a multivariate region at point $x$. On the other hand, if $\rho_+$ is smaller than $\epsilon_1$, we check for regions with small contributions of $\rho_\ell$ to $\rho_+$. In particular regions with $\rho_\ell$ smaller than a given threshold $\epsilon_2 \in (0,1)$ are coarsened (i.e., merged with a neighboring region) if doing so does not increase $\rho_+$ beyond $\epsilon_1$. This process is repeated until no further merges can be completed. 

To emphasize, $\epsilon_1$ represents a threshold for a sufficiently small rejection rate while $\epsilon_2$ represents a tolerance for low contribution. In practice, a sufficiently small rejection rate need not be far from 1 if rejections are computationally inexpensive; a value such as $\epsilon_1 = 0.85$ may yield good performance.

Algorithm \ref{alg:self-tune} represents rejection sampling with self-tuning described in the manner above. It makes use of Algorithm \ref{alg:merge} to handle merges. Note that only certain pairs of regions are eligible to be merged, and eligibility is determined by the particular type of region. For example, the regions implemented in @sec-api-realconstregion and @sec-api-intconstregion are based on univariate intervals; here, only adjacent intervals having a common endpoint can be merged together. There may be multiple merges possible with a given partition, where different choices on line  \ref{step:merge} of Algorithm \ref{alg:merge} yield different configurations after the loop on line \ref{step:merge-loop} of Algorithm \ref{alg:self-tune}. For the sake of simplicity, the current implementation performs the first eligible merge in each traversal of the loop rather than considering all possible merges.

@saevws-2025 present several examples where tuning can effectively be limited to an initial part of the MCMC chain and discontinued for the remainder of the chain. Here the number of rejections per iteration can be monitored to ensure that proposals remain adequate as the chain progresses. If rejections increase dramatically after tuning is halted, a longer tuning period may be needed.

\begin{algorithm}
\caption{Generate a single draw from $f$ using self-tuned VWS with tolerances $\epsilon_1$ and $\epsilon_2$.}
\label{alg:self-tune}
\begin{algorithmic}[1]
\Function{RejectionTune}{$h$ by reference, $\epsilon_1$, $\epsilon_2$}
\Do
\State Draw $u$ from $\text{Uniform}(0,1)$ and $x$ from $h$.
\State $\text{Accept} \gets \ind\{ u \leq f_0(x) / h_0(x) \}$ where $f_0(x) = w(x) g(x)$ and $h_0(x) = \overline{w}(x) g(x)$.
\State Let $\rho_{+}$ be the current value of bound \eqref{eqn:bound}.
\If{$\text{Accept} = 0$ and $\rho_{+} < \epsilon_1$}
\Do
\label{step:merge-loop}
\State $\text{Repeat} \leftarrow \text{MergeOnce}(h, \epsilon_1, \epsilon_2)$.
\doWhile{$\text{Repeat} = 1$}
\ElsIf{$\text{Accept} = 0$ and $\rho_{+} \geq \epsilon_1$}
\State Refine $h$ at $x$.
\EndIf
\doWhile{$\text{Accept} = 0$}
\State \Return $x$ and $h$.
\EndFunction
\end{algorithmic}
\end{algorithm}

\begin{algorithm}
\caption{Attempt to merge a single pair of regions in $h$. Return $1$ if a merge was completed and $0$ otherwise.}
\label{alg:merge}
\begin{algorithmic}[1]
\Function{MergeOnce}{$h$ by reference, $\epsilon_1$, $\epsilon_2$}
\State Let $\mathscr{D}_1, \ldots, \mathscr{D}_N$ be the regions corresponding to current mixture components of $h$.
\State Let $\rho_{+}$ be the current value of bound \eqref{eqn:bound} and $\rho_1, \ldots, \rho_N$ be contribution of each region.
\For{$j = 1, \ldots, N$}
\If{$\rho_j < \epsilon_2$}
\State Let $\ell$ ($\neq j$) be the index of a region that can be merged with $\mathscr{D}_j$.
\label{step:merge}
\State Let $\tilde{h}$ be equal to $h$ with regions $\mathscr{D}_j$ and $\mathscr{D}_\ell$ merged.
\State Let $\tilde{\rho}_+$ be the bound \eqref{eqn:bound} for $\tilde{h}$.
\If{$\tilde{\rho}_+ < \epsilon_1$}
\State Modify $h$ to $\tilde{h}$.
\State \Return $1$.
\EndIf
\EndIf
\EndFor
\State \Return $0$.
\EndFunction
\end{algorithmic}
\end{algorithm}

# Overview of Package {#sec-overview}

The `vws` package aims to support the methodology which was described in the previous section. The present section will describe tools in the package which can be used to formulate a problem, construct a proposal, and generate samples. Some useful preliminaries are also presented.

## Design {#sec-overview-design}

Use of the `vws` package centers on two classes. The `fmm_proposal` class represents a finite mixture of the form \eqref{eqn:fmm-proposal} and encapsulates the operations needed for rejection sampling. The `Region` class represents region $\mathscr{D}_j$ and the operations required for VWS. `Region` itself is an abstract base class which defines an interface needed by the framework. All problem-specific logic is coded within a subclass of `Region`: either by using one of the provided subclasses or by coding a new subclass. The subclasses `real_const_region` and `int_const_region` provide implementations of the constant majorizer from @sec-vws-constant, for continuous and discrete univariate distributions, respectively. The linear majorizer described in @sec-vws-linear requires more customization than the constant majorizer; here, users should create a subclass of `Region`. An `fmm_proposal` object is created from $N \geq 1$ `Region` objects that represent the partition $\mathscr{D}_1, \ldots, \mathscr{D}_N$ of $\Omega$.

The `rejection` function takes an object `h` of class `fmm_proposal` and carries out rejection sampling to obtain $n$ draws. The `rejection` function returns a vector of accepted draws along with information about the number of rejections. @fig-software-design displays a diagram of the high-level design just described.

![Design of `vws` package.](software-design.pdf){#fig-software-design}

## Typical Workflow {#sec-overview-workflow}

The following display outlines a typical workflow for VWS rejection sampling in C++. Some of the lines are labeled with numbers; descriptions of these labels are given after the code.

```{cpp}
// [[Rcpp::depends(vws, fntl)]]                  // <1>
#include "vws.h"                                 // <2>
#include "MyRegion.h"

// [[Rcpp::export]]                              // <3>
Rcpp::List sample(unsigned int n)
{
    MyRegion supp( ... );                        // <4>
    vws::fmm_proposal<double, MyRegion> h(supp); // <5>
    h.refine(N - 1, 0.01);                       // <6>
    auto out = vws::rejection(h, n);             // <7>

    return Rcpp::List::create(                   // <8>
        Rcpp::Named("draws") = out.draws,
        Rcpp::Named("rejects") = out.rejects
    );
}
```

1. This annotation ensures that `vws` and `fntl` packages are linked during compilation. `fntl` is not explicitly used in this example, but it is needed for linking by the `vws` library regardless.
2. Include the main header for the `vws` package so that the API is accessible.
3. This annotation exports the `sample` function for use in R. This is not needed if the `sample` function will be used only in C++.
4. Construct a region of type `MyRegion` which we have presumably defined in `MyRegion.h`. The symbol `...` here is a placeholder for any arguments needed by the constructor.
5. Construct an `fmm_proposal` based on a single region `supp`. Note that it is a template class that needs two template arguments: the first (`double`) indicates the data type of the distribution and the second (`MyRegion`) indicates the type of region it will contain.
6. Refine the proposal using Algorithm \ref{alg:refine} up to $N-1$ times or until tolerance `0.01` is achieved for \eqref{eqn:bound}.
7. Carry out rejection sampling with the refined proposal `h`.
8. Pack elements of the struct `out` into an `Rcpp::List` and return them to the caller.

Note that the `fmm_proposal` class, the `rejection` function, and other elements of the `vws` API are accessed in the `vws` namespace. Details for the API are given in @sec-api. Once the `sample` function is successfully exported as an R function, we may invoke it as usual. Suppose the code in the previous display is saved is in a file named `sample.cpp`.

```{r}
#| eval: false
#| prompt: true
Rcpp::sourceCpp("sample.cpp")
out = sample(n = 100)
head(out$draws)
head(out$rejects)
```

::: {.callout-tip title="Working Examples"}
A complete VWS example is given in @sec-vmf-const-default; this demonstrates the use of `real_const_region` to implement a sampler with a minimal amount of coding. @sec-vmf-const-custom is an improvement which requires more coding; here we provide functions to determine values of the constant majorizer and minorizer analytically. Finally, @sec-vmf-linear demonstrates how to subclass `Region` to implement a linear majorizer, which requires more intricate coding. Subsequent examples are given in @sec-privacy and @sec-bessel, but with less explanation.
:::

## Lambdas {#sec-overview-lambdas}

Lambdas in C++ are functions which can be defined within the course of a program, regarded as objects which can be passed to other functions, and capture variables in the environment so that they do not need to be passed as arguments. They should feel very familiar to R programmers who likely use R functions in this way. Here is a snippet with several examples of lambdas.

```{cpp}
double z = 3;
typedef function<double(double)> dfdd;
dfdd f1 = [&](double x, double y) -> double { return x*y*z; };
dfdd f2 = [&](double x, double y) { return x*y*z; };
dfdd f3 = [=](double x, double y) { return x*y*z; };
dfdd f4 = [](double x, double y) { return x*y*3; };
```

These four functions carry out the same operation and may be invoked in the usual way.

```{cpp}
double out = f1(1, 2);
```

The typedef `dfdd` is given as a shorthand for a Standard Template Library (STL) [`function`][stl-func] which takes two `double` variables as arguments and returns a `double`. An STL `function` can be passed to other functions like any other variable. Function `f1` explicitly states the return type of the lambda using `-> double` while it is omitted in the others to be inferred by the compiler. The variable `z` is *captured* by `f1`, `f2`, and `f3`; in other words, it is "baked in" to the function definition without needing to be passed as an argument. Functions `f1` and `f2` capture `z` by reference while `f3` captures `z` by value, indicated with the syntaxes `[&]` and `[=]` respectively. Function `f4` uses empty square brackets `[]` to indicate that it does not capture any outside variables.

[stl-func]: https://en.cppreference.com/cpp/utility/functional/function

The `fntl` package defines several function typedefs used in the `vws` package, and also provides access to numerical tools such as integration, differentiation, and optimization which take STL `function`s as arguments. \remarkqed

## Computations on the Log-Scale {#sec-logscale}

Many of the computations in this vignette and within the `vws` package are carried out on the logarithmic scale. This is to avoid loss of precision when working with very large or small magnitude numbers. Several functions to support log-scale arithmetic are given in @sec-api-logscale. \remarkqed

## Optimization {#sec-overview-optimization}

To help automate computation of constant majorizers and minorizers for use with `real_const_region` and `int_const_region`, the default is to use numerical optimization to find appropriate constants on each region in the partition. The particular optimization method is presented in @sec-api-opt. This can be somewhat of a computational burden; it can sometimes be avoided when there is a closed form solution or other informed way to arrive at maximizing / minimizing constants. The API for `real_const_region` and `int_const_region` accepts user-specified lambdas to specify alternative methods. \remarkqed

# C++ API {#sec-api}

This section documents functions, classes, and other components in the C++ API. Note that some functions in the C++ API are also exposed as R functions; see the `vws` package manual for these.

## Typedefs and Aliases {#sec-api-typedefs}

The following typedef is used to represent functions that take `double` and `bool` arguments and return a `double`. It is used as the type for weight functions in particular, with `x` the argument of the weight function and `log = true` if the result is to be returned on the log-scale.

```{cpp}
typedef std::function<double(double x, bool log)> dfdb;
```

The following typedef represents a function that optimizes (i.e., maximizes or minimizes) a weight function on a given interval.

```{cpp}
typedef std::function<double(
	const dfdb& w,  // <1>
	double lo,      // <2>
	double hi,      // <3>
	bool log        // <4>
)> optimizer;
```

1. A weight function.
2. Lower bound of the interval; may be `R_NegInf`.
3. Upper bound of the interval; may be `R_PosInf`.
4. Logical; `log = true` specifies that the result should be the optimal value $\log w(x^*)$, given on the log-scale. Otherwise, the result should be the optimal value $w(x^*)$ on the original scale.

The following typedef represents a function that computes the midpoint of a given interval. Note that this type name can clash with the `midpoint` function in the STL; therefore, it is referenced using its namespace `vws::midpoint` within the `vws` package.

```{cpp}
typedef std::function<double(double a, double b)> midpoint;
```

The initial release of the `vws` package used a "camelcase" convention for class names, which has since changed to lowercase with underscores. The following aliases are defined for backward compatibility.

```{cpp}
template<class T, class R>
using FMMProposal = fmm_proposal<T,R>;

template<class T>
using Region = region<T>;

using RealConstRegion = real_const_region;
using IntConstRegion = int_const_region;
using UnivariateHelper = univariate_helper;
```

## Finite Mixture Proposal {#sec-api-proposal}

The class `fmm_proposal` represents a VWS proposal in the form of \eqref{eqn:fmm-proposal}.

### Class Definition {.unlisted}

`fmm_proposal` has two template arguments: `T` is the data type for the underlying distribution (e.g., `T = double` for a real-valued univariate random variable) and `R` is the type of region which the proposal will be based upon.

```{cpp}
template <class T, class R>
class fmm_proposal { ... };
```

### Constructor {.unlisted}

The primary constructor takes a vector of regions of type `R`. There must be at least one region in the vector. These regions are expected to be a partition of the support $\Omega$.

```{cpp}
fmm_proposal(const std::vector<R>& regions)
```

A second constructor takes a single region. This is a special case of the previous one which is provided for convenience.

```{cpp}
fmm_proposal(const R& region)
```

A copy constructor is also provided.

```{cpp}
fmm_proposal(const fmm_proposal& p)
```

### Distribution Methods {.unlisted}

The following functions make use of the proposal as a distribution. They are especially utilized in rejection sampling.

```{cpp}
std::vector<T> r(unsigned int n = 1) const;                           // <1>
std::pair<std::vector<T>, std::vector<unsigned int>>                  // <2>
	r_ext(unsigned int n = 1) const;
double d(const T& x, bool normalize = true, bool log = false) const;  // <3>
double w_major(const T& x, bool log = true) const;                    // <4>
double d_target_unnorm(const T& x, bool log = true) const;            // <5>
```

1. Draw $n$ variates of type `T` from the proposal.
2. Draw $n$ variates of type `T` from the proposal and retain the indices of the regions used in each draw. The result is an STL pair whose first element is the vector of draws and second element is the vector of indices.
3. Evaluate the density $h$ on the given $x$. If `normalize = false`, $h_0(x)$ is computed; otherwise $h(x)$ is computed.
4. Evaluate the majorized weight function $\overline{w}(x)$.
5. Evaluate the unnormalized target $f_0(x) = w(x) g(x)$.

For the `log` argument, the value on the log-scale is returned if `true`; otherwise, the value on the original scale is returned.

### Accessors {.unlisted #sec-api-proposal-accessors}

The following accessors are provided.

```{cpp}
Rcpp::NumericVector xi_upper(bool log = true) const;        // <1>
Rcpp::NumericVector xi_lower(bool log = true) const;        // <2>
Rcpp::LogicalVector bifurcatable() const;                   // <3>
Rcpp::IntegerVector mergeable(unsigned int i) const;        // <4>
Rcpp::NumericVector pi(bool log = false) const;             // <5>
Rcpp::NumericVector bound_contrib(bool log = false) const;  // <6>
double bound(bool log = false) const;                       // <7>
double nc(bool log = false) const;                          // <8>
unsigned int size() const;                                  // <9>
```

1. Get the constants $\overline{\xi}_1, \ldots, \overline{\xi}_N$.
2. Get the constants $\underline{\xi}_1, \ldots, \underline{\xi}_N$.
3. Get a vector of $N$ logical values indicating whether the corresponding regions can be bifurcated. For example, when the support `T` is `int`, a region $(a,b]$ containing one integer should not be bifurcated because one of the two resulting regions will not contain any points of the support.
4. Get a vector of zero-based indices corresponding to regions which are mergeable to the region with index `i`. This vector does not include `i` itself and may be empty.
5. Get the mixing proportions $\pi_1, \ldots, \pi_N$.
6. Get the contributions $\rho_1, \ldots, \rho_N$ to bound \eqref{eqn:bound} corresponding to the $N$ regions.
7. Get the overall rejection bound \eqref{eqn:bound}.
8. Get the normalizing constant $\psi_N$
9. Get the number of regions $N$.

Methods above with the a `log` argument return values on the log-scale when `log = true`; values are returned on the original scale otherwise.

### Iterators {.unlisted}

The following methods provide iterators to the set of regions in the proposal: const for read-only access and non-const for read and write access. Using them can be more efficient than the accessors in @sec-api-proposal because they do not make a copy of the data.

```{cpp}
std::set<R>::const_iterator begin() const;
std::set<R>::const_iterator end() const;
std::set<R>::iterator begin() const;
std::set<R>::iterator end() const;
```

### Refining the Proposal {.unlisted}

Two functions are provided to refine the proposal from $\mathscr{D}_1, \ldots, \mathscr{D}_N$ into a finer partition. Both variants return a vector $\rho_+^{(0)}, \ldots, \rho_+^{(N_0)}$ which represents values of the bound \eqref{eqn:bound} at each refinement step; $N_0$ is the number of steps taken, $\rho_+^{(j)}$ represents the value at refinement $j = 1, \ldots, N_0$, and $\rho_+^{(0)}$ represents the initial value. The first variant partitions at a given vector of knots.

```{cpp}
Rcpp::NumericVector refine(
	const std::vector<T>& knots, // <1>
	bool log = true              // <2>
)
```

1. A vector of knots.
2. If `log = true` return bound values on the log-scale.

The second variant uses the rule of thumb for sequential knot selection from @vws-2026. Refining will halt when \eqref{eqn:bound} reduces below `tol`; which is possible when `tol` is positive; otherwise, refining halts after `N` steps.

```{cpp}
Rcpp::NumericVector refine(
	unsigned int N,                        // <1>
	double tol = 0,                        // <2>
	bool greedy = false,                   // <3>
	unsigned int report = fntl::uint_max,  // <4>
	bool log = true                        // <5>
)
```

1. Maximum number of refinements to make.
2. Tolerance for \eqref{eqn:bound}; refinement will halt if $\rho_+$ reaches this.
3. If `greedy = true`, the region with the largest $\rho_\ell$ is always selected for partitioning; otherwise, regions are selected with probabilities proportional to $\rho_1, \ldots, \rho_N$.
4. If `log = true` return bound values $\rho_+^{(0)}, \ldots, \rho_+^{(N)}$ on the log-scale.
5. The period at which progress is reported to the console. E.g., use `period = 2` to report progress every two selections.

The two `refine` functions depend on methods `is_bifurcatable` and `bifurcate` of the region. These determine how to bifurcate the region and if it is eligible to be bifurcated. For a subclass `MyRegion` of `Region` whose elements are of type `T`, these methods should be defined as follows.

```{cpp}
bool is_bifurcatable() const;                             // <1>
std::pair<MyRegion,MyRegion> bifurcate() const;           // <2>
std::pair<MyRegion,MyRegion> bifurcate(const T& x) const; // <3>
```
1. Indicator of whether the current region can be bifurcated.
2. Bifurcate the current region, but without specifying a point of reference as an argument.
3. Bifurcate the current region using `x` as a point of reference.

The `seq` function is provided as a convenience to generate equally-spaced knots for univariate real-valued intervals.

```{cpp}
std::vector<double> seq(double lo, double hi, unsigned int N,
	bool endpoints = false)
```

### Coarsening the Proposal {.unlisted}

The `merge` function merges the pair of regions with indices `i` and `j` so that number of regions $N$ in the proposal is reduced by one. The return value is the updated value of bound \eqref{eqn:bound} after merging; this value is returned on the log-scale if argument `log = true`.

```{cpp}
double merge(unsigned int i, unsigned int j, bool log = true)
```

The `merge` function depends on methods `is_mergeable` and `merge` of the region. These determine how to merge two regions and whether they are eligible to merged. For a subclass `MyRegion` of `Region`, these methods should be defined as follows.

```
bool is_mergeable(const MyRegion& x) const; // <1>
MyRegion merge(const MyRegion& x) const;    // <2>
```
1. Indicator of whether the region in argument `x` can be merged with this region.
2. Construct a region by merging the region in argument `x` with this region.

### Summary Methods {.unlisted}

Several methods are provided to summarize the regions in the proposal. @tbl-summary describes the columns in the summary and rows correspond to the $N$ regions.

```{cpp}
Rcpp::DataFrame summary() const;       // <1>
void print(unsigned int n = 5) const;  // <2>
```

1. Get a data frame with the summary.
2. Print summary to the console.

| Column | Description |
|-----|-----|
`Region` | String describing region $j$.
`log_xi_upper` | $\log \overline{\xi}_j$
`log_xi_lower` | $\log \underline{\xi}_j$
`log_volume` | $\log \rho_j$

: Data frame returned by summary method of `fmm_proposal`. {#tbl-summary tbl-colwidths="\[30,70\]"}

## Region Base Class {#sec-api-region}

`Region` is an abstract base class whose interface represents the problem-specific logic that must be coded to implement VWS. Users create a subclass of this method to construct a proposal for a given problem. However, for the most common application of VWS - univariate support with a constant majorizer - users may be able to leverage the provided subclasses in Sections [-@sec-api-realconstregion] and [-@sec-api-intconstregion].

The class has one template argument `T`, which is the data type for the underlying distribution.

```{cpp}
template <class T>
class Region { ... };
```

The interface consists of the following public methods. These are abstract and must be implemented in a subclass.

```{cpp}
virtual double d_base(const T& x, bool log = false) const = 0; // <1>
virtual std::vector<T> r(unsigned int n) const = 0;            // <2>
virtual bool s(const T& x) const = 0;                          // <3>
virtual double w(const T& x, bool log = true) const = 0;       // <4>
virtual double w_major(const T& x, bool log = true) const = 0; // <5>
virtual double xi_upper(bool log = true) const = 0;            // <6>
virtual double xi_lower(bool log = true) const = 0;            // <7>
virtual double bound(bool log) const                           // <8>
virtual std::string description() const = 0;                   // <9>
```

1. Evaluate the density function $g$ of the base distribution.
2. Generate a vector of $n$ draws from $g_j$ specific to this region.
3. Indicator of whether $x$ is in the support for $g_j$ specific to this region.
4. The weight function $w$.
5. Majorized weight function $\overline{w}_j$ for this region.
6. The quantity $\overline{\xi}_j$ for this region.
7. The quantity $\underline{\xi}_j$ for this region.
8. Compute the portion $(\overline{\xi}_j - \underline{\xi}_j) / \overline{\xi}_j$ of bound \eqref{eqn:bound} specific to this region.
9. A string that describes this region; used for printing to the console.

The argument `log = true` in the methods above requests values to be returned on the log-scale.

Subclasses of `Region` are expected to provide several additional functions which are not specified in the interface. Suppose `MyRegion` is a subclass of `Region` whose elements are variables of type `type` (e.g., `double`); its definition should include the following.

```{cpp}
class MyRegion : public Region<type>
{
public:
    std::pair<MyRegion,MyRegion> bifurcate() const;              // <1>
    std::pair<MyRegion,MyRegion> bifurcate(const type& x) const; // <2>
    virtual bool is_bifurcatable() const;                        // <3>
    virtual bool is_mergeable(const region<T>& x) const;         // <4>

    MyRegion singleton(const type& x) const;                     // <5>
    bool operator<(const MyRegion& x) const;                     // <6>
    bool operator==(const MyRegion& x) const;                    // <7>
    const MyRegion& operator=(const MyRegion& x);                // <8>
...
};
```
1. Bifurcate the current region into two regions; this version takes no arguments and should include logic to decide where to make the split.
2. Bifurcate the current region into two regions; the argument `x` is used as the location to make the split.
3. Indicator of whether region can be bifurcated.
4. Indicator of whether region can be merged with region `x`. For example, if regions are intervals, this condition may be equivalent to testing whether the current regions is adjacent to `x`.
5. Return a region on the singleton set $(x,x]$ with the information in the current object.
6. Comparison operator; used to order regions of this type.
7. Equality operator; check whether two regions of this type are equal.
8. Assignment operator; assign information from `x` to this object.

## Region on Real-Valued Support with Constant Majorizer {#sec-api-realconstregion}

`real_const_region` is a subclass of `Region`, defined in @sec-api-region, specifically for univariate problems with continuous support where $\overline{w}(x) = \sum_{j=1}^N \overline{w}_j \ind\{x \in \mathscr{D}_j\}$ is constructed from constants $\overline{w}_1, \ldots, \overline{w}_N$. It assumes a constant minorizer $\underline{w}(x) = \sum_{j=1}^N \underline{w}_j \ind\{x \in \mathscr{D}_j\}$ is constructed from constants $\underline{w}_1, \ldots, \underline{w}_N$. The $\overline{w}_j$ and $\underline{w}_j$ are obtained using numerical optimization by default; however, user-supplied methods can be supplied to identify these values when it is possible.

### Constructors {.unlisted #sec-api-realconstregion-constructors}

`real_const_region` has the following constructors. The first creates a region based on the interval $(a,b]$; the second creates a region based on the singleton set $\{a\}$, which is intended primarily for internal use.

```{cpp}
real_const_region(
	double a,                                 // <1>
	double b,                                 // <2>
	const dfdb& w,                            // <3>
	const univariate_helper& helper,          // <4>
    const optimizer& maxopt = maxopt_default, // <5>
    const optimizer& minopt = minopt_default, // <6>
	const midpoint& mid = midpoint_default    // <7>
)

real_const_region(
	double a,                                 // <1>
	const dfdb& w,                            // <3>
	const univariate_helper& helper,          // <4>
    const optimizer& maxopt = maxopt_default, // <5>
    const optimizer& minopt = minopt_default, // <6>
	const midpoint& mid = midpoint_default    // <7>
)
```

1. Lower limit of interval that defines this region.
2. Upper limit of interval that defines this region.
3. Weight function $w$ for the target distribution.
4. A container with operations of the base distribution $g$.
5. A function of type `optimizer` that maximizes `w` on the given region.
6. A function of type `optimizer` that minimizes `w` on the given region.
7. A function of type `midpoint` to compute the midpoint of region's interval.

The default optimizers, `maxopt_default` and `minopt_default`, use a hybrid numerical optimization method in @sec-api-opt to optimize `w`.

The function `midpoint_default` is a version of the arithmetic midpoint with special handling for infinite endpoints. It is used primarily to bifurcate a given region into two smaller regions. In particular, the current implementation is
<!-- -->
\begin{align*}
\text{midpoint}(a,b) =
\begin{cases}
0, & \text{if $a = -\infty$ and $b = \infty$}, \\
b \cdot 2^{-\text{sign}(b)} - 1, & \text{if $a = -\infty$ and $b < \infty$}, \\
a \cdot 2^{\text{sign}(a)} + 1, & \text{if $a > -\infty$ and $b = \infty$}, \\
(a + b) / 2, & \text{otherwise},
\end{cases}
\end{align*}
<!-- -->
where $\text{sign}(x) = \ind(x > 0) - \ind(x < 0)$. If both endpoints are infinite, zero is returned. When one endpoint is infinite, the finite endpoint is reduced in magnitude by a factor of $1/2$ and brought in by an additional unit of one. The exponential factor helps to avoid large numbers of wasted regions during knot selection. The additive factor accommodates endpoints that are near zero; without it, $\text{midpoint}(a,b)$ only brings the finite endpoint closer to zero which likely to create wasteful bifurcations.

### Methods {.unlisted}

The `midpoint` method returns a midpoint for the current region.

```{cpp}
double real_const_region::midpoint() const
```

The `bifurcate` method returns two disjoint regions whose union is the current region. The result is given as an STL pair. The first version bifurcates at the midpoint of the current region, determined by the `midpoint` method. The second version partitions at the given $x$.

```{cpp}
std::pair<real_const_region,real_const_region> bifurcate() const;
std::pair<real_const_region,real_const_region> bifurcate(const double& x) const;
```

The `is_bifurcatable` function always returns `true` in the case of a continuous support.

```{cpp}
bool is_bifurcatable() const;
```

The `merge` method produces a region consisting of the current region and the argument `x` merged together.

```{cpp}
real_const_region merge(const real_const_region& x) const;
```

The `is_mergeable` function is an indicator of whether the current region can be merged with the argument `x`. This equates to `true` if and only if the two regions are adjacent intervals whose endpoints meet at a common value (on either side).

```{cpp}
bool is_mergeable(const real_const_region& x) const;
```

Several functions are provided to update an existing region.

```{cpp}
void set_w(const dfdb& w);                         // <1>
void set_helper(const univariate_helper& helper);  // <2>
void set_maxopt(const optimizer& maxopt);          // <3>
void set_minopt(const optimizer& minopt);          // <4>
void set_mid(const vws::midpoint& mid);            // <5>
void init();                                       // <6>
```
1. Update the weight function.
2. Update the univariate helper which represents the base distribution.
3. Update the maximization routine for the majorizer.
4. Update the minimization routine for the minorizer.
5. Update the function to find the midpoint.
6. Recompute saved quantities such as the majorization and minorization constants needed for rejection sampling.

::: {.callout-important title="Important"}
The `init` method should be applied after using the set methods above and before the region is used in any other capacity. More than one set method may be used before calling `init`.
:::

The `singleton` method returns a singleton interval $(x, x]$, using the current object's member data.

```{cpp}
real_const_region singleton(const double& x) const
```

The following methods determine an ordering of the current region and an another region specified as argument `x`. Region $(a_1, b_1]$ is considered "less than" $(a_2, b_2]$ if $b_1 \leq a_2$. The regions are considered equal if $a_1 = a_2$ and $b_1 = b_2$. Note that other elements such as $w$ are $g$ are not explicitly checked, and are assumed to be the same.

```{cpp}
bool operator<(const real_const_region& x) const
bool operator==(const real_const_region& x) const
```

The following method assigns the current region to be equal to the argument `x`.

```{cpp}
const real_const_region& operator=(const real_const_region& x)
```

## Region on Integer-Valued Support with Constant Majorizer {#sec-api-intconstregion}

The class `int_const_region` is a subclass of `real_const_region` that is specialized to integer-valued distributions.

```{cpp}
class int_const_region : public real_const_region { ... };
```

### Constructors {.unlisted}

`int_const_region` has the following constructors.

```{cpp}
int_const_region(
	double a,                                 // <1>
	double b,                                 // <2>
	const dfdb& w,                            // <3>
	const univariate_helper& helper,          // <4>
    const optimizer& maxopt = maxopt_default, // <5>
    const optimizer& minopt = minopt_default, // <6>
	const midpoint& mid = midpoint_default    // <7>
)

int_const_region(
	double a,                                 // <1>
	const dfdb& w,                            // <3>
	const univariate_helper& helper,          // <4>
    const optimizer& maxopt = maxopt_default, // <5>
    const optimizer& minopt = minopt_default, // <6>
	const midpoint& mid = midpoint_default    // <7>
)
```

1. Lower limit of interval that defines this region.
2. Upper limit of interval that defines this region.
3. Weight function $w$ for the target distribution.
4. A container with operations of the base distribution $g$.
5. A function of type `optimizer` that maximizes `w` on the given region.
6. A function of type `optimizer` that minimizes `w` on the given region.
7. A function of type `midpoint` to compute the midpoint of region's interval.

The arguments here are analogous to those in @sec-api-realconstregion-constructors. Note that `a` and `b` need not be integers here.

### Methods {.unlisted}

The following methods are specialized from `real_const_region`.

```{cpp}
std::pair<int_const_region,int_const_region> bifurcate() const;
std::pair<int_const_region,int_const_region> bifurcate(const double& x) const;
```

The `is_bifurcatable` function returns `false` for regions whose width is smaller than $1$.

```{cpp}
bool is_bifurcatable() const;
```

The `singleton` method returns a singleton interval $(x, x]$, using the current object's weight function, base distribution, etc.

```{cpp}
int_const_region singleton(const double& x) const;
```

The following methods determine an ordering of the current region and an another region specified as argument `x`. Region $(a_1, b_1]$ is considered "less than" $(a_2, b_2]$ if $b_1 \leq a_2$. The regions are considered equal if $a_1 = a_2$ and $b_1 = b_2$. Note that other elements such as $w$ are $g$ are not explicitly checked, and are assumed to be the same.

```{cpp}

bool operator<(const int_const_region& x) const
bool operator==(const int_const_region& x) const
```

The following method assigns the current region to be equal to the argument `x`.

```{cpp}
const int_const_region& operator=(const int_const_region& x)
```

## Univariate Helper {#sec-api-helper}

The class `univariate_helper` is intended for use with `real_const_region` and `int_const_region`. It encapsulates several operations needed from the base distribution $g$. These operations are specified as lambdas in the constructor.

```{cpp}
univariate_helper(
	const fntl::density& d,  // <1>
	const fntl::cdf& p,      // <2>
	const fntl::quantile& q  // <3>
)
```

1. A function to evaluate density $g$.
2. A function to evaluate the cumulative distribution function (CDF) $G$.
3. A function to evaluate the quantile function $G^{-}$.

The following methods on `univariate_helper` utilize the lambdas specified above.

```{cpp}
double d(double x, bool log = false) const;                     // <1>
double p(double q, bool lower = true, bool log = false) const;  // <2>
double q(double p, bool lower = true, bool log = false) const;  // <3>
```

1. Evaluate the density function at argument $x$. Result is on the log-scale if `log = true`.
2. Evaluate the CDF at argument $q$. Result is on the log-scale if `log = true`. Result represents $P(X \leq q)$ if `lower = true` and $P(X > q)$ otherwise.
3. Evaluate the quantile function at argument $p$. Assume $p$ is specified on the log-scale if `log = true`. Request $p$ quantile if `lower = true` and $1-p$ quantile otherwise.

## Rejection Sampling {#sec-api-rejection}

### Without Tuning {#sec-api-rejection-basic}

The following functions perform rejection sampling using a VWS proposal described in @sec-vws.

```{cpp}
template <typename T, typename R>
rejection_result<T> rejection(
	const fmm_proposal<T,R>& h,        // <1>
	unsigned int n,                    // <2>
	const rejection_args& args         // <3>
)

template <typename T, typename R>
rejection_result<T> rejection(
	const fmm_proposal<T,R>& h,        // <1>
	unsigned int n = 1                 // <2>
)
```
1. An `fmm_proposal` object to use as the proposal.
2. The number of desired draws.
3. Additional arguments for rejection sampling. Default values are assumed in the second form.

Template arguments `T` and `R` correspond to the given proposal `h` and are described in @sec-api-proposal.

### With Tuning {#sec-api-rejection-tune}

The following functions perform rejection sampling with tuning as described in @sec-self-tune.

```{cpp}
template <typename T, typename R>
rejection_result<T> rejection_tune(
	fmm_proposal<T,R>& h,       // <1>
	unsigned int n,             // <2>
	const rejection_args& args  // <3>
)

template <typename T, typename R>
rejection_result<T> rejection_tune(
	const fmm_proposal<T,R>& h, // <1>
	unsigned int n = 1          // <2>
)
```
1. An `fmm_proposal` object to use as the proposal.
2. The number of desired draws.
3. Additional arguments for rejection sampling. Default values are assumed in the second form.

Template arguments `T` and `R` correspond to the given proposal `h` and are described in @sec-api-proposal.

### Additional Arguments {#sec-api-rejection-args}

Additional arguments to `rejection` and `rejection_tune` are provided via the following struct.

```{cpp}
struct rejection_args {
	unsigned int max_rejects = std::numeric_limits<unsigned int>::max(); // <1>
	unsigned int report = std::numeric_limits<unsigned int>::max();      // <2>
	double ratio_ub = std::exp(1e-5);                                    // <3>
	fntl::error_action action = fntl::error_action::STOP;                // <4>
	double tol_suff = 0.85;                                              // <5>
	double tol_merge = 0.01;                                             // <6>
	bool metrics = true;                                                 // <7>

	rejection_args() { };                                                // <8>
	rejection_args(SEXP obj);                                            // <9>
	operator SEXP() const;                                               // <10>
};
```

1. Maximum number of rejections to tolerate overall (among all $n$ attempted draws) before bailing out.
2. Determines period at which progress is logged to the console.
3. Maximum allowed value for the ratio $f_0(x) / h_0(x)$, which may be slightly larger than 1 due to numerical precision. An error is thrown if the ratio is larger than this value.
4. The action to take if `max_rejects` rejections is exceeded; see @tbl-action. The definition of the enum `fntl::error_action` is given in Section 2 of @fntl.
5. Tolerance $\epsilon_1$ in Algorithm \ref{alg:self-tune} used to determine whether rejection rate is sufficiently below 1. This is used only in `rejection_tune`.
6. Tolerance $\epsilon_2$ in Algorithm \ref{alg:self-tune} used to determine whether a region's contribution is low enough to be merged with another region.
This is used only in `rejection_tune`.
7. Determines whether to save additional metrics. If `true`, non-essential metrics in `rejection_result` will be populated. Otherwise, only `draws` will be populated with the generated variates.
8. Constructor that takes no arguments and initializes to default values.
9. Convert an `Rcpp::List` to a `rejection_args` struct.
10. Return a `Rcpp::List` from a `rejection_args` struct.

| Action | Effect |
|-----|-----|
`fntl::error_action::NONE`    | Error condition is ignored.
`fntl::error_action::MESSAGE` | A message is emitted without halting.
`fntl::error_action::WARNING` | A warning is emitted without halting.
`fntl::error_action::STOP`    | Halts when an error condition is encountered.

: Value of `action` used in `rejection_args.action` and its effect in `rejection`. {#tbl-action tbl-colwidths="\[35,65\]"}

### Results {#sec-api-rejection-results}

The return value of `rejection` and `rejection_tune` is a struct of the following type.

```{cpp}
template <typename T>
struct rejection_result
{
	std::vector<T> draws;               // <1>
	std::vector<unsigned int> rejects;  // <2>
	std::vector<unsigned int> tunes;    // <3>
	std::vector<unsigned int> regions;  // <4>
	std::vector<double> log_bounds;     // <5>

	operator SEXP() const;
};
```

1. Vector of draws.
2. Vector of rejection counts; the $i$th element represents number of rejections observed before accepting the $i$th draw.
3. Vector of tuning counts; the $i$th element represents number of tuning adjustments made before accepting the $i$th draw.
4. Vector of region counts; the $i$th element represents number of regions in the proposal when accepting the $i$th draw.
5. Return a `Rcpp::List` from a `rejection_args` struct.

The `tunes`, `regions`, and `log_bounds` fields are only relevant to `rejection_tune`; they are left empty in `rejection`. Furthermore, all fields aside from `draws` are left empty when `metrics = false` in `rejection_args`.

If the maximum number of rejections are reached and `rejection_args` specifies `action = fntl::error_action::STOP`, the length of the vectors `draws` and `rejects` will be less than $n$.

## Optimization on an Interval {#sec-api-opt}

The function `optimize_hybrid` is a hybrid optimization method for univariate functions $f(x) : [a,b] \rightarrow \mathbb{R}$ with bounds $x \in [a,b]$ whose endpoints may be finite or infinite. Uses Brent's method if both bounds are finite and BFGS otherwise. In the latter case, the bounds are enforced via a transformation.

```{cpp}
optimize_hybrid_result optimize_hybrid(
	const fntl::dfd& f,       // <1>
    double init,              // <2>
	double lower,             // <3>
	double upper,             // <4>
	bool maximize,            // <5>
	unsigned maxiter = 100000 // <6>
)
```

1. Objective function.
2. Initial value used with BFGS.
3. Lower bound $a$ of domain $[a,b]$.
4. Upper bound $b$ of domain $[a,b]$.
5. Logical; if `true`, optimization will be a maximization. Otherwise it is a minimization.
6. Maximum number of iterations.

The result is a `optimize_hybrid_result` struct which is defined as follows.

```{cpp}
struct optimize_hybrid_result {
	double par;            // <1>
	double value;          // <2>
	std::string method;    // <3>
	int status;            // <4>

	operator SEXP() const; // <5>
};
```

1. Final value of the optimization variable $x$.
2. Final value of the objective function $f(x)$.
3. Description of the method used to obtain the result; see @tbl-method for possible values.
4. Corresponds to a code from BFGS if it is used as `method`; otherwise $0$.
5. Return an `Rcpp::List` from a `optimize_hybrid_result` struct.

| Message | Description |
|-----|-----|
| `"Brent"` | Brent optimization method was used. |
| `"BFGS"` | BFGS method was used. |
| `"Lower Limit Inf"` | For maximization, lower limit taken as `par` because its value is `inf`. |
| `"Upper Limit Inf"` | For maximization, upper limit taken as `par` because its value is `inf`. |
| `"Lower Limit NegInf"` | For minimization, lower limit taken as `par` because its value is `-inf`. |
| `"Upper Limit NegInf"` | For minimization, upper limit taken as `par` because its value is `-inf`. |
| `"Max at Lower Limit"` | Numerical maximization used, but larger value found at lower limit. |
| `"Max at Upper Limit"` | Numerical maximization used, but larger value found at upper limit. |
| `"Min at Lower Limit"` | Numerical minimization used, but smaller value found at lower limit. |
| `"Min at Upper Limit"` | Numerical minimization used, but smaller value found at upper limit. |

: Possible values for `method` field of `optimize_hybrid_result` struct. {#tbl-method tbl-colwidths="\[28,72\]"}

## Log-Scale Arithmetic {#sec-api-logscale}

As mentioned in @sec-logscale, calculations in the package are carried out on the log-scale, where possible, to avoid issues from floating point numbers with very small or very large magnitudes. Users may want to follow this convention when implementing their own sampling problems. Several included functions help to avoid explicit exponentiation.

The following computes the scalar $f(x) = \log\{ \sum_{i=1}^n \exp(x_i) \}$ from a vector $x \in \mathbb{R}^n$ using the method described in [this StackExchange post](https://stats.stackexchange.com/a/381937).

```{cpp}
double log_sum_exp(const Rcpp::NumericVector& x)
```

The following functions carry out addition on the log scale using the property
<!-- -->
\begin{displaymath}
\log(e^x + e^y) = t + \log\{1 + \exp[s - t]\},
\quad s = \min(x,y), \quad t = \max(x,y).
\end{displaymath}
<!-- -->
The first form takes scalar arguments $x$ and $y$. The second and third forms take $x, y \in \mathbb{R}^n$ and produce a vector with elements $\log(e^{x_i} + e^{y_i})$, $i = 1, \ldots, n$.

```{cpp}
double log_add2_exp(double x, double y)

std::vector<double> log_add2_exp(
	const std::vector<double>& x,
	const std::vector<double>& y
)

Rcpp::NumericVector log_add2_exp(
	const Rcpp::NumericVector& x,
	const Rcpp::NumericVector& y
)
```

The following carry out subtraction on the log scale using the property
<!-- -->
\begin{displaymath}
\log(e^x - e^y) = x + \log\{1 - \exp[y - x]\}.
\end{displaymath}
<!-- -->
When $x$ is smaller than $y$, the results is assumed to be `nan`. The first form takes scalar arguments $x$ and $y$. The second and third forms take $x, y \in \mathbb{R}^n$ and produce a vector with elements $\log(e^{x_i} - e^{y_i})$, $i = 1, \ldots, n$. 

```{cpp}
double log_sub2_exp(double x, double y)

std::vector<double> log_sub2_exp(
	const std::vector<double>& x,
	const std::vector<double>& y
)

Rcpp::NumericVector log_sub2_exp(
	const Rcpp::NumericVector& x,
	const Rcpp::NumericVector& y
)
```

## Generating from a Discrete Distribution {#sec-api-discrete}

We make use of the Gumbel trick [e.g., @Huijben-EtAl-2023] to draw from a discrete distribution with probabilities $p_1, \ldots, p_k$. This approach allows probabilities to be specified on the log-scale without the need to exponentiate or normalize them so that they sum to $1$. The Gumbel trick generates a draw $x$ from the desired discrete distribution via
<!-- -->
$$
X = \argmax \{ Z_1 + \log p_1, \ldots, Z_k + \log p_k \},
\quad Z_1, \ldots, Z_k \iid \text{Gumbel}(0,1),
$$
<!-- -->
where $\text{Gumbel}(0,1)$ is a standard Gumbel distribution with density $f(x) = e^{-(x + e^{x})} \ind(x > 0)$.

### Categorical Distribution {.unlisted}

The following functions generate a draw of $X$. The first form generates a single variate and the second form generates a sample of size $n$.

```{cpp}
unsigned int r_categ(
	const Rcpp::NumericVector& p,  // <2>
	bool log = false,              // <3>
	bool one_based = false         // <4>
)

Rcpp::IntegerVector r_categ(
	unsigned int n,                // <1>
	const Rcpp::NumericVector& p,  // <2>
	bool log = false,              // <3>
	bool one_based = false         // <4>
)
```

1. Desired sample size.
2. Vector of probabilities $p_1, \ldots, p_k$.
3. Logical; if `true`, argument `p` is interpreted as $\log p_1, \ldots, \log p_k$. Otherwise, it is interpreted as $p_1, \ldots, p_k$ on the original scale.
4. Logical; if `true`, support is assumed to be $\{ 1, \ldots, k \}$, where $k$ is length of the given `p`. Otherwise it is assumed to be $\{ 0, \ldots, k-1 \}$. The former is useful to generate indices in C++ while the latter is useful for indices in R.

### Gumbel Distribution {.unlisted}

The following functions are provided for the Gumbel distribution with location parameter $\mu$ and scale $\sigma$. They compute the density, CDF, quantile, and variate generation, respectively. The first group operates on scalars and the second are vectorized versions which operate on an independent and identically distributed sample.

```{cpp}
double d_gumbel(
	double x,                      // <1>
	double mu = 0,                 // <5>
	double sigma = 1,              // <6>
	bool log = false               // <7>
)

double p_gumbel(
	double q,                      // <2>
	double mu = 0,                 // <5>
	double sigma = 1,              // <6>
	bool lower = true,             // <8>
	bool log = false               // <7>
)

double q_gumbel(
	double p,                      // <3>
	double mu = 0,                 // <5>
	double sigma = 1,              // <6>
	bool lower = true,             // <8>
	bool log = false               // <7>
)

double r_gumbel(
	double mu = 0,                 // <5>
	double sigma = 1               // <7>
)

Rcpp::NumericVector d_gumbel(
	const Rcpp::NumericVector& x,  // <1>
	double mu = 0,                 // <5>
	double sigma = 1,              // <6>
	bool log = false               // <7>
)

Rcpp::NumericVector p_gumbel(
	const Rcpp::NumericVector& q,  // <2>
	double mu = 0,                 // <5>
	double sigma = 1,              // <6>
	bool lower = true,             // <8>
	bool log = false               // <7>
)

Rcpp::NumericVector q_gumbel(
	const Rcpp::NumericVector& p,  // <3>
	double mu = 0,                 // <5>
	double sigma = 1,              // <6>
	bool lower = true,             // <8>
	bool log = false               // <7>
)

Rcpp::NumericVector r_gumbel(
	unsigned int n,                // <4>
	double mu = 0,                 // <5>
	double sigma = 1               // <6>
)
```
1. Argument to the density.
2. Argument to the CDF.
3. Argument to the quantile function.
4. Number of draws.
5. Location parameter.
6. Scale parameter.
7. Logical; if `true`, probability arguments and return values are on the log-scale
8. Logical; if `true`, probability arguments and return values are lower tailed in the form $P[X \leq x]$; otherwise, $P[X > x]$.

# Example: Von Mises Fisher Distribution {#sec-vmf}

Let us consider generating from the von Mises Fisher (VMF) distribution as in @vws-2026. VMF is useful in modeling directional data whose support is the sphere $\mathbb{S}^{d-1} = \{ \vec{v} \in \mathbb{R}^d : \vec{v}^\top \vec{v} = 1 \}$ [e.g., @Mardia-Jupp-1999]. A random variable $\vec{V}$ with distribution $\text{VMF}_d(\vec{\mu}, \kappa)$ has density
<!-- -->
\begin{align*}
f_{\text{VMF}}(\vec{v}) = \frac{
\kappa^{d / 2 - 1}
}{
(2 \pi)^{d/2} I_{d/2 - 1}(\kappa)
} \exp(\kappa \cdot \vec{\mu}^\top \vec{v}) \cdot \ind\{\vec{v} \in \mathbb{S}^{d-1}\},
\end{align*}
<!-- -->
with modified Bessel function of the first kind
<!-- -->
\begin{math}
I_{\nu}(x) = \sum_{m=0}^\infty \{m! \cdot \Gamma(m + \nu + 1)\}^{-1} (\frac{x}{2})^{2m + \nu}
\end{math}
<!-- -->
and gamma function $\Gamma(x) = \int_0^\infty t^{x-1} e^{-t} dt$. Parameters $\vec{\mu} \in \mathbb{S}^{d-1}$ and $\kappa > 0$ determine the orientation on the sphere and the concentration, respectively. First consider $\vec{\mu}_0 = (1, 0, \ldots, 0)$. A random variable $\vec{V}_0 \sim \text{VMF}_d(\vec{\mu}_0, \kappa)$ can be obtained from the transformation
<!-- -->
\begin{align}
\vec{V}_0 = \left( X, \sqrt{1 - X^2} \cdot \vec{U} \right),
\label{eqn:vmf-tx}
\end{align}
<!-- -->
where $\vec{U}$ is a uniform random variable on the sphere $\mathbb{S}^{d-2}$ and $X$ has density
<!-- -->
\begin{align}
f(x) =
\frac{
(\kappa / 2)^{d/2 - 1} (1 - x^2)^{(d-3)/2} \exp(\kappa x)
}{
\sqrt{\pi} \cdot I_{d/2 - 1}(\kappa) \cdot \Gamma((d-1)/2)
} \cdot \ind(-1 < x < 1).
\label{eqn:vmf-target}
\end{align}
<!-- -->
To obtain a draw from $\vec{V} \sim \text{VMF}_d(\vec{\mu}, \kappa)$ with an arbitrary $\vec{\mu}$, we can rotate $\vec{V} = \vec{Q} \vec{V}_0$ using an orthonormal matrix $\vec{Q}$ whose first column is $\vec{\mu}$. To generate from $\vec{V}_0$ in \eqref{eqn:vmf-tx}, we may draw $\vec{U} = \vec{Z} / \sqrt{\vec{Z}^\top \vec{Z}}$ from $\vec{Z} \sim \text{N}(\vec{0}, \vec{I}_{d-1})$ and $X$ independently from \eqref{eqn:vmf-target}. In the following, we consider the use of VWS to draw the univariate random variable $X$. Before proceeding, we identify a useful distribution.

::: {#def-texp}
Denote $X \sim \text{Exp}_{(a,b)}(\kappa)$ as a "doubly truncated" Exponential random variable with density
<!-- -->
\begin{align*}
g(x) = \frac{\kappa e^{\kappa x}}{e^{\kappa b} - e^{\kappa a}} \cdot \ind(a < x < b),
\end{align*}
<!-- -->
where $-\infty < a < b < \infty$ and rate $\kappa$ may be any real number. The CDF and quantile function corresponding to $g$ are
<!-- -->
\begin{align*}
&G(x) = \frac{e^{\kappa x} - e^{\kappa a}}{e^{\kappa b} - e^{\kappa a}} \cdot,
\quad x \in (a, b), \\
%
&G^{-1}(\varphi) = \frac{1}{\kappa} \log\left[e^{\kappa a} + \varphi (e^{\kappa b} - e^{\kappa a}) \right], \quad \varphi \in [0,1].
\end{align*}
\remarkqed
:::

Returning to target \eqref{eqn:vmf-target}, let us decompose $f$ into
<!-- -->
\begin{align*}
f(x) \propto
\underbrace{(1 - x^2)^{(d-3)/2}}_{w(x)}
\underbrace{\exp(\kappa x) \cdot \ind(-1 < x < 1)}_{g_0(x)},
\end{align*}
<!-- -->
excluding terms from the normalizing constant, so that $w$ is the weight function and $g_0$ is proportional to the density of $\text{Exp}_{(-1,1)}(\kappa)$,
<!-- -->
\begin{align*}
g(x) = \frac{\kappa e^{\kappa x}}{e^\kappa - e^{-\kappa}} \cdot \ind(-1 < x < 1).
\end{align*}
<!-- -->
@sec-vmf-const-default obtains a VWS sampler with this decomposition using a constant majorizer. @sec-vmf-const-custom replaces the default numerical optimization with custom code, which reduces the amount of computational overhead. @sec-vmf-linear considers a linear majorizer which is substantially more involved but also obtains substantially lower rejection rates with a moderate number of regions. Codes for this example are in the folder `examples/vmf`. C++ functions for the $\text{Exp}_{(a,b)}(\kappa)$ distribution from @def-texp are given in the file `examples/vmf/texp.h`.

::: {#rem-vmf-caveat}
A caveat of this decomposition is that, in the $d < 3$ case, $w(x) \rightarrow \infty$ as $x$ approaches $\pm 1$. One way to avoid this is by truncating the support to $(\alpha_0, \alpha_N] = (-1 + \epsilon, 1 - \epsilon]$ for a small $\epsilon > 0$. Rejection sampling can proceed using the truncated support if the exclusion of $(-1, -1 + \epsilon] \cup (1 - \epsilon, 1]$ is known to have a negligible impact on the result. Otherwise, @vws-2026 mention another strategy where the support is initially truncated and gradually expanded as rejections are encountered. In this document, we assume $(\alpha_0, \alpha_N] = (-1, 1]$ for $d > 3$ and a fixed truncation $(\alpha_0, \alpha_N] = (-1 + \epsilon, 1 - \epsilon]$ for $d < 3$.
:::

```{r}
#| echo: false
#| prompt: true
Rcpp::sourceCpp("examples/vmf/functions.cpp")
```

## Constant Majorizer with Numerical Optimization {#sec-vmf-const-default}

We now give our first example demonstrating the `vws` package. We consider a constant majorizer for $w$ which uses the default numerical optimization routine to identify appropriate constants $\overline{w}_j$ and $\underline{w}_j$. Because this is our first example, the source file for the sampler (`examples/vmf/vmf-v1.cpp`) is displayed in its entirety as follows.

``` {.cpp include="examples/vmf/vmf-v1.cpp" code-line-numbers="true"}
```

1. Ensure that this code links with the `vws` and `fntl` packages during compilation.
2. Include the header for C++ framework in `vws`.
3. Include `texp.h`, which provides functions described in @def-texp.
4. Define a C++ function which invokes the sampler and export it for use in R.
5. Prepare a struct with extra arguments for rejection sampling.
6. Define the weight function as a lambda. We are careful to avoid `nan` values that can occur when $x = \pm 1$. Computations are carried out on the log-scale to avoid numerical loss of precision.
7. Specify the density, CDF, and quantile function of the base distribution $\text{Exp}_{(-1,1)}(\kappa)$. The density, CDF, and quantile function types are defined in @sec-api-typedefs.
8. Create a "helper" object as a container for the distribution functions.
9. Construct a `real_const_region` which contains all problem-specific logic of the sampler. We construct one initial region which contains the entire support $(-1,1]$.
10. Construct an `fmm_proposal` based on our initial region `supp`. The first template argument specifies that the data type of the support is `double`; the second specifies that regions (which include the logic of the sampler) are of type `real_const_region`.
11. Request the proposal `h` to refine itself `N-1` times using Algorithm \ref{alg:refine} so that there are `N` regions.
12. Carry out rejection sampling with proposal `h`.
13. Assemble an `Rcpp::List` to return to the caller. It contains the draws in element `draws`, a vector of rejection counts in `rejects` where the $i$th element represents the number of rejections for the $i$th draw, and a vector in element `lbdd` with the $N$ bounds $\rho_+^{(1)}, \ldots, \rho_+^{(N)}$, where $\rho_+^{(j)}$ is the value achieved with $j$ regions.

The name of the function `r_vmf_pre_v1` reflects that this is our first version of the sampler for target \eqref{eqn:vmf-target}, which is a precursor to transformation \eqref{eqn:vmf-tx} to obtain a VMF random variable. The following R snippet builds the C++ code and invokes the sampling function.

```{r}
#| prompt: true
Rcpp::sourceCpp("examples/vmf/vmf-v1.cpp")
out1 = r_vmf_pre_v1(n = 1000, kappa = 5, d = 4, N = 50, tol = 0.10)
head(out1$draws)
```

@fig-vmf-const-refine plots the bound for the rejection probability during the refinement process, which is captured in the variable `out$lbdd`. @fig-vmf-const-draws plots the empirical distribution of the draws and compares them to the density.

```{r}
#| out-width: 60%
#| fig-cap: |
#|   Refinement for VMF example with constant majorizer.
#| label: fig-vmf-const-refine
#| echo: false

df = data.frame(bdd = exp(out1$lbdd)) %>%
	mutate(step = row_number() - 1)
ggplot(df) +
	geom_line(aes(x = step, y = bdd)) +
	scale_x_continuous(breaks = df$step[df$step %% 2 == 0]) +
	scale_y_continuous(n.breaks = 10) +
	xlab("Step") +
	ylab("Bound") +
	theme_minimal()
```

```{r}
#| out-width: 60%
#| fig-cap: |
#|   Empirical distribution of draws (solid) versus target (dashed) for VMF
#|   example with constant majorizer.
#| label: fig-vmf-const-draws
#| echo: false

data.frame(x = out1$draws) %>%
	ggplot() +
	geom_density(aes(x = x), col = "black") +
	geom_function(fun = d_target, args = list(kappa = 5, d = 4), lty = 2) +
	ylab("Empirical Density") +
	theme_minimal()
```

## Constant Majorizer with Custom Optimization {#sec-vmf-const-custom}

It is not difficult to find the minimum and maximum of the function $w$ on a given interval for this example. We can reduce some of the computational burden by providing functions to compute these two values.

We have
<!-- -->
\begin{align*}
&\log w(x) = \frac{d-3}{2} \log(1 - x^2), \\
&\frac{d}{dx} \log w(x) = -(d-3) \frac{x}{1 - x^2}, \\
&\frac{d^2}{dx^2} \log w(x) = -(d-3) \frac{1 + x^2}{(1 - x^2)^2}.
\end{align*}
<!-- -->
Therefore, it is seen that $\log w(x)$ is concave when $d > 3$, convex when $d = 2$, and constant otherwise. When $d > 3$, $\frac{d}{dx} \log w(x)$ is positive for $x \in (-1, 0)$, negative for $x \in (0, 1)$, and has root $x = 0$; therefore, $\log w(x)$ is unimodal on $(-1, 1)$ with a maximum at $x = 0$. When $d = 2$, the point $x = 0$ is instead a minimum of $\log w(x)$. Finally, $w$ is a constant in the case $d = 3$. The following code implements these maximization and minimization routines as lambdas.

```{cpp}
vws::optimizer opt1 = [&](const vws::dfdb& w, double lo, double hi, bool log)
{
	double x = 0;
	if (hi < 0) {
		x = hi;
	} else if (lo > 0){
		x = lo;
	}
	double out = w(x, true);
	return log ? out : std::exp(out);
};
```

```{cpp}
vws::optimizer opt2 = [&](const vws::dfdb& w, double lo, double hi, bool log)
{
	double w_lo = w(lo, true);
	double w_hi = w(hi, true);
	double out = std::min(w_lo, w_hi);
	return log ? out : std::exp(out);
};
```

The following snippet creates pointers `maxopt` and `minopt` for the maximizer and minimizer functions. The condition $d > 3$ is checked to determine which of `opt1` and `opt2` is the maximizer and which is the minimizer.

```{cpp}
vws::optimizer* maxopt;
vws::optimizer* minopt;
if (d >= 3) {
	maxopt = &opt1;
	minopt = &opt2;
} else {
	minopt = &opt1;
	maxopt = &opt2;
}
```

Finally, we create the initial region `supp`. The `maxopt` and `minopt` functions are provided to the constructor as additional arguments. Note that we dereference our pointers and pass the function objects themselves (which are taken as references by the constructor).

```{cpp}
vws::univariate_helper helper(df, pf, qf)
vws::real_const_region supp(-1, 1, w, helper, *maxopt, *minopt)
```

The following R snippet builds the C++ code and invokes the sampling function.

```{r}
#| prompt: true
Rcpp::sourceCpp("examples/vmf/vmf-v2.cpp")
out2 = r_vmf_pre_v2(n = 1000, kappa = 5, d = 4, N = 50, tol = 0.10)
head(out2$draws)
```

## Linear Majorizer {#sec-vmf-linear}

We noted in @sec-vmf-const-custom that $w$ is log-convex when $d < 3$, log-concave when $d > 3$, and a constant otherwise. Therefore, we can majorize $w$ with exponentiated linear functions of the form $\overline{w}_j(x) = \exp\{ \overline{\beta}_{0j} + \overline{\beta}_{1j} x \}$. This yields the expression
<!-- -->
\begin{align*}
\overline{\xi}_j
= \int_{\alpha_{j-1}}^{\alpha_j} \overline{w}_j(x) g(x) dx
= \frac{\kappa \exp\{ \overline{\beta}_{0j} \}}{(\kappa + \overline{\beta}_{1j}) (e^{\kappa \alpha_N} - e^{\kappa \alpha_0})} \left\{
\exp\{(\kappa + \overline{\beta}_{1j}) \alpha_j\} - \exp\{(\kappa + \overline{\beta}_{1j}) \alpha_{j-1} \}
\right\}.
\end{align*}
<!-- -->
The proposal $h$ is then a finite mixture $h(x) = \sum_{j=1}^N \pi_j g_j(x)$ with
<!-- -->
\begin{align*}
g_j(x)
&= \overline{w}_j(x) g(x) \ind\{ x \in \mathscr{D}_j \} / \overline{\xi}_j \\
&= \frac{
(\kappa + \overline{\beta}_{1j}) \exp\{(\kappa + \overline{\beta}_{1j}) x\}
}{
\exp\{(\kappa + \overline{\beta}_{1j}) \alpha_j\} - \exp\{(\kappa + \overline{\beta}_{1j}) \alpha_{j-1} \}
}
\cdot \ind(\alpha_{j-1} < x \leq \alpha_j),
\end{align*}
<!-- -->
the density of $\text{Exp}_{(\alpha_{j-1}, \alpha_j]}(\kappa + \overline{\beta}_{1j})$. @rem-expansion-point is used to select the expansion point $c$ to determine coefficients for the majorizer in the log-concave case, with
<!-- -->
\begin{align*}
M_j(s) = \int_{\alpha_{j-1}}^{\alpha_j} e^{sx} g(x) dx
= \frac{e^{s\alpha_j} - e^{s\alpha_{j-1}}}{s(\alpha_j - \alpha_{j-1})}.
\end{align*}
<!-- -->
To compute \eqref{eqn:bound}, we assume the "trivial" minorizer $\underline{w}_j(x) = w(x)$ so that
<!-- -->
\begin{align*}
\underline{\xi}_j
= \int_{\alpha_{j-1}}^{\alpha_j} w(x) g(x) dx
= \int_{\alpha_{j-1}}^{\alpha_j} \frac{(1-x^2)^{(d-3)/2} \kappa e^{\kappa x}}{e^{\kappa \alpha_N} - e^{\kappa \alpha_0}} dx,
\end{align*}
<!-- -->
which we compute using numerical integration. The proposal is implemented with the `vws` package by inheriting from the abstract `Region` base class (@sec-api-region) and implementing each of the functions using the expressions above. We name the resulting subclass `LinearVWSRegion`; its complete code is given in `examples/vmf/LinearVWSRegion.h`. The code in `examples/vmf/vmf-v3.cpp` instantiates a proposal with regions of type `LinearVWSRegion`, invokes rejection sampling, and returns an `Rcpp::List` with the results. This code is displayed below.

``` {.cpp include="examples/vmf/vmf-v3.cpp" code-line-numbers="true"}
```

The following R snippet builds the C++ code and invokes the sampling function.

```{r}
#| prompt: true
Rcpp::sourceCpp("examples/vmf/vmf-v3.cpp")
out3 = r_vmf_pre_v3(n = 1000, kappa = 5, d = 4, N = 50, tol = 0.01)
head(out3$draws)
```

::: {#rem-vmf-endpoints}
This code will fail if we attempt to use it with $d = 2$ as described in @rem-vmf-endpoints. Here, $\log w(x) \rightarrow \infty$ as $x \rightarrow \pm 1$ so that the first and last regions cannot be bounded by an exponentiated linear form. This can be averted by truncating the support to $(-1 + \epsilon, 1 - \epsilon]$ for a small $\epsilon > 0$.
\remarkqed
:::

@fig-vmf-const-draws plots the empirical distribution of the draws and compares them to the density. @fig-vmf-refine plots the bound for the rejection probability for this sampler after each step of refining, along with those from the previous two versions. A substantial improvement in efficiency is seen here; fewer regions are needed to achieve a small rejection probability. Although versions 1 and 2 of the sampler are based on the same proposal, the changes to their bounds are seen to differ due to the randomness in knot selection.

```{r}
#| out-width: 60%
#| fig-cap: |
#|   Empirical distribution of draws (solid) versus target (dashed) for VMF
#|   example with linear majorizer.
#| label: fig-vmf-linear-draws
#| echo: false

data.frame(x = out3$draws) %>%
	ggplot() +
	geom_density(aes(x = x), col = "black") +
	geom_function(fun = d_target, args = list(kappa = 5, d = 4), lty = 2) +
	ylab("Empirical Density") +
	theme_minimal()
```

```{r}
#| out-width: 60%
#| fig-cap: |
#|   Refinement for VMF example with constant majorizer.
#| label: fig-vmf-refine
#| echo: false

df1 = data.frame(sampler = "v1", bdd = exp(out1$lbdd)) %>% 
	mutate(N = row_number())
df2 = data.frame(sampler = "v2", bdd = exp(out2$lbdd)) %>% 
	mutate(N = row_number())
df3 = data.frame(sampler = "v3", bdd = exp(out3$lbdd)) %>% 
	mutate(N = row_number() + 1)
df = bind_rows(df1, df2, df3) %>% mutate(N = as.integer(N))

ggplot(df) +
	geom_line(aes(x = N, y = bdd, color = sampler), lwd = 1.02, alpha = 1) +
	geom_point(aes(x = N, y = bdd, pch = sampler, color = sampler),
		cex = 3, alpha = 1) +
	scale_x_continuous(breaks = df$N[df$N %% 2 == 0]) +
	scale_y_continuous(n.breaks = 10) +
	xlab("N") +
	ylab("Bound") +
	theme_minimal()
```

# Example: Privacy Protected Data {#sec-privacy}

Suppose an official statistics agency collects sensitive data from respondents in a survey or census. The agency then seeks to disseminate a tabulation to data users while protecting the privacy of the respondents. The agency can consider adding random noise $\gamma$ to the tabulation $Y$ and releasing $Z = Y + \gamma$. This setting is the basis of a modeling scenario considered by @Direct-Sampling-DAS-2021, @DP-Simulation-2022, and @DAS-Methods-2025 where the distribution for $\gamma$ has been selected by the agency and the data have been disseminated; the objective in these works is to carry out inference on the unobserved $Y$ given observed $Z = z$. The field of differential privacy studies the design of noise mechanisms (including distributions of $\gamma$ in the present setting) to satisfy mathematical criteria for privacy [e.g., @Dwork-Roth-2014]. @Abowd-EtAl-2022 describe work by the U.S. Census Bureau to implement differential privacy in the release of data from the decennial census.

Our present motivation is to focus on a specific target that arises in analysis of the released $z$. @sec-privacy-within-gibbs considers a data analysis situation with $n$ subjects and a regression underlying $y$. To reduce the amount of computation under repeated use within a Gibbs sampler, we will consider proposals that persist over the MCMC chain and get updated with the current state.

Suppose $Y$ and $\gamma$ are independently distributed with $Y \sim \text{Lognormal}(\mu, \sigma^2)$ and $\gamma \sim \text{N}(0, \lambda^2)$. The variance $\lambda^2$ is assumed to be known and provided with the noisy data. Such transparency into the noise mechanism is often featured in differential privacy.

Suppose the target distribution is the conditional $[Y \mid Z = z, \mu, \sigma^2]$ whose density is given by
<!-- -->
\begin{align}
f(y \mid z, \mu, \sigma^2)
%
&\propto \frac{1}{\lambda \sqrt{2\pi}} \exp\left\{
-\frac{1}{2\lambda^2} (z - y)^2
\right\} \cdot
\frac{1}{y\sigma \sqrt{2\pi}} \exp\left\{
-\frac{1}{2\sigma^2} (\log y - \mu)^2
\right\} \ind(y > 0)
\nonumber \\
%
&\propto
\underbrace{\frac{1}{\lambda \sqrt{2\pi}} \exp\left\{
-\frac{1}{2\lambda^2} (z - y)^2
\right\}}_{g(y)} \cdot
\underbrace{\frac{1}{y} \exp\left\{
-\frac{1}{2\sigma^2} (\log y - \mu)^2
\right\} \ind(y > 0)}_{w(y)}.
\label{eqn:privacy}
\end{align}
<!-- -->
We have decomposed $f$ into weight function
<!-- -->
\begin{math}
w(y) = \frac{1}{y} \exp\left\{ -\frac{1}{2\sigma^2} (\log y - \mu)^2 \right\}
\ind\{ y \in (0,\infty) \},
\end{math}
<!-- -->
from the Lognormal component---dropping some of the terms in its normalizing constant---and base distribution whose density $g$ is $\text{N}(z, \lambda^2)$. The conditional $[\vec{y} \mid \sigma^2, \vec{\mu}, \vec{z}]$ in \eqref{eqn:privacy} would be encountered in a Gibbs sampler for the posterior distribution $[\vec{y}, \vec{\mu}, \sigma^2 \mid \vec{z}]$, based on an observed sample $\vec{z} = (z_1, \ldots, z_n)$ with augmented data $\vec{y} = (y_1, \ldots, y_n)$. The remaining conditionals $[\vec{\mu} \mid \sigma^2, \vec{y}, \vec{z}]$ and $[\sigma^2 \mid \vec{\mu}, \vec{y}, \vec{z}]$ may be straightforward to generate if a convenient prior distribution is selected; therefore, we will focus on \eqref{eqn:privacy}.


Before proceeding, let us fix the following values for the parameters.

```{r}
#| prompt: true
mu = 5
sigma = sqrt(0.5)
lambda = 10
```

Jointly draw values of $Y$ and $Z$ from the model; $Z$ will be observed by users of the data while $Y$ is unobserved and the objective for inference.

```{r}
#| prompt: true
source("examples/privacy/functions.R")
y_true = rlnorm(1, mu, sigma)
z = rnorm(1, y_true, lambda)
vws::printf("y_true = %g and z = %g\n", y_true, z)
```

Coding the target density is helpful to evaluate the distribution of the draws. To compute the normalizing constant, we use Hermite quadrature via the `gauss.quad` function in the `statmod` package [@Smyth-2005]. The integral $\psi = \int_{-\infty}^\infty q(x) e^{-x^2} dx$ is approximated as $\psi \approx \sum_{j=1}^Q \omega_j q(x_j)$ using quadrature points $x_1, \ldots, x_Q$ and weights $\omega_1, \ldots, \omega_Q$; to identify the function $q$, we have
<!-- -->
\begin{align*}
\psi &= \int_{-\infty}^\infty w(y) g(y) dy \\
%
&= \int_{-\infty}^\infty \ind\{y > 0\} \cdot \frac{1}{y}
\exp\left\{ -\frac{1}{2\sigma^2} (\log y - \mu)^2 \right\}
\frac{1}{\lambda \sqrt{2\pi}} \exp\left\{ -\frac{1}{2\lambda^2} (z - y)^2 \right\} dy \\
%
&= \int_{-\infty}^\infty \ind\left( z > \sqrt{2} \lambda x \right)
\cdot \frac{1}{z - \sqrt{2} \lambda x}
\exp\left\{ -\frac{1}{2\sigma^2} \left[\log\left( z - \sqrt{2} \lambda x \right) - \mu \right]^2 \right\}
\frac{1}{\sqrt{\pi}} e^{-x^2} dx,
\end{align*}
<!-- -->
using the transformation $y = z - \sqrt{2} \lambda x$, so that
<!-- -->
\begin{align*}
q(x) = \ind\left( z > \sqrt{2} \lambda x \right)
\cdot \frac{1}{z - \sqrt{2} \lambda x}
\exp\left\{ -\frac{1}{2\sigma^2} \left[\log\left( x - \sqrt{2} \lambda x \right) - \mu \right]^2 \right\}
\frac{1}{\sqrt{\pi}}
\end{align*}
<!-- -->
is used with `gauss.quad`.

We will consider three variations of a VWS sampler, progressing from easier-to-implement to more computationally efficient. @sec-privacy-const-default considers a constant majorizer where the constant for each region is obtained by numerical optimization. @sec-privacy-const-custom replaces numerical optimization with code for a closed-form solution. @sec-privacy-linear makes use of a linear majorizer which subclasses the abstract `Region` class. All codes for this example are in the folder `examples/privacy`. The function `d_target` defined in `examples/privacy/functions.R` computes the target density.

## Constant Majorizer with Numerical Optimization {#sec-privacy-const-default}

We may adapt the code from @sec-vmf-const-default to the target density in the present problem. The present weight function may be coded as follows.

```{cpp}
const vws::dfdb& w = [&](double x, bool log = true) {
	double out = R_NegInf;
	if (x > 0) {
		double sigma2 = std::pow(sigma, 2)
		out = -std::log(x) - std::pow(std::log(x) - mu, 2) / (2 * sigma2);
	}
	return log ? out : std::exp(out);
};
```

The base distribution's density, CDF, and quantile function are coded as lambdas and packaged into a `univariate_helper` object. Here we can make use of the implementations of the normal distribution in R's API.

```{cpp}
fntl::density df = [&](double x, bool log = false) {
	return R::dnorm(x, z, lambda, log);
};
fntl::cdf pf = [&](double q, bool lower = true, bool log = false) {
	return R::pnorm(q, z, lambda, lower, log);
};
fntl::quantile qf = [&](double p, bool lower = true, bool log = false) {
	return R::qnorm(p, z, lambda, lower, log);
};

vws::univariate_helper helper(df, pf, qf);
```

The following R snippet builds the C++ code and invokes the sampling function.

```{r}
#| prompt: true
Rcpp::sourceCpp("examples/privacy/ln-norm-v1.cpp")
out1 = r_ln_norm_v1(n = 1000, z, mu, sigma, lambda, N = 50, tol = 0.10)
head(out1$draws)
```

@fig-privacy-const-refine plots the bound for the rejection probability during the refinement process, which is captured in the variable `out$lbdd`. @fig-privacy-const-draws plots the empirical distribution of the draws and compares them to the density. It displays an interval based on the $0.025$ and $0.975$ quantiles of the distribution $[Y \mid Z = z]$ approximated from the empirical quantiles of the draws. The value of the observed $z$ and the latent $y$ are also highlighted for reference.

```{r}
#| out-width: 60%
#| fig-cap: |
#|   Refinement for lognormal-normal example with constant majorizer.
#| label: fig-privacy-const-refine
#| echo: false

df = data.frame(bdd = exp(out1$lbdd)) %>%
	mutate(step = row_number() - 1)
ggplot(df) +
	geom_line(aes(x = step, y = bdd)) +
	scale_x_continuous(breaks = df$step[df$step %% 2 == 0]) +
	scale_y_continuous(n.breaks = 10) +
	xlab("Step") +
	ylab("Bound") +
	theme_minimal()
```

```{r}
#| out-width: 60%
#| fig-cap: |
#|   Empirical distribution of draws (solid black) versus target (dashed black)
#|   for lognormal-normal example with constant majorizer. Observed value of
#|   $z$ (blue line) and latent value of $y$ (red line) are displayed along
#|   with a 95\% interval (blue ribbon) based on draws from $[Y \mid Z = z]$.
#| label: fig-privacy-const-draws
#| echo: false

interval_lo = quantile(out1$draws, probs = 0.025)
interval_hi = quantile(out1$draws, probs = 0.975)

data.frame(x = out1$draws) %>%
	ggplot() +
	geom_density(aes(x = x), col = "black") +
	geom_function(fun = d_target, lty = 2,
		args = list(mu = mu, sigma = sigma, z = z, lambda = lambda)) +
	annotate("rect", xmin = interval_lo, xmax = interval_hi, ymin = 0,
		ymax = Inf, alpha = 0.1, fill = "blue") +
	geom_vline(xintercept = z, col = "blue", lwd = 1.05) +
	geom_vline(xintercept = y_true, col = "red", lwd = 1.05) +
	ylab("Empirical Density") +
	theme_minimal()
```

## Constant Majorizer with Custom Optimization {#sec-privacy-const-custom}

The log-weight function can be both maximized and minimized in closed form. Coding it explicitly reduces computational overhead and avoids convergence issues from numerical optimization. This section will demonstrate how to override the default `optimize` method of `real_const_region`.

We have first derivative
<!-- -->
\begin{align*}
\frac{d}{dy} \log w(y) = -\frac{1}{y}\left(
1 + \frac{\log y - \mu}{\sigma^2}
\right),
\end{align*}
<!-- -->
for $y \in (0, \infty)$. Let $y^* = \exp(\mu - \sigma^2)$; it is seen that $\frac{d}{dy} \log w(y)$ is positive when $y < y^*$, negative when $y > y^*$, and takes value zero at $y = y^*$. Therefore, $\log w(y)$ is unimodal and $y^*$ maximizes $\log w(y)$ with
<!-- -->
\begin{align*}
\log w(y^*)
&= -\log[\exp(\mu - \sigma^2)]  - \frac{(\log [\exp(\mu - \sigma^2)] - \mu)^2}{2 \sigma^2} \\
&= -\mu + \sigma^2 - \frac{(\sigma^2)^2}{2 \sigma^2} \\
&= -\mu + \sigma^2 / 2.
\end{align*}
<!-- -->
Therefore, on a region $\mathscr{D}_j = (\alpha_{j-1}, \alpha_j]$ where both endpoints are smaller than $y^*$, the maximum of $\log w(y)$ is $\log w(\alpha_j)$ and the minimum is $\log w(\alpha_{j-1})$. On the other hand, for a region where both endpoints are larger than $y^*$, the maximum of $\log w(y)$ is $\log w(\alpha_{j-1})$ and the minimum is $\log w(\alpha_j)$. @fig-privacy-weight displays a plot of $\log w(y)$ with our selected $\mu$ and $\sigma^2$ values.

```{r}
#| prompt: true
y_star = exp(mu - sigma^2)
w_star = -mu + sigma^2 / 2
vws::printf("Maximizer y = %g obtains value log w(%g) = %g.\n",
	y_star, y_star, w_star)
```

```{r}
#| out-width: 60%
#| fig-cap: |
#|   Weight function for Lognormal-Normal example on the log-scale, with the
#|   maximizer $y^* = \exp(\mu - \sigma^2)$ highlighted.
#| label: fig-privacy-weight
#| echo: false

xlim = y_star + c(-10,10)
w_one = function(y, log = T) {
	out = -Inf
	if (y > 0) { out = -log(y) - (log(y) - mu)^2 / (2*sigma^2) }
	if (log) { return(out) } else { return(exp(out)) }
}
w = Vectorize(w_one, vectorize.args = "y")
ggplot() +
	geom_function(fun = w, xlim = xlim) +
	geom_vline(xintercept = y_star, lty = 2) +
	geom_hline(yintercept = w_star, lty = 2) +
	xlab("y") +
	ylab(expression("log w(y)")) +
	theme_minimal()
```

The maximizer may be coded as follows.

```{cpp}
const vws::optimizer& maxopt = [&](const vws::dfdb& w, double lo,
	double hi, bool log)
{
	double y_star = exp(mu - sigma2);

	double y = y_star;
	if (y_star > hi) {
		y = hi;
	} else if (y_star < lo) {
		y = lo;
	}

	double out = w(y, true);
	return log ? out : exp(out);
};
```

Here is code for the minimizer.

```{cpp}
const vws::optimizer& minopt = [&](const vws::dfdb& w, double lo,
	double hi, bool log)
{
	double lwa = w(lo, true);
	double lwb = w(hi, true);
	double out = std::min(lwa, lwb);
	return log ? out : exp(out);
};
```

We can construct the proposal using a single `real_const_region` that represents the support; the `maxopt` and `minopt` arguments are specified here.

```{cpp}
vws::real_const_region supp(0, R_PosInf, w, helper, maxopt, minopt);
```

The following R snippet builds the C++ code and invokes the sampling function.

```{r}
#| prompt: true
Rcpp::sourceCpp("examples/privacy/ln-norm-v2.cpp")
out2 = r_ln_norm_v2(n = 1000, z, mu, sigma, lambda, N = 50, tol = 0.10)
head(out2$draws)
```

## Linear Majorizer {#sec-privacy-linear}

We can obtain a linear majorizer by noting the convexity of $\log w(y)$. Its second derivative is seen to be
<!-- -->
\begin{align*}
&\frac{d^2}{dy^2} \log w(y) = -\frac{1}{y^2} \left(
1 + \frac{\log y - \mu - 1}{\sigma^2}
\right),
\end{align*}
<!-- -->
which has root $y_0 = \exp\{ \mu - \sigma^2 + 1 \}$. The weight function $w$ is log-concave for $y < y_0$ and log-convex for $y > y_0$. This is plotted in @fig-privacy-weight-convexity. We will assume that there are two initial regions $\mathscr{D}_1 = (0, y_0]$ and $\mathscr{D}_2 = (y_0, \infty]$ so that all partitions considered thereafter consist of regions on which $\log w(y)$ is entirely concave or convex.

```{r}
#| prompt: true
y0 = exp(mu - sigma^2 + 1)
printf("Convexity changes at y = %g.\n", y0)
```

```{r}
#| out-width: 60%
#| fig-cap: |
#|   Weight function for Lognormal-Normal example on the log-scale,
#|   highlighting $y_0 = \exp(\mu - \sigma^2 + 1)$ where there is a change in
#|   convexity.
#| label: fig-privacy-weight-convexity
#| echo: false

xlim = c(0, 2000)
ggplot() +
	geom_function(fun = w, xlim = xlim) +
	geom_vline(xintercept = y0, lty = 2) +
	xlab("y") +
	ylab(expression("log w(y)")) +
	theme_minimal()
```

As in @sec-vmf-linear, we can majorize $w$ with exponentiated linear functions of the form $\overline{w}_j(y) = \exp\{ \overline{\beta}_{0j} + \overline{\beta}_{1j} y \}$. Here, we also assume a minorizer of the form $\underline{w}_j(y) = \exp\{ \underline{\beta}_{0j} + \underline{\beta}_{1j} y \}$. Before proceeding, the following remark gives an integral that will be used several times.

::: {#rem-normal-mgf-integral}
Let $\phi(\cdot \mid \mu, \sigma^2)$ and $\Phi(\cdot \mid \mu, \sigma^2)$ be the density and CDF of $Y \sim \text{N}(\mu, \sigma^2)$, respectively. If $a < b$ are scalars (possibly infinite), then
<!-- -->
\begin{align*}
\int_a^b e^{ty} \phi(y \mid \mu, \sigma^2) dy
&= \exp(\mu t + t^2 \sigma^2 / 2) \left\{ 
\Phi(b \mid \mu + t \sigma^2, \sigma^2) - \Phi(a \mid \mu + t \sigma^2, \sigma^2)
\right\} \\
&= \exp(\mu t + t^2 \sigma^2 / 2)
\Prob(a < T \leq b).
\end{align*}
<!-- -->
where $T \sim \text{N}(\mu + t \sigma^2, \sigma^2)$. The special case $a = -\infty$ and $b = \infty$ yields the moment-generating function $M(t) = \exp(\mu t + t^2 \sigma^2 / 2)$ of $Y$. \remarkqed
:::

@rem-normal-mgf-integral yields the expressions
<!-- -->
\begin{align*}
\overline{\xi}_j
&= \int_{\alpha_{j-1}}^{\alpha_j} \overline{w}_j(y) g(y) dy \\
%
&= \exp\{ \overline{\beta}_{0j} \} \int_{\alpha_{j-1}}^{\alpha_j} 
\exp\{ \overline{\beta}_{1j} y \}
\frac{1}{\lambda \sqrt{2\pi}}
\exp\left\{ -\frac{1}{2 \lambda^2} (y - z)^2 \right\} dy \\
%
&= \exp\left\{
\overline{\beta}_{0j} + z \overline{\beta}_{1j}
+ \frac{1}{2} \overline{\beta}_{1j}^2 \lambda^2
\right\}
\Prob(\alpha_{j-1} < \overline{T} \leq \alpha_j),
\end{align*}
<!-- -->
where $\overline{T} \sim \text{N}(z + \overline{\beta}_{1j} \lambda^2, \lambda^2)$, and similarly,
<!-- -->
\begin{align*}
\underline{\xi}_j
&= \exp\left\{
\underline{\beta}_{0j} + z \underline{\beta}_{1j}
+ \frac{1}{2} \underline{\beta}_{1j}^2 \lambda^2
\right\}
\Prob(\alpha_{j-1} < \underline{T} \leq \alpha_j),
\end{align*}
<!-- -->
where $\underline{T} \sim \text{N}(z + \underline{\beta}_{1j} \lambda^2, \lambda^2)$. The proposal $h$ is a finite mixture $h(y) = \sum_{j=1}^N \pi_j g_j(y)$ with
<!-- -->
\begin{align*}
g_j(y)
&= \overline{w}_j(y) g(y) \ind\{ y \in \mathscr{D}_j \} / \overline{\xi}_j \\
%
&= \exp\{ \overline{\beta}_{0j} + \overline{\beta}_{1j} y \}
\frac{1}{\lambda \sqrt{2\pi}}
\exp\left\{ -\frac{1}{2 \lambda^2} (y - z)^2 \right\}
\cdot \ind(\alpha_{j-1} < y \leq \alpha_j) / \overline{\xi}_j  \\
%
&= \frac{
\frac{1}{\lambda \sqrt{2\pi}}
\exp\left\{ -\frac{1}{2 \lambda^2} (y - z - \lambda^2 \overline{\beta}_{1j})^2 \right\}
}{
\Prob(\alpha_{j-1} < \overline{T} \leq \alpha_j)
}
\cdot \ind(\alpha_{j-1} < y \leq \alpha_j),
\end{align*}
<!-- -->
which is the density of $\text{N}(z + \lambda^2 \overline{\beta}_{1j}, \lambda^2)$ truncated to the interval $(\alpha_{j-1}, \alpha_j]$. @rem-expansion-point is used to select the expansion point $c$ to determine coefficients for the majorizer in the log-concave case and the minorizer in the log-convex case, with
<!-- -->
\begin{align*}
M_j(s) = \exp(z s + s^2 \lambda^2 / 2) 
\frac{
\Phi(\alpha_j \mid z + s \lambda^2, \lambda^2) - \Phi(\alpha_{j-1} \mid z + s \lambda^2, \lambda^2)
}{
\Phi(\alpha_j \mid z, \lambda^2) - \Phi(\alpha_{j-1} \mid z, \lambda^2)
}.
\end{align*}
<!-- -->
Adapting the code from @sec-vmf-linear to the present problem, we implement `LinearVWSRegion` as a subclass of the abstract `Region` base class (@sec-api-region). See the file `examples/privacy/LinearVWSRegion.h`. The code in `examples/privacy/ln-norm-v3.cpp` instantiates a proposal from this region type and proceeds with rejection sampling. This code is displayed below.

``` {.cpp include="examples/privacy/ln-norm-v3.cpp" code-line-numbers="true"}
```

The following R snippet builds the C++ code and invokes the sampling function.

```{r}
#| prompt: true
Rcpp::sourceCpp("examples/privacy/ln-norm-v3.cpp")
out3 = r_ln_norm_v3(n = 1000, z, mu, sigma, lambda, lo = 1e-8, 1e8,
	N = 50, tol = 0.10)
head(out3$draws)
```

@fig-privacy-linear-draws plots the empirical distribution of the draws and compares them to the density. @fig-privacy-refine plots the bound for the rejection probability for this sampler after each step of refining, along with those from the previous two versions. A substantial improvement in efficiency is seem here; fewer regions are needed to achieve a small rejection probability. Although versions 1 and 2 of the sampler are based on the same proposal, the changes to their bounds are seen to differ due to randomness in knot selection.

```{r}
#| out-width: 60%
#| fig-cap: |
#|   Empirical distribution of draws (solid) versus target (dashed) for
#|   lognormal-normal example with constant majorizer.
#| label: fig-privacy-linear-draws
#| echo: false

data.frame(x = out1$draws) %>%
	ggplot() +
	geom_density(aes(x = x), col = "black") +
	geom_function(fun = d_target, lty = 2,
		args = list(mu = mu, sigma = sigma, z = z, lambda = lambda)) +
	ylab("Empirical Density") +
	theme_minimal()
```

```{r}
#| out-width: 60%
#| fig-cap: |
#|   Refinement for lognormal-normal example with constant majorizer.
#| label: fig-privacy-refine
#| echo: false

df1 = data.frame(sampler = "v1", bdd = exp(out1$lbdd)) %>% 
	mutate(N = row_number())
df2 = data.frame(sampler = "v2", bdd = exp(out2$lbdd)) %>% 
	mutate(N = row_number())
df3 = data.frame(sampler = "v3", bdd = exp(out3$lbdd)) %>% 
	mutate(N = row_number() + 1)
df = bind_rows(df1, df2, df3) %>% mutate(N = as.integer(N))

ggplot(df) +
	geom_line(aes(x = N, y = bdd, color = sampler), lwd = 1.02, alpha = 1) +
	geom_point(aes(x = N, y = bdd, pch = sampler, color = sampler),
		cex = 3, alpha = 1) +
	scale_x_continuous(breaks = df$N[df$N %% 2 == 0]) +
	scale_y_continuous(n.breaks = 10) +
	xlab("N") +
	ylab("Bound") +
	theme_minimal()
```

# Example: Bessel Count Distribution {#sec-bessel}

Where Sections [-@sec-vmf] and [-@sec-privacy] considered continuous distributions, we now generate from a discrete target. The Bessel distribution [@Yuan-Kalbfleisch-2000] has density
<!-- -->
\begin{align*}
f(x) = \frac{(\lambda/2)^{2x+\nu}}{I_v(\lambda) \cdot x! \cdot \Gamma(x + \nu + 1)}
\cdot \ind\{ x \in \mathbb{N} \}, \quad
\nu > -1, \quad
\lambda > 0,
\end{align*}
<!-- -->
where $\mathbb{N} = \{ 0, 1, 2, \ldots \}$ is the set of nonnegative integers and where 
<!-- -->
\begin{align*}
I_v(\lambda) = \sum_{x=0}^\infty \frac{(\lambda/2)^{2x+\nu}}{x! \cdot \Gamma(x + \nu + 1)}
\end{align*}
<!-- -->
is a modified Bessel function of the first kind. Let $X \sim \text{Bessel}(\lambda, \nu)$ denote a random variable with density $f$. @Devroye-2002 develops a rejection sampling method to generate draws of $X$ using properties of the distribution. Here we implement several VWS samplers which do not require as much insight. Our approach is to decompose the density into
<!-- -->
\begin{align*}
f(x)
&= \frac{(\lambda/2)^{2x+\nu}}{I_v(\lambda) \cdot x! \cdot \Gamma(x + \nu + 1)}
\cdot \ind\{ x \in \mathbb{N} \} \\
%
&\propto 
\underbrace{\frac{1}{\Gamma(x + \nu + 1)}}_{w(x)}
\underbrace{\frac{(\lambda/2)^{2x} e^{-\lambda^2 / 4}}{x!} \ind\{ x \in \mathbb{N} \}}_{g(x)},
\end{align*}
<!-- -->
disregarding several terms in the normalizing constant, so that $g$ is the density of a $\text{Poisson}(\lambda^2 / 4)$ distribution and the weight function is specified by $\log w(x) = -\log \Gamma(x + \nu + 1)$.

Before proceeding, let us fix values of $\lambda$ and $\nu$. Let us also load some useful R functions from the script `examples/bessel/bessel.R`; primarily, `d_bessel` which computes the target density.

```{r}
#| prompt: true
source("examples/bessel/bessel.R")
lambda = 10
nu = 2
```

```{r}
#| echo: false
source("examples/common/plots.R")
xseq = seq(0, 15)
fseq = d_bessel(xseq, lambda, nu)
lwseq = w(xseq, log = T)
```

Three variations of VWS samplers are considered in subsections. @sec-bessel-const-default uses a constant majorizer obtained by numerical optimization. @sec-bessel-const-custom replaces numerical optimization with a closed-form solution. @sec-bessel-linear uses a linear majorizer and subclasses the abstract `Region` class. All codes for this example are in the folder `examples/bessel`. The function `d_bessel` defined in `examples/bessel/bessel.R` computes the target density.

## Constant Majorizer with Numerical Optimization {#sec-bessel-const-default}

The weight function may be coded as a lambda in C++ as follows.

```{cpp}
const vws::dfdb& w = [&](double x, bool log = true) {
	double out = -std::lgamma(x + nu + 1);
	return log ? out : std::exp(out);
};
```

The base distribution’s density, CDF, and quantile function may be supplied as a `univariate_helper` object using Poisson functions in the R API.

```{cpp}
fntl::density df = [&](double x, bool log = false) {
	return R::dpois(x, mean, log);
};
fntl::cdf pf = [&](double q, bool lower = true, bool log = false) {
	return R::ppois(q, mean, lower, log);
};
fntl::quantile qf = [&](double p, bool lower = true, bool log = false) {
	return R::qpois(p, mean, lower, log);
};

vws::univariate_helper helper(df, pf, qf);
```

The following R snippet builds the C++ code and invokes the sampling function.

```{r}
#| prompt: true
Rcpp::sourceCpp("examples/bessel/bessel-v1.cpp")
out1 = r_bessel_v1(n = 1000, lambda, nu, N = 50, tol = 0.10)
head(out1$draws)
```

@fig-bessel-const-refine plots the bound for the rejection probability during the refinement process, which is captured in the variable `out$lbdd`. @fig-bessel-const-draws plots the empirical distribution of the draws and compares them to the target probability mass function (pmf).

```{r}
#| out-width: 60%
#| fig-cap: |
#|   Refinement for lognormal-normal example with constant majorizer.
#| label: fig-bessel-const-refine
#| echo: false

df = data.frame(bdd = exp(out1$lbdd)) %>%
	mutate(step = row_number() - 1)
ggplot(df) +
	geom_line(aes(x = step, y = bdd)) +
	scale_x_continuous(breaks = df$step[df$step %% 2 == 0]) +
	scale_y_continuous(n.breaks = 10) +
	xlab("Step") +
	ylab("Bound") +
	theme_minimal()
```

```{r}
#| out-width: 60%
#| fig-cap: |
#|   Empirical distribution of draws (bars) versus target (points) for Bessel
#|   example with constant majorizer.
#| label: fig-bessel-const-draws
#| echo: false

plot_pmf(out1$draws) +
	geom_point(data = data.frame(x = xseq, y = fseq), aes(x,y)) +
	scale_x_continuous(breaks = xseq[xseq %% 2 == 0])
```

## Constant Majorizer with Custom Optimization {#sec-bessel-const-custom}

We can make several improvements to the default numerical majorizer by coding it explicitly. Note that $\log w(x) = -\log \Gamma(x + \nu + 1)$ is a decreasing function in $x$ so that numerical optimization is not necessary. Therefore,
<!-- -->
\begin{displaymath}
\log w(\alpha_j) \leq \log w(x) \leq \log w(\alpha_{j-1}), \quad
x \in \mathscr{D}_j.
\end{displaymath}
<!-- -->
If the endpoints are not integers, we can improve the bounds by restricting consideration to integers, using $\lfloor \alpha_{j-1} \rfloor + 1$ and $\lfloor \alpha_j \rfloor$ to obtain the maximizer and minimizer, respectively. The following functions put this into practice.

```{cpp}
const vws::optimizer& maxopt = [](const vws::dfdb& w, double lo,
	double hi, bool log)
{
	if (lo < 0 && hi < 0) { Rcpp::stop("Did not code this case"); }
	double x = (lo < 0) ? 0 : std::floor(lo) + 1;
	double out = w(x, true);
	return log ? out : std::exp(out);
};

const vws::optimizer& minopt = [](const vws::dfdb& w, double lo,
	double hi, bool log)
{
	if (lo < 0 && hi < 0) { Rcpp::stop("Did not code this case"); }
	double x = std::isinf(hi) ? hi : std::floor(hi);
	double out = w(x, true);
	return log ? out : std::exp(out);
};
```

We can construct the proposal using a single `int_const_region` that represents the support; the `maxopt` and `minopt` arguments are specified here. Also note that we have specified a lower endpoint $\alpha_0$ below zero to ensure that the support $(\alpha_0, \alpha_N]$ contains zero.

```{cpp}
vws::int_const_region supp(-0.1, R_PosInf, w, helper, maxopt, minopt);
```

The following R snippet builds the C++ code and invokes the sampling function.

```{r}
#| prompt: true
Rcpp::sourceCpp("examples/bessel/bessel-v2.cpp")
out2 = r_bessel_v2(n = 1000, lambda, nu, N = 50, tol = 0.10)
head(out2$draws)
```

## Linear Majorizer {#sec-bessel-linear}

Our decomposition of $f$ into the given $w$ and $g$ lends itself well to a linear majorizer because the function $\log w(x) = -\log \Gamma(x + \nu + 1)$ is log-concave. In fact, a theorem attributed to @Bohr-Mollerup-1922 states that the function $q(x) = \Gamma(x)$ is uniquely characterized as the only function on $(0, \infty)$ which is log-convex with $q(1) = 1$ and $q(x+1) = x q(x)$. The following remark is also useful to formulate the linear majorizer.

::: {#rem-incgamma}
By rearranging equation 6.5.13 of @Abramowitz-Stegun-1972, the cdf of $T \sim \text{Poisson}(\lambda)$ can be written as
<!-- -->
\begin{displaymath}
F(x \mid \lambda)
= e^{-\lambda} \sum_{k=0}^n \frac{\lambda^k}{k!}
= \frac{\Gamma(x + 1, \lambda)}{\Gamma(x + 1)},
\quad x \in \{ 1, 2, \ldots \},
\end{displaymath}
<!-- -->
where $\Gamma(c,z) = \int_z^\infty t^{c-1} e^{-t} dt$ and $\Gamma(c) = \Gamma(c,0)$. Therefore, for $a, b \in \mathbb{R}$ with $0 \leq a \leq b$ and $T \sim \text{Poisson}(\lambda)$,
<!-- -->
\begin{align*}
\Prob(a < T \leq b)
&= \Prob(\lfloor a \rfloor + 1 \leq X \leq \lfloor b \rfloor) \\
&= \sum_{x = \lfloor a \rfloor + 1}^{\lfloor b \rfloor}
\frac{e^{-\lambda} \lambda^x}{x!} \\
&= \sum_{x = 0}^{\lfloor b \rfloor} \frac{e^{-\lambda} \lambda^x}{x!}
- \sum_{x = 0}^{\lfloor a \rfloor} \frac{e^{-\lambda} \lambda^x}{x!} \\
&= \frac{\Gamma(\lfloor b \rfloor + 1, \lambda)}{\Gamma(\lfloor b \rfloor + 1)}
- \frac{\Gamma(\lfloor a \rfloor + 1, \lambda)}{\Gamma(\lfloor a \rfloor + 1)}.
\end{align*}
<!-- -->
To verify the first equality, if $a$ is an integer then
<!-- -->
\begin{align*}
\Prob(a < T \leq b)
&= \Prob(T \in  \{a + 1, \ldots, \lfloor b \rfloor \} ) \\
&= \Prob(T \in  \{ \lfloor a \rfloor + 1, \ldots, \lfloor b \rfloor \})
= \Prob(\lfloor a \rfloor + 1 < T \leq \lfloor b \rfloor);
\end{align*}
<!-- -->
otherwise, if $a$ is not an integer then
<!-- -->
\begin{align*}
\Prob(a < T \leq b)
&= \Prob(T \in  \{ \lceil a \rceil, \ldots, \lfloor b \rfloor \} ) \\
&= \Prob(T \in  \{ \lfloor a \rfloor + 1, \ldots, \lfloor b \rfloor \})
= \Prob(\lfloor a \rfloor + 1 < T \leq \lfloor b \rfloor).
\end{align*}
<!-- -->
\remarkqed
:::

We have the expression
<!-- -->
\begin{align*}
\overline{\xi}_j
&= \int_{\alpha_{j-1}}^{\alpha_j} \overline{w}_j(x) g(x) d\nu(x) \\
%
&= \exp\{ \overline{\beta}_{0j} \}
\sum_{k \in (\alpha_{j-1}, \alpha_j]}
\exp\{ \overline{\beta}_{1j} x \}
\frac{(\lambda^2 / 4)^x e^{-\lambda^2 / 4} }{x!} \\
%
%&= \exp\{ \overline{\beta}_{0j} \}
%\exp\{ -\lambda^2 / 4 \}
%\sum_{k \in (\alpha_{j-1}, \alpha_j]}
%\frac{1}{x!}
%\left[ e^{ \frac{\lambda^2}{4}\overline{\beta}_{1j} } \right]^x \\
%
&= \exp\left(
\overline{\beta}_{0j} - \frac{\lambda^2}{4}
+ \frac{\lambda^2}{4} e^{ \overline{\beta}_{1j} } \right)
\Prob(\alpha_{j-1} < \overline{T} \leq \alpha_j),
\end{align*}
<!-- -->
where $\overline{T} \sim \text{Poisson}(e^{ \overline{\beta}_{1j}} \lambda^2 / 4)$, and similarly,
<!-- -->
\begin{align*}
\underline{\xi}_j
%
&= \exp\left(
\overline{\beta}_{0j} - \frac{\lambda^2}{4}
+ \frac{\lambda^2}{4} e^{ \overline{\beta}_{1j} } \right)
\Prob(\alpha_{j-1} < \underline{T} \leq \alpha_j).
\end{align*}
<!-- -->
where $\underline{T} \sim \text{Poisson}(e^{ \underline{\beta}_{1j}} \lambda^2 / 4)$. The proposal $h$ is a finite mixture $h(x) = \sum_{j=1}^N \pi_j g_j(x)$ with
<!-- -->
\begin{align*}
g_j(x)
&= \overline{w}_j(x) g(x) \ind\{ x \in \mathscr{D}_j \} / \overline{\xi}_j \\
%
&= 
\frac{
\exp\{ \overline{\beta}_{0j} + \overline{\beta}_{1j} x \}
(\lambda^2 / 4)^x e^{-\lambda^2 / 4} / x!
}{
\exp\left(
\overline{\beta}_{0j} - \frac{\lambda^2}{4}
+ \frac{\lambda^2}{4} e^{ \overline{\beta}_{1j} } \right)
\Prob(\alpha_{j-1} < \overline{T} \leq \alpha_j)
} \ind\{ x \in \mathscr{D}_j \} \\
%
&=
\frac{1}{x!}
\left(\frac{\lambda^2}{4} e^{ \overline{\beta}_{1j} } \right)^x
\exp\left( -\frac{\lambda^2}{4} e^{ \overline{\beta}_{1j} } \right)
\frac{1}{
\Prob(\alpha_{j-1} < \overline{T} \leq \alpha_j)
}
\ind(\alpha_{j-1} < x \leq \alpha_j),
\end{align*}
<!-- -->
which is the density of $\text{Poisson}(\frac{\lambda^2}{4} e^{ \overline{\beta}_{1j}})$ truncated to the interval $(\alpha_{j-1}, \alpha_j]$. @rem-expansion-point is used to select the expansion point $c$ to determine coefficients for the majorizer in the log-concave case, with
<!-- -->
\begin{align*}
M_j(s)
&= \int_{\alpha_{j-1}}^{\alpha_j} e^{sx} g(x) d\nu(x) \\
&= \exp( -\lambda^2 / 4 ) \exp(e^s \lambda^2 / 4 )
\sum_{x \in (\alpha_{j-1}, \alpha_j]}
\frac{1}{x!} \exp\{ -e^s \lambda^2 / 4 \}
(e^s \lambda^2 / 4)^x \\
&= \exp\left\{ -\frac{\lambda^2}{4} + e^s \frac{\lambda^2}{4} \right\}
\Prob(\alpha_{j-1} < T \leq \alpha_j),
\end{align*}
<!-- -->
where $T \sim \text{Poisson}(e^s \lambda^2 / 4)$.

The following R snippet builds the C++ code and invokes the sampling function.

```{r}
#| prompt: true
Rcpp::sourceCpp("examples/bessel/bessel-v3.cpp")
out3 = r_bessel_v3(n = 1000, lambda, nu, lo = -0.1, hi = 1e5, N = 50,
	tol = 0.10)
head(out3$draws)
```

@fig-bessel-refine plots the bound for the rejection probability for this sampler after each step of refining, along with that of the previous two versions. A substantial improvement in efficiency is seen here; fewer regions are needed to achieve a small rejection probability. Although versions 1 and 2 of the sampler are based on the same proposal, the changes to their bounds are seen to differ due to randomness in knot selection. @fig-bessel-majorizers and @fig-bessel-proposals display the three versions of majorized weight functions and proposal distributions, respectively, after refinement. It is seen that the majorizer from @sec-bessel-const-default does not provide a tight bound, even with many knots, because our numerical optimization is not specific to the integers in the support. This is rectified by the manual optimization from @sec-bessel-const-custom. The linear majorizer provides a tight bound with fewer knots.

Because the mass of our selected $\text{Bessel}(\lambda, \nu)$ distribution is focused on a relatively small set of integers, the second and third versions of the proposal can be made practically equivalent to the target with a sufficient number of regions. After being refined to this point, candidates are rarely discarded in rejection sampling.

```{r}
#| out-width: 60%
#| fig-cap: |
#|   Refinement for Bessel example with three majorizers.
#| label: fig-bessel-refine
#| echo: false

df1 = data.frame(sampler = "v1", bdd = exp(out1$lbdd)) %>% 
	mutate(N = row_number())
df2 = data.frame(sampler = "v2", bdd = exp(out2$lbdd)) %>% 
	mutate(N = row_number())
df3 = data.frame(sampler = "v3", bdd = exp(out3$lbdd)) %>% 
	mutate(N = row_number())
df = bind_rows(df1, df2, df3) %>% mutate(N = as.integer(N)) %>% filter(N <= 18)

ggplot(df) +
	geom_line(aes(x = N, y = bdd, color = sampler), lwd = 1.02, alpha = 1) +
	geom_point(aes(x = N, y = bdd, pch = sampler, color = sampler),
		cex = 3, alpha = 1) +
	scale_x_continuous(breaks = df$N[df$N %% 2 == 0]) +
	scale_y_continuous(n.breaks = 10) +
	xlab("N") +
	ylab("Bound") +
	theme_minimal()
```

```{r}
#| fig-cap: |
#|   Majorizers (on the log-scale) for Bessel example for three sampler
#|   versions. The blue curve is the majorizer with a blue dot marking the
#|   right endpoint of each region. The black triangle displays the value of
#|   $\log w(x)$ at integers $x \in \mathbb{N}$.
#| fig-subcap:
#|   - v1.
#|   - v2.
#|   - v3.
#| label: fig-bessel-majorizers
#| layout-ncol: 3
#| fig-width: 3
#| fig-height: 2.5
#| echo: false

ggplot(out1$df_weight) +
	geom_segment(aes(x = lo, xend = hi, y = lwmax), col = "blue") +
	geom_point(aes(x = hi, y = lwmax), col = "blue") +
	geom_point(data = data.frame(x = xseq, w = lwseq), aes(x, w), pch = 2) +
	coord_cartesian(xlim = c(NA, 15), ylim = c(min(lwseq), NA)) +
	scale_x_continuous(breaks = xseq[xseq %% 2 == 0]) +
	xlab("x") +
	ylab("Logarithm of Weight") +
	theme_minimal()

ggplot(out2$df_weight) +
	geom_segment(aes(x = lo, xend = hi, y = lwmax), col = "blue") +
	geom_point(aes(x = hi, y = lwmax), col = "blue") +
	geom_point(data = data.frame(x = xseq, w = lwseq), aes(x, w), pch = 2) +
	coord_cartesian(xlim = c(NA, 15), ylim = c(min(lwseq), NA)) +
	scale_x_continuous(breaks = xseq[xseq %% 2 == 0]) +
	xlab("x") +
	ylab("Logarithm of Weight") +
	theme_minimal()

out3$df_weight %>%
	mutate(w_lo = beta0_max + beta1_max * lo) %>%
	mutate(w_hi = beta0_max + beta1_max * hi) %>%
	ggplot() +
	geom_segment(aes(x = lo, xend = hi, y = w_lo, yend = w_hi), col = "blue") +
	geom_point(aes(x = hi, y = w_hi), col = "blue") +
	geom_point(data = data.frame(x = xseq, y = lwseq), aes(x, y), pch = 2) +
	coord_cartesian(xlim = c(NA, 15), ylim = c(min(lwseq), NA)) +
	scale_x_continuous(breaks = xseq[xseq %% 2 == 0]) +
	xlab("x") +
	ylab("Logarithm of Weight") +
	theme_minimal()
```

```{r}
#| fig-cap: |
#|   Proposal density (bars) for Bessel example for three sampler versions.
#|   Triangles display the values of $f(x)$, $x \in \mathbb{N}$.
#| fig-subcap:
#|   - v1.
#|   - v2.
#|   - v3.
#| label: fig-bessel-proposals
#| layout-ncol: 3
#| fig-width: 3
#| fig-height: 2
#| echo: false

out1 = r_bessel_v1(0, lambda, nu, N = 50, tol = 0.1, x = xseq)
out2 = r_bessel_v2(0, lambda, nu, N = 50, tol = 0.1, x = xseq)
out3 = r_bessel_v3(0, lambda, nu, lo = -0.1, hi = 1e5, N = 50, tol = 0.1, x = xseq)

data.frame(x = xseq, h = out1$hx) %>%
	ggplot() +
	geom_bar(aes(x, h), stat = "identity", fill = NA, col = "black") +
	geom_point(aes(x, fseq), pch = 2) +
	scale_x_continuous(breaks = xseq[xseq %% 2 == 0]) +
	xlab("x") +
	ylab("Density") +
	theme_minimal()

data.frame(x = xseq, h = out2$hx) %>%
	ggplot() +
	geom_bar(aes(x, h), stat = "identity", fill = NA, col = "black") +
	geom_point(aes(x, fseq), pch = 2) +
	scale_x_continuous(breaks = xseq[xseq %% 2 == 0]) +
	xlab("x") +
	ylab("Density") +
	theme_minimal()

data.frame(x = xseq, h = out3$hx) %>%
	ggplot() +
	geom_bar(aes(x, h), stat = "identity", fill = NA, col = "black") +
	geom_point(aes(x, fseq), pch = 2) +
	scale_x_continuous(breaks = xseq[xseq %% 2 == 0]) +
	xlab("x") +
	ylab("Density") +
	theme_minimal()
```

# Persistent Proposals in R {#sec-persistent}

The samplers developed in previous sections each expose a function in R which can be used to generate samples. A new proposal is constructed each time the function is called and it exists only for the duration of the call. For the sake of efficiency or to maintain the state within a proposal, it may be desirable to have a persistent proposal that can be gradually reused and refined. This section demonstrates several Rcpp mechanisms to enable this workflow: first using external pointers then using using modules. The latter is based on object-oriented programming but is recommended over the former as it is less susceptible to crashing R. We demonstrate the two approaches by revisiting the example from @sec-vmf-const-default. For more information on external pointers and modules, refer to the vignette "Exposing C++ functions and classes with Rcpp modules" in the `Rcpp` package [@Rcpp] and section "External pointers and weak references" of the guide "Writing R Extensions" [@R-ext].

The approach from @sec-persistent-module using Rcpp modules is used in @sec-privacy-within-gibbs to generate from an unfamiliar conditional within a Gibbs sampler.

## External Pointers {#sec-persistent-xptr}

An Rcpp external pointer (`XPtr`) can be used to construct objects in C++ and return them to R so that they persist and can be accessed in subsequent calls to C++. This is especially useful with complicated data structures that are not easy to directly represent as R objects. A VWS proposal is an example of such as structure. The full C++ code is given in `examples/vmf/vmf-v1-xptr.cpp`. Here we describe the components of that file.

The following typedefs provide more readable names for the proposal and its `XPtr`.

```{cpp}
typedef vws::fmm_proposal<double, vws::real_const_region> t_proposal;
typedef Rcpp::XPtr<t_proposal> t_proposal_xptr;
```

The following function creates an `XPtr` to a proposal. There are several notable changes from the implementation in @sec-vmf-const-default.

```{cpp}
// [[Rcpp::export]]
t_proposal_xptr vmf_pre_v1_xptr(double kappa, double d)                  // <1>
{
	vws::dfdb w =
	[=](double x, bool log = true) {                                     // <2>
		double out = R_NegInf;
		if (std::fabs(x) < 1){
			out = 0.5 * (d - 3) * std::log1p(-std::pow(x, 2));
		}
		return log ? out : std::exp(out);
	};

	fntl::density df = [=](double x, bool log = false) {
		return d_texp(x, kappa, -1, 1, log);
	};

	fntl::cdf pf = [=](double q, bool lower = true, bool log = false) {
		return p_texp(q, kappa, -1, 1, lower, log);
	};

	fntl::quantile qf = [=](double p, bool lower = true, bool log = false) {
		return q_texp(p, kappa, -1, 1, lower, log);
	};

	vws::univariate_helper helper(df, pf, qf);
	vws::real_const_region supp(-1, 1, w, helper);

	auto p = new t_proposal(supp);                                       // <3>
	auto out = t_proposal_xptr(p, true);
	set_tag(out, "t_proposal");                                          // <4>
	return out;
}
```
1. The function returns an `XPtr` to a proposal rather than the result of sample generation.
2. Capture by value with the lambdas here using `[=]`. Variables `kappa` and `d` may not be available in the environment when evaluating the lambdas later. Capture by value ensures that variables can still be accessed by the lambdas.
3. Create a proposal object on the heap using `new` and wrap it in an `XPtr` to be returned to R.
4. The function `set_tag` is a utility we will define; more explanation is provided below.

The following functions make use of the proposal `XPtr`. The `refine` function
modifies the proposal by applying up to `N` refinement steps of Algorithm \ref{alg:refine} until toleration `tol` is achieved. The `draw` function draws
a sample of size `n`.

```{cpp}
// [[Rcpp::export]]
Rcpp::NumericVector refine(
	t_proposal_xptr h,                    // <1>
	unsigned int N, double tol = 0)
{
	assert_tag(h, "t_proposal");          // <2>
    return h->refine(N - 1, tol);         // <3>
}

// [[Rcpp::export]]
Rcpp::List draw(t_proposal_xptr h, unsigned int n,
	unsigned int max_rejects = 10000, unsigned int report = 10000)
{
	assert_tag(h, "t_proposal");

    vws::rejection_args args;
    args.max_rejects = max_rejects;
    args.report = report;

    auto out = vws::rejection(*h, n, args);

    return Rcpp::List::create(
        Rcpp::Named("draws") = out.draws,
        Rcpp::Named("rejects") = out.rejects
    );
}
```
1. The first argument is expected to be the `XPtr` obtained from `vmf_pre_v1_xptr`.
2. The function `assert_tag` is a utility we will define; more explanation is provided below.
3. The `refine` member function is called on the object pointed to by `h`; this modifies it in place.

Something to be aware of with `XPtr` arguments is that type safety is not as well-established as with other Rcpp argument types. For example, consider attempting the operation `h->refine(N - 1, tol)` with `h` as an `XPtr` to something other than a `t_proposal`. At the time of this writing, this crashes R rather than the operation being prevented with an error message. To help avoid this issue, we use the tagging mechanism for `XPtr` which is in the Rcpp API. Define the following functions for convenience: `set_tag` applies a string tag to an `XPtr` and `assert_tag` throws an exception if an `XPtr` has no tag or has a tag different from the expected one.

```{cpp}
template <typename T>
void assert_tag(Rcpp::XPtr<T> x, const char* tag)
{
	if (Rf_isNull(R_ExternalPtrTag(x))) {
		Rcpp::stop("XPtr does not have a tag");
	}

	if (Rcpp::as<Rcpp::String>(R_ExternalPtrTag(x)) != Rcpp::String(tag)) {
		Rcpp::stop("XPtr is not tagged as %s", tag);
	}
}

template <typename T>
void set_tag(Rcpp::XPtr<T> x, const char* tag)
{
	R_SetExternalPtrTag(x, Rcpp::wrap(Rcpp::String(tag)));
}
```

Now let us demonstrate use of these functions from R.

```{r}
#| prompt: true
Rcpp::sourceCpp("examples/vmf/vmf-v1-xptr.cpp")
h = vmf_pre_v1_xptr(kappa = 5, d = 4)
lbdd = refine(h, N = 50, tol = 0.10)
out = draw(h, n = 1000, max_rejects = 10000)
rm(h)
print(lbdd)
head(out$draws)
```

For demonstration purposes, let us define a C++ function that returns an `XPtr` to a vector (i.e., something which is not a proposal object). When an `XPtr` to a proposal is expected but an `XPtr` to a vector is encountered, ensure that an appropriate exception is thrown.

```{cpp}
// [[Rcpp::export]]
Rcpp::XPtr<std::vector<double>> test_vector()
{
	std::vector<double>* p = new std::vector<double>();
	p->push_back(10);
	p->push_back(10);
	p->push_back(12);
	return Rcpp::XPtr<std::vector<double>>(p);
}
```

```{r}
#| prompt: false
tryCatch({
	h = test_vector()
	refine(h, N = 50, tol = 0.10)
}, error = function(e) {
	print(e)
})
```

::: {#rem-capture-value}
In the current implementation of `vws`, lambdas `w`, `df`, `pf`, and `qf` are copied to all regions used in the sampler. The overhead created by copying will likely not be significant in most instances, but could be an issue when very large objects are captured by value. In such cases, users may consider using pointers rather than values to capture these objects.
\remarkqed
:::

## Modules {#sec-persistent-module}

The external pointer approach in @sec-persistent-xptr lets us operate on a proposal object from R through C++ functions. We can take this one step further with object-oriented programming so that the proposal appears as a class in R that includes its methods. The complete C++ code for this example is in `examples/vmf/vmf-v1-module.cpp`.

The following typedef provides a more readable name for the proposal.

```{cpp}
typedef vws::fmm_proposal<double, vws::real_const_region> t_proposal;
```

We now define a class named `rcpp_proposal` which will be exposed as an Rcpp module. This class inherits from `t_proposal`. (I.e., `rcpp_proposal` "is a" `t_proposal`.) Its constructor, `refine` operation, and `draw` operation are analogous to those in @sec-persistent-xptr.

```{cpp}
class rcpp_proposal : t_proposal
{
public:
	rcpp_proposal(double kappa, double d)
	: t_proposal(supp(kappa, d))
	{
	}

	Rcpp::List draw(unsigned int n, unsigned int max_rejects = 10000,
		unsigned int report = 10000);

	Rcpp::NumericVector refine(unsigned int N, double tol = 0) {  // <1>
		return t_proposal::refine(N, tol);
	}

private:
	vws::real_const_region supp(double kappa, double d);            // <2>
};
```
1. `t_proposal::refine` is the method in the base class; it is possible to expose this directly from the module, but we opt to call it from `rcpp_proposal::refine` for simplicity.
2. The `supp` function handles most of the constructor work; it returns a single region `supp` with all of the necessary operations. The `rcpp_proposal` constructor creates an instance of itself from region.

The following code exposes `rcpp_proposal` as an Rcpp module.

```{cpp}
RCPP_MODULE(vws_module) {
	Rcpp::class_<rcpp_proposal>("rcpp_proposal")
	.constructor<double,double>()
	.method("draw", &rcpp_proposal::draw)
	.method("refine", &rcpp_proposal::refine)
	;
}
```

See `examples/vmf/vmf-v1-module.cpp` for the implementations for the `draw` and `refine` methods.

Now let us demonstrate use of these functions from R.

```{r}
#| prompt: true
Rcpp::sourceCpp("examples/vmf/vmf-v1-module.cpp")
h = new(rcpp_proposal, kappa = 5, d = 4)
lbdd = h$refine(N = 50 - 1, tol = 0.10)
out = h$draw(n = 1000, max_rejects = 1000, max_rejects = 10000)
print(lbdd)
head(out$draws)
```

::: {#rem-capture-value-2}
The lambdas `w`, `df`, `pf`, and `qf` in this code capture `kappa` and `d` by value as in @sec-persistent-xptr, which can lead to unwanted overhead from copying mentioned in @rem-capture-value. In the present case, we can consider storing a single copy of any large objects as member variables of the `rcpp_proposal` class and capturing them by reference (rather than by value) in the lambdas.
\remarkqed
:::

# Example: Privacy Protected Data and Sampling within Gibbs {#sec-privacy-within-gibbs}

This section presents an example using the self-tuned VWS method described in @sec-self-tune. Consider again the setting from @sec-privacy and suppose there are $n$ tabulations. Let $y_i$ be the $i$th sensitive tabulation, $z_i$ be its privacy-protected counterpart, and $\vec{x}_i \in \mathbb{R}^d$ be a non-sensitive covariate which is associated with $y_i$. Consider the following hierarchical Bayesian model:
<!-- -->
\begin{align*}
&z_i \sim \text{N}(y_i, \lambda_i^2), \\
&y_i \sim \text{Lognormal}(\vec{x}_i^\top \vec{\beta}, \sigma^2), \\
&\vec{\beta} \sim \text{N}(\vec{0}, \sigma_\beta^2 \vec{I}), \\
&\sigma^2 \sim \text{IG}(a_\sigma, b_\sigma),
\end{align*}
<!-- -->
for $i = 1, \ldots, n$, where $\sigma_\beta^2$, $a_\sigma$, and $b_\sigma$ are hyperparameters which are assumed to be known. The variances $\lambda_i^2$ are assumed to be known and provided with the data. Denote $\vec{\Lambda} = \text{Diag}(\lambda_1^2, \ldots, \lambda_n^2)$ and $\log \vec{y} = (\log y_1, \ldots, \log y_n)$. To obtain a Gibbs sampler, the joint distribution of all random variables is
<!-- -->
\begin{align}
f(\vec{z}, \vec{y}, \vec{\beta}, \sigma^2) =
f_\text{N}(\vec{z} \mid \vec{y}, \vec{\Lambda}) \cdot
f_\text{N}(\log \vec{y} \mid \vec{X} \vec{\beta}, \sigma^2 \vec{I}) \cdot
f_\text{N}(\vec{\beta} \mid \vec{0}, \sigma_\beta^2\vec{I}) \cdot
f_\text{IG}(\sigma^2 \mid a_\sigma, b_\sigma) \cdot
\prod_{i=1}^n \frac{1}{y_i},
\label{eqn:privacy-joint}
\end{align}
<!-- -->
where $f$ with subscripts "N", "LN", and "IG" denotes the normal, lognormal, and inverse gamma densities respectively. From \eqref{eqn:privacy-joint}, we obtain the following conditionals after some algebra:

1. $[\vec{\beta} \mid \sigma^2, \vec{y}, \vec{z}] \sim \text{N}(\vec{\vartheta}, \vec{\Omega}^{-1})$ where $\vec{\Omega} = \sigma^{-2}\vec{X}^\top \vec{X} + \sigma_\beta^{-2} \vec{I}$ and $\vec{\vartheta} = \sigma^{-2} \vec{\Omega}^{-1} \vec{X}^\top (\log \vec{y})$.

2. $[\sigma^2 \mid \vec{\beta}, \vec{y}, \vec{z}] \sim \text{IG}(a', b')$ where $a' = a_\sigma + n/2$ and $b' = b_\sigma + \lVert \log \vec{y} -\vec{X} \vec{\beta} \rVert^2 / 2$.

3. $[\vec{y} \mid \vec{\beta}, \sigma^2, \vec{z}] = \prod_{i=1}^n [y_i \mid \vec{\beta}, \sigma^2, \vec{z}]$ where
<!-- -->
\begin{math}
f(y_i \mid \vec{\beta}, \sigma^2, \vec{z}) \propto f_\text{N}(z_i \mid y_i, \lambda_i^2) \cdot
f_\text{LN}(y_i \mid \vec{x}_i^\top \vec{\beta}, \sigma^2).
\end{math}

Steps 1 and 2 can make use of well-known normal and inverse gamma variate generation methods. For step 3, we will use VWS with a constant majorizer and minorizer as in @sec-privacy-const-custom; these can be readily obtained without numerical optimization. 

@sec-privacy-within-gibbs-r presents a version of the code where the Gibbs sampler is written in R with Rcpp modules to interact with VWS proposals written in C++, as in @sec-persistent-module. @sec-privacy-within-gibbs-cpp presents an alternative version where the Gibbs sampler is written in C++ rather than R. The latter is substantially faster because it avoids the overhead of repeated calls between R and C++; however, it requires more C++ programming and cannot be debugged as readily as the former. Several timings are mentioned; they were measured in Windows using a laptop with a Intel Core i7-1185G7 3.0GHz CPU with four CPU cores.

## Gibbs Sampler in R {#sec-privacy-within-gibbs-r}

This version of the code implements the Gibbs sampler in R and utilizes VWS via Rcpp. We consider two implementations of VWS. One uses persistent proposals that are updated with current parameters from the Gibbs sampler. It makes use of the self-tuning method from @sec-self-tune. The other implementation uses the more basic method of constructing new proposals each time they are needed. Code for this example is given in the folder `examples/privacy-gibbs/v1-R`.

The function `draw_ln_norm` in `draw-ln-norm.cpp` draws a variate from the conditional distribution $[y_i \mid \text{---}]$ using the current values of $z_i$, $\lambda_i$, $x_i^\top \vec{\beta}$, and $\sigma^2$. The remaining arguments are specific to VWS: `tol` and `N` correspond to the tolerance and maximum number of knots in Algorithm \ref{alg:refine}; `max_rejects` is the maximum number of rejections permitted in the `rejection` function. The return value is a list with the contents of a `rejection_result<double>`.

```{cpp}
Rcpp::List draw_ln_norm(
	double z,
	double lambda,
	double xbeta,
	double sigma,
	double tol,
	unsigned int N,
	unsigned int max_rejects
);
```

The class `ln_norm_proposal` in `ln-norm-proposal.cpp` is an `fmm_proposal` which defines a few additional functions that are exposed to R and contains state that should persist between calls from R. The interface of the class is as follows.

```{cpp}
class ln_norm_proposal : public vws::fmm_proposal<double, vws::real_const_region>
{
public:
	ln_norm_proposal(double z, double lambda, double xbeta, double sigma); // <1>
	void update(double xbeta, double sigma);                               // <2>
	vws::rejection_result<double> draw(unsigned int max_rejects) const;    // <3>
	vws::rejection_result<double> draw_tune(double tol_suff,
		double tol_merge, unsigned int max_rejects);                       // <4>
	unsigned int size() const;                                             // <5>
};
```
1. Constructor from initial values of $z_i$, $\lambda_i$, $x_i^\top \vec{\beta}$, and $\sigma^2$.
2. Update the proposal with values of $x_i^\top \vec{\beta}$ and $\sigma^2$. These values change at each iteration of the Gibbs sampler. The proposal must update its representation of weight function $w$ with these values and recompute the majorizer and minorizer to be an envelope for rejection sampling.
3. Generate a single draw using the `rejection` function.
4. Generate a single draw and self-tune using the `rejection_tune` function. 
5. Return the number of regions in the proposal.

The following snippet in `ln-norm-proposal.cpp` defines an Rcpp Module which exposes functions to R.

```{cpp}
RCPP_MODULE(vws_module) {
	Rcpp::class_<ln_norm_proposal>("ln_norm_proposal")
	.constructor<double,double,double,double>()
	.method("draw", &ln_norm_proposal::draw)
	.method("draw_tune", &ln_norm_proposal::draw_tune)
	.method("update", &ln_norm_proposal::update)
	.method("size", &ln_norm_proposal::size)
	;
}
```

The file `gibbs.R` contains R code for the Gibbs sampler and some helper functions. See the Roxygen comments that accompany each function for a brief description.

The script `driver.R` generates a dataset from the model and illustrates several calls to the Gibbs sampler. The following snippet uses transient proposals via `draw_ln_norm` (`method = "vws-basic"`), with some details omitted.

```{r}
#| eval: false
inner = control_inner(tol_suff = 0.85, method = "vws-basic")
control = control_gibbs(R = 3000, burn = 1000, report = 100, inner = inner,
	save_latent = 1:10)
out = gibbs(z, lambda, X, control = control)
```

Note that the setting `tol_suff = 0.85` ensures that each proposal is tuned to have a rejection rate no larger than 85% before rejection sampling. The following snippet illustrates self-tuned VWS using `ln_norm_proposal` (`method = "vws-tune"`).

```{r}
#| eval: false
inner = control_inner(tol_suff = 0.85, tol_merge = 0.01, tune = 100,
	method = "vws-tune")
control = control_gibbs(R = 3000, burn = 1000, report = 100, inner = inner,
	save_latent = 1:10)
out = gibbs(z, lambda, X, control = control)
```

The settings in `control_inner` select the self-tuning method with adjustments made during the first 100 iterations of the MCMC, using $\epsilon_1 = 0.85$ and $\epsilon_2 = 0.01$. See the file `driver.R` script for details.

A particular run of both methods yielded elapsed times of 30 seconds for `vws-basic` and 25 seconds for `vws-tune`. Much of the gain in efficiency due to self-tuning is negated by the overhead from frequent calls between R and C++. This could be mitigated by passing the list of $n$ proposals to C++ at each iteration of the Gibbs sampler. The following section avoids the overhead further by implementing the Gibbs sampler in C++.

## Gibbs Sampler in C++ {#sec-privacy-within-gibbs-cpp}

This version of the code implements the Gibbs sampler in C++ so that it is not necessary to expose VWS proposals to R. Code for this example is given in the folder `examples/privacy-gibbs/v2-cpp`.

The class `ln_norm_proposal` in `ln-norm-proposal.h` defines a VWS proposal as follows.

```{cpp}
class ln_norm_proposal : public vws::fmm_proposal<double, vws::real_const_region>
{
public:
	ln_norm_proposal(double z, double lambda, double xbeta, double sigma); // <1>
	void update(double xbeta, double sigma);                               // <2>
};
```
1. Constructor from initial values of $z_i$, $\lambda_i$, $x_i^\top \vec{\beta}$, and $\sigma^2$.
2. Update the proposal with values of $x_i^\top \vec{\beta}$ and $\sigma^2$. These values change at each iteration of the Gibbs sampler. The proposal must update its representation of weight function $w$ with these values and recompute the majorizer and minorizer to be an envelope for rejection sampling.

This version of `ln_norm_proposal` is simpler than the one in @sec-privacy-within-gibbs-r. When programming in C++, we can directly make use of the fact that it inherits from `fmm_proposal` as a base class. For example, we can call the `fmm_proposal::r` method to generate draws. We can also pass objects of class `ln_norm_proposal` to the `rejection` or `rejection_tune` functions.

The files `gibbs.cpp` and `gibbs.R` contain code for the Gibbs sampler and helper functions: `gibbs.cpp` implements the Gibbs sampler itself in C++ using RcppArmadillo [@RcppArmadillo] for matrix algebra, while helper functions are implemented in R in `gibbs.R`. Refer to the Roxygen comments in `gibbs.R` for details.

The script `driver.R` generates a dataset from the model and runs two versions of the Gibbs sampler, `"vws-basic"` and `"vws-tune"`, as in @sec-privacy-within-gibbs-r. Here, a run of each yielded elapsed times of 28 seconds for `vws-basic` and 6.2 seconds for `vws-tune`. This represents a reduction of $(28 - 6.2) / 28 = 77.9\%$ compared to $(30 - 25) / 30 = 16.67\%$ seen in @sec-privacy-within-gibbs-r.

\clearpage

# References

