Generalized BART with bartisan

A Sum of Trees

A binary regression tree is a pair \((T, M)\), where \(T\) is the tree structure, meaning the splitting rules at the interior nodes, and \(M = \{\mu_1, \dots, \mu_B\}\) collects one parameter for each of its \(B\) leaves. The tree defines a function \(g(x; T, M)\) that sends \(x\) down the tree according to the rules and returns the parameter in the leaf it reaches. A single tree is a step function.

BART sums many of them (Chipman et al. 2010):

\[f(x) = \sum_{m=1}^{M} g(x; T_m, M_m).\]

Each tree is kept small by its prior, so no single one explains much. The sum is flexible while each term is a weak learner, which makes the fit stable. Interactions come free: a path through a tree that splits on \(x_1\) and then on \(x_2\) is an interaction between them, and nothing had to be specified for it to appear.

In ordinary BART, the sum is the conditional mean of a Gaussian outcome,

\[Y = f(x) + \varepsilon, \qquad \varepsilon \sim \mathrm{N}(0, \sigma^2).\]

The Generalization

bartisan keeps the sum of trees and replaces what sits on top of it. The forest supplies an additive predictor \(\eta(x)\), and the outcome follows whatever density the family names,

\[\eta(x) = \sum_{m=1}^{M} g(x; T_m, M_m), \qquad Y \mid x \sim p\{y \mid \eta(x), \theta\},\]

with \(\theta\) collecting any parameters of the density that are not functions of \(x\), such as a residual standard deviation or a shape. These are the nuisance parameters, and they are drawn alongside the trees.

Some families need more than one function of \(x\). A location-scale model has a mean and a log standard deviation; a zero-inflated count has a rate and a probability of a structural zero. Those get one forest each:

\[\eta_h(x) = \sum_{m=1}^{M_h} g(x; T_{hm}, M_{hm}), \qquad h = 1, \dots, H,\]

and the sampler cycles over all the forests. vignette("families") says how many forests each family uses.

The Prior

The model prior has three pieces, all with defaults chosen so that the fit is regularized without being tuned (Chipman et al. 2010).

Tree structure. A node at depth \(d\) is split rather than left as a leaf with probability

\[\gamma (1 + d)^{-\beta}, \qquad \gamma = 0.95, \; \beta = 2,\]

which are the gamma and beta arguments of bartisan_control().

This falls off quickly, so trees are shallow: most have two or three leaves. It is the main reason a forest of 50 trees does not overfit.

Leaf values. Each \(\mu_b\) is normal with mean zero and a standard deviation that shrinks as trees are added, so that the sum has a sensible scale whatever \(M\) is. The k argument in bartisan_control() sets how hard this shrinks: larger values of k pull the fit toward the intercept-only model. The default is a compromise that works across a wide range of problems, and it is rarely worth changing.

Nuisance parameters. Each family sets its own, for example a half-Cauchy on a scale parameter.

The prior does the work that cross-validation does in other tree ensembles. There is no learning rate, no depth limit, and no early stopping.

Soft Rules (gate)

A rule in classic BART is a step: an observation with \(x_j \le c\) goes left and everything else goes right. The fitted function is then a step function, which is a poor description of a smooth truth.

Soft rules replace the step with a gate (i.e., a smooth function of the distance from the cutpoint) (Linero and Yang 2018). An observation goes left with weight

\[\psi\left(\frac{x_j - c}{\tau}\right),\]

where \(\psi\) is a cumulative distribution function and \(\tau\) is a bandwidth, and right with the remaining weight. An observation near a cutpoint reaches both of that node’s children, so the fitted function is smooth rather than a step. As \(\tau \to 0\) the gate becomes a step and the hard rule is recovered.

The gate argument chooses \(\psi\): "smoothstep" (the default), "smootherstep", "logistic", or "hard". They differ in how far the smoothing reaches. The logistic has unbounded support, so under it every observation really does reach every leaf with some weight. The two polynomial gates are the cumulative distribution functions of symmetric Beta kernels and reach exactly zero and one outside a finite window, so an observation far enough from a cutpoint takes one side only and the subtree on the other side is never visited. That skipped work turns out to be worth less than it sounds; they are faster because a polynomial needs no exp(), and they keep that advantage even at a bandwidth wide enough to truncate nothing at all. The bandwidth argument sets the prior mean of \(\tau\), which is drawn rather than fixed, and it means the same amount of smoothing under every gate, the kernels being matched on their standard deviation rather than on their width.

Soft rules cost about three times as much per iteration but are usually worth it:

friedman <- function(n) {
  x <- as.data.frame(matrix(runif(n * 10), n, 10))
  names(x) <- paste0("x", 1:10)
  x$eta <- 10 * sin(pi * x$x1 * x$x2) + 20 * (x$x3 - 0.5)^2 +
    10 * x$x4 + 5 * x$x5
  x
}

train <- friedman(n)
test  <- friedman(1000)
train$y <- train$eta + rnorm(n)

timed <- function(gate) {
  started <- proc.time()[["elapsed"]]
  fit <- bartisan(y ~ . - eta, data = train, family = gaussian(),
                  control = bartisan_control(num_trees = 20, num_burn = 300,
                                             num_draws = 300, gate = gate))
  data.frame(rules = gate,
             test_rmse = sqrt(mean((predict(fit, newdata = test) - test$eta)^2)),
             seconds = round(proc.time()[["elapsed"]] - started, 1))
}

rbind(timed("smoothstep"), timed("hard"))
#>        rules test_rmse seconds
#> 1 smoothstep     0.401     0.3
#> 2       hard     1.169     0.2

The true function has a standard deviation of about 4.9, so both are fitting real structure and the soft fit is the more accurate.

The Predictor Transform (x_transform)

Numeric predictors are mapped to \([0, 1]\) before any rule sees them, and the map is not cosmetic. A cutpoint is drawn uniformly on a node’s live range in the transformed coordinate, so the map decides where a rule can split, and under soft rules it is also the scale on which \(\tau\) is measured. The x_transform argument of bartisan_control() chooses it. The three options fail in different places, which is the whole of the choice between them.

"range" rescales linearly, \(T(x) = (x - \min x) / (\max x - \min x)\), and so preserves the original spacing. That makes it the most accurate of the three on a well-behaved predictor, and, being the only affine map among them, the only one a derivative passes through cleanly. Its failure is a predictor whose bulk sits inside a sliver of its range, which a few extreme values are enough to produce: a cutpoint drawn uniformly on the range then almost never lands where the structure is. In a simulation with 1% of a predictor at 300 times the scale of the rest and a truth that oscillates within the bulk, it was 4 times worse than the alternatives at \(n = 500\) and 6 times worse at \(n = 2000\), with 95% intervals covering .24 rather than .95. The deterioration with \(n\) is the signature: the extremes grow with the sample, so the bulk occupies an ever smaller share of the range.

"quantile" maps each predictor through its own empirical distribution function, as SoftBart::softbart() does. It cannot fail that way, because it places cutpoints by rank and knows nothing about spacing. It fails instead by being a step function, so the fit is a step function of the original predictor: there is no derivative to take, and a relationship that is straight in the predictor becomes a staircase wherever the data have gaps. On two tight clusters with a linear truth it was nearly three times worse than the alternatives.

"smoothcdf", the default, smooths that same distribution function and so does neither. The estimate is

\[\hat F_h(x) = \frac{1}{n} \sum_{i=1}^n K\left(\frac{x - x_i}{h}\right),\]

where \(K\) is the integrated Epanechnikov kernel and \(h = 2.1 \hat\sigma n^{-1/3}\), with \(\hat\sigma\) the smaller of the standard deviation and the interquartile range over 1.349 so that a long tail sets the bandwidth from the bulk rather than from itself. The \(n^{-1/3}\) rate is the one for a distribution function rather than the \(n^{-1/5}\) that is right for a density (Tenreiro 2006), a rate Azzalini (1981) had pointed out for a second-order approximation. The result is stretched affinely so that the smallest and largest observations land on 0 and 1, as they do under the other two maps; normalizing on the smoothed values instead would leave a few percent of the coordinate outside the data, where no cutpoint can be useful. Cutpoints then land where the data are, as under "quantile", and the map is strictly increasing and differentiable, as under "range". Across eight data-generating processes it was never the most accurate and never far from it, where each of the other two was badly wrong somewhere.

It is not a free lunch for slopes. Writing the fit as \(f(T(x))\), a derivative in the original coordinate is \(f'(T(x)) T'(x)\): an affine \(T\) has a known constant derivative, so only \(f'\) is estimated, where a smoothed distribution function contributes an estimated density and the two errors multiply. vignette("effects") measures what that costs and recommends refitting with "range" when a slope is the quantity being reported.

The kernel is Epanechnikov rather than Gaussian because its support is bounded, which makes the estimate cheap to compute exactly. Written directly, \(\hat F_h\) on a grid of \(m\) points is \(nm\) kernel evaluations, tens of millions per predictor on a large fit. Almost all of that work is avoidable, because almost every weight is a 0 or a 1 rather than something in between: an observation further below a grid point than the kernel reaches has passed it entirely and contributes 1, one further above contributes 0, and only those within \(h\) need evaluating. With the data and the grid both sorted, the two ends of that window advance monotonically, so each is found by a pointer that never goes back and the sweep costs \(O(n + m)\).

Sparsity (sparsity)

By default (sparsity = FALSE), every predictor has the same probability of being the variable a rule splits on, which is classic BART. Setting sparsity = TRUE draws that variable from a categorical distribution whose probabilities have a Dirichlet prior instead (Linero 2018):

\[s \sim \mathrm{Dirichlet}(\alpha/p, \dots, \alpha/p),\]

with \(\alpha\) itself drawn. Small \(\alpha\) concentrates the probability on a few predictors, so the forest can stop splitting on the rest entirely. This is DART.

The difference matters for variable selection. Chipman et al. (2010) note that counting splits works poorly when there are many trees, “because the redundancy offered by so many trees tends to mix many irrelevant predictors in with the relevant ones”, and recommend reducing the number of trees so that predictors compete. The Dirichlet prior addresses the same problem directly, by letting irrelevant predictors be dropped rather than crowded out.

Both effects are visible on the Friedman function, where x1 to x5 matter and x6 to x10 are noise. The table reports the average posterior probability of being used, over three replicates at \(n = 500\):

Under the default, the advice to use fewer trees is essential if the table is to be read for selection: at 50 trees the noise predictors are used in 95% of draws and are indistinguishable from the real ones. With the sparsity prior on, the noise predictors stay near zero at every tree count, and reducing the trees buys little. The recommendation is a workaround for the absence of the prior rather than a property of variable selection in general. However, sparsity needs signal, so on a small sample with a weak fit, nothing separates, and the honest reading of a flat table is that the data cannot rule anything out; and it says nothing about effect size, since a predictor can be split on constantly and move the prediction very little. vignette("importance") covers both.

The prior is off by default because it is a poor choice when an effect is the estimand: a predictor it drops from the forest has a contrast of exactly zero in that draw, so the posterior of the effect carries a point mass at zero and a weak effect is attenuated. ?bartisan_control explains the trade-off and vignette("effects") covers what it does to an estimate.

Bayesian Backfitting

The trees are updated one at a time. To update tree \(m\), subtract the contribution of every other tree from the predictor, leaving a partial residual that only tree \(m\) has to explain, then draw \((T_m, M_m)\) given that residual. Cycling over all trees, and then over the nuisance parameters, is one sweep. The draws saved after warmup are the posterior sample.

Moving the Trees

A tree is proposed by a small local change: growing a leaf into two, pruning a pair of leaves back into one, or changing a splitting rule. The proposal is accepted or rejected by a Metropolis-Hastings step. Because growing and pruning change the number of parameters, this is a reversible-jump sampler.

The acceptance ratio needs the likelihood of the tree with the leaf parameters integrated out,

\[p(\text{residual} \mid T) = \int p(\text{residual} \mid T, M) \, p(M) \, dM\]

The General Case, and the Laplace Approximation

For a Gaussian outcome with a normal leaf prior, that integral is a normal integral and has a closed form. This confines classic BART to Gaussian likelihoods, or to likelihoods that can be made Gaussian by augmentation.

For anything else, the integral has no closed form. Linero (2025) replaces it with a Laplace approximation: expand the log integrand around its mode, integrate the resulting quadratic, and use that both to score the tree and to propose the leaf values. A Metropolis correction then makes the chain target the exact posterior. The approximation only has to be good enough to be accepted often; it does not have to be right.

The practical consequence is the reason the package exists. A family needs only the log density of one observation and its first two derivatives with respect to the predictor. Anything that can supply those can be fit, which is why custom_family() takes an R function and works.

The sampler recognizes three shapes of the leaf-level target and uses the cheapest machinery each admits:

A family that lands higher in this table is several times faster than one that does not. Data augmentation exists to move a family higher in this table.

Data Augmentation (augment)

Some likelihoods move up the table once an extra latent variable is imagined and drawn alongside everything else. The variable is not part of the model; it is machinery, and the posterior over everything else is unchanged.

Several such rewritings are in use here, and they do not all arrive at the same shape. A probit model is a Gaussian one on a latent scale truncated by the observed category (Albert and Chib 1993), so drawing the latent value makes the target exactly quadratic. A logistic model is Gaussian conditional on a Pólya-Gamma variable (Polson et al. 2013), which acts as a per-observation precision, and an ordinal complementary log-log model is Gaussian conditional on an exponential waiting time. A censored outcome becomes an uncensored one once the unobserved value is imputed from its truncated distribution, which for the log-normal and log-logistic survival families leaves a quadratic target. A negative binomial is a Poisson whose rate carries a gamma prior, so drawing the rate leaves the exponential shape rather than the quadratic one, which is a smaller gain. And a zero-inflated count is blocked by a mixture rather than by a link, so what is drawn there is the indicator of which component produced each observation; that separates the two forests, and for a zero-inflated negative binomial a second rewriting goes on top of the first.

augment = TRUE, which is the default wherever a rewriting exists, turns these on. They change the sampler and not the model: the posterior is the same, reached faster. vignette("families") says which families have one, and the section on speed below measures what each is worth.

Nuisance Parameters

Anything in the density that is not a function of \(x\) is drawn once per sweep, after the trees. A residual standard deviation, a negative binomial dispersion, the cutpoints of an ordinal model, and the baseline hazard of a proportional hazards model are all handled this way. They come back in fit$aux and are summarized by summary().

A likelihood written with custom_family() can have them too, by naming them in aux_names. There is no prior argument: a parameter with a restricted range is handled by writing the transform into the density, exactly as it would be for a real predictor.

The Dirichlet Process Mixture

Two families, dpm() and dpm_aft(), do not assume a shape for the error at all. They model it as a mixture of normals with a Dirichlet process prior on the mixing distribution (Escobar and West 1995; George et al. 2019), which lets the number of components grow with the data. The mixture is constrained to have mean zero so that the forest still carries the conditional mean.

Conditional on which component each observation belongs to, the model is Gaussian and the target is quadratic, so the extra flexibility is cheap. error_density() reports the fitted error distribution, which is the diagnostic the family exists to provide. Because the mixture can collapse to a single component, it costs little when a single normal was right.

Missing Predictor Values

Missing predictor values are handled natively, by treating missingness as something to split on (Twala et al. 2008; for BART, Kapelner and Bleich 2015). A rule can send missing values left, send them right, or split on missingness itself, and which of these is used is part of the posterior. No imputation is required and rows are not dropped.

d_miss <- train
d_miss$x1[sample(n, 60)] <- NA

fit_miss <- bartisan(y ~ . - eta, data = d_miss, family = gaussian(),
                     control = ctrl)

sqrt(mean((predict(fit_miss, newdata = test) - test$eta)^2))
#> [1] 0.61

The fit uses the incomplete variable rather than discarding it. Note that this is a prediction tool and not a substitute for thinking about why values are missing: what it estimates is the mean of the response given the predictors and the pattern of missingness, which is the quantity prediction calls for, and which is not the regression or causal effect defined on complete data that multiple imputation targets.

Several Chains (chains)

The chains argument runs the sampler more than once from different starting points. The default is one, which is enough for estimates. Running several makes a between-chain comparison possible, which is the only way to see whether the sampler has converged, and diagnose() is where that comparison is computed. It is not part of the fit: the per-observation statistics cost more than the sampling does, so they run when they are asked for.

fit_chains <- bartisan(y ~ . - eta, data = train, family = gaussian(),
                       chains = 4, control = ctrl)

diagnose(fit_chains)$table
#>                              quantity rhat rhat_late ess_bulk ess_tail ess_frac
#> 1                              loglik 1.71     1.974     6.54     19.7  0.00545
#> 2                           aux.sigma 1.12     1.158    24.11    459.6  0.02009
#> 3                          splits.eta 2.15     2.420     5.27     23.7  0.00439
#> 4 eta.eta (average over observations) 1.00     0.997  1183.17   1090.8  0.98598
#> 5  eta.eta (worst 5% of observations) 1.62     1.714     6.78     20.3  0.00565
#> 6      bandwidth (average over trees) 1.23     1.310    12.90     39.7  0.01075
#> 7       bandwidth (worst 5% of trees) 2.03     2.666     5.50     12.9  0.00458
#>   rhat_bad late_bad
#> 1        1        1
#> 2        1        1
#> 3        1        1
#> 4        0        0
#> 5        1        1
#> 6        1        1
#> 7        1        1

This table is read selectively. aux.sigma and loglik are close to 1, as they are at convergence, and they respond to longer chains in the usual way.

The eta row does not, and it is worth knowing why before it causes alarm. It summarizes the worst 5% of the observations rather than the single worst, so it is a high percentile by construction, and forests mix slowly on their fitted values: two chains can visit quite different sets of trees that imply very similar predictions, which inflates a between-chain statistic without meaning that the answers disagree. Adding draws does not reliably bring it down. Values above 1.05 on this row are ordinary for a BART fit, and are characteristic of the method rather than of this implementation.

What to check instead is the quantity that will be reported, on its own draws. diagnose() takes the output of estimate_effect() for exactly this, and posterior::summarise_draws() takes anything else that can be arranged as draws by chains. A high eta row is not on its own a reason to discard a fit, but neither is a low one a reason to trust an effect read off it: a contrast can mix badly where the function it contrasts mixes well, so it has to be diagnosed rather than inferred.

Chains run sequentially unless a future plan is set, in which case they run in parallel:

future::plan("multisession", workers = 4)

fit <- bartisan(y ~ ., data = d, chains = 4)
diagnose(fit)

The chain is the only parallel axis the sampler has, since a sweep conditions on the one before it, and it also makes a convergence diagnostic possible. Outside the sampler, a plan is also used when one is set: the convergence pass splits its per-observation columns across workers, and estimate_effect() predicts the treatment levels in parallel. Any future backend works, including mirai’s, and one set.seed() reproduces the whole run whatever the backend.

diagnose() is the one call for whether the fit can be reported. It computes split-R-hat and the bulk and tail effective sample sizes for every scalar the sampler draws, for the fitted function twice over, once averaged across observations and once over its worst 5% of them, and for the size of the forest itself, then says which of them fall short and what to change. Handed the output of estimate_effect() instead of a fit, it reports on the contrast and on the two potential outcomes it is a contrast of. Some of its choices are worth knowing about. It repeats R-hat on the second half of the draws alone, which distinguishes a warmup that ended too early from chains that have each settled somewhere different: discarding the early draws has the same effect as a larger num_burn, so if that fixes R-hat then warmup was the problem, and the advice says so rather than listing everything a reader might try. And it leaves the leaf scale out, because that parameter mixes badly in every implementation (1.12 in dbarts and 1.16 in stochtree against 1.19 here, on the same data) and its disagreement never reaches the fitted function.

A within-chain drift test would answer the warmup question more directly, and is not offered, because it cannot be made to work at a forest’s autocorrelation. Three versions were calibrated against stationary series, where there is nothing to find by construction: taking each half’s Monte Carlo error from that half alone fires 31% of the time at an autocorrelation of .995 against a nominal 5%, taking it from the whole chain holds specificity under 3% but then misses a six-standard-deviation trend three times in four, and batch means catch everything and fire 92% of the time on a chain that has converged. Re-running a statistic that is already calibrated sidesteps the choice.

vignette("diagnostics") covers all of this from the reader’s side, including which statistic to look at when the fitted values are the estimand.

Choosing the Settings

Most of bartisan_control() does not need to be touched. In rough order of how often changing one is worth it:

Everything else exists so that the checks in the test suite can be written, or because a specific family needs it.

The single decision that matters more than all of these is the family. vignette("families") covers it.

Gains From Data Augmentation (augment)

A major contributor to bartisan’s speed is data augmentation. Where the rewriting reaches a Gaussian conditional the Laplace approximation stops being an approximation at all, which is the largest of the gains: the augmentation uses Albert and Chib (1993) for a probit link and for an ordinal probit, Pólya-Gamma augmentation (Polson et al. 2013) for a logit link, an ordinal logit, or a multinomial, and imputation of the censored failure times for the log-normal and log-logistic survival models. The negative binomial and the zero-inflated families are rewritten too, but to a Poisson rather than to a Gaussian, and the two paragraphs after the table say what that costs them.

Each of them is faster, from 1.2 times for the negative binomial under soft rules to 31 for an ordinal probit with hard rules, with a median of 8.7. What each costs in mixing is the question, and the answer is not that they all cost something: over 15 replicate datasets the worst-mixing quantity’s effective sample size came out a median of 0.89 times what the direct likelihood gives, and of the seventeen cases below, seven have an 80% interval lying entirely below 1, three lie entirely above it, and seven include it. Effective draws per second, which is the ratio to judge, favors augmentation in fifteen of the seventeen, the two negative-binomial rows being the exceptions:

Every cell is the median over replicates of the ratio within a replicate, which is the reading the design supports: both arms see the same data, so the paired ratio is the quantity with the least noise in it, and a median rather than a mean because a wall-clock timing can be disturbed by work outside the fit. The bracketed range is an 80% bootstrap interval over the replicates, and it is in the table because without it the column invites conclusions it cannot support: at three replicates the same quantity varied by a median factor of 4.5 within a single cell, since the worst-quantity effective sample size is a minimum over many quantities and so is a high-variance thing to estimate.

The augment argument is on by default for every family that has a rewriting: the binomial, ordinal, multinomial, negative binomial, zero-inflated, and survival families. The negative binomial’s is a Poisson whose rate comes from a gamma, which needs one gamma draw per observation rather than a Pólya-Gamma one, and it is the one family whose two rows differ in kind rather than degree. Under hard rules the target reaches the exponential form and the rewriting buys time at a real cost in mixing; under soft rules it cannot, and the rewriting instead buys a little time at no measurable cost in mixing. Neither comes out ahead on effective draws per second, at 1.0 and 0.8 times, which makes the negative binomial the one family whose rewriting does not pay for itself; setting augment = FALSE is worth trying whenever a negative binomial fit’s diagnostics look poor.

The survival families are the clearest case of what the rewriting is for. Right-censoring makes their likelihood expensive: an observed failure contributes a density in the additive predictor and a censored one contributes a survival function, and the two have different shapes, so nothing about the target is exploitable. Imputing each censored failure time above its censoring time replaces the survival term with a density, and then every observation contributes the same quadratic shape. weibull_aft() is absent from the table because it needs none of this, its likelihood being already exponential in the sense above, censoring included, which also makes it the slowest of the three at the default soft gate, where that form cannot be used.

The zero-inflated families are the odd ones out, because what blocks them is a mixture rather than a link: the zero contributes \(\log[\pi + (1 - \pi) P_0]\), so neither the count predictor nor the inflation predictor has a shape to exploit. Introducing the indicator of which component produced each observation separates them, and the count forest then sees a plain Poisson or negative binomial, and the inflation forest a Bernoulli logistic likelihood that Pólya-Gamma handles.

The ordinal probit is the largest of these because the target it replaces is the most expensive one in the package: two cumulative-normal evaluations per observation per pass, and no exploitable shape, so every trial value of a leaf parameter costs its own pass. It comes with one subtlety worth naming. The cutpoints are not drawn from their conditional given the latent normals, which would be uniform between the two order statistics that bracket each one, pinning them to an interval of width \(O(1/n)\) and mixing worse the more data there is. They are drawn from the ordinal likelihood with the latent normals integrated out, and the latent normals are redrawn immediately afterward: a partially collapsed Gibbs sampler (van Dyk and Park 2008), which is the standard remedy (Cowles 1996).

The ordinal logit has no latent normal, its latent variable being logistic, and the usual route to one goes through the Kolmogorov-Smirnov density, which needs a sampler of its own. However, Polson et al. (2013)’s Theorem 1 at \(a = 1\), \(b = 2\) says that the standard logistic density is \(\frac{1}{4} E[\exp(-w x^2 / 2)]\) with \(w \sim \mathrm{PG}(2, 0)\), so a logistic residual is a normal whose precision is Pólya-Gamma; and because \(\mathrm{PG}(b, c)\) is \(\mathrm{PG}(b, 0)\) tilted by \(\exp(-c^2 w / 2)\), the precision’s conditional given a residual \(r\) is exactly \(\mathrm{PG}(2, \lvert r \rvert)\). That is an integer-parameter draw, which the exact Devroye sampler already in the package covers, so the logit link gets the same quadratic target with nothing approximate anywhere. The Poisson and gamma families need no rewriting, their targets being already in a shape the sampler collapses to a single pass.