Pith. sign in

REVIEW 3 major objections 5 minor 29 references

The Lasso Distribution: Properties, Sampling Methods, and Applications in Bayesian Lasso Regression

T0 review · 3 major / 5 minor · reviewed 2026-08-07 · deepseek-v4-flash

Challenge this review Re-run · record.json Download PDF Read on arXiv ↗

Pith's one-line read The Bayesian lasso's per-coefficient full conditional is a new exponential-family distribution with closed-form moments and a numerically stable inverse-CDF sampler.

desk verdict Careful and mostly correct distributional work, but the inverse-CDF sampler's branch rule is printed backwards; fix that and this is a solid computational statistics paper. read the letter →

arxiv 2506.07394 v3 pith:6JSSN5RE submitted 2025-06-09 stat.CO

Mohammad Javad Davoudabadi , Jonathon Tidswell , Samuel Muller , Garth Tarr , John T. Ormerod This is my paper
classification stat.CO MSC 62E1562F1562J0765C10
keywords LassodistributionBayesianGibbssamplertruncatednormalmixtureinverse-CDFsamplingMillsratioexponentialfamilyhigh-dimensionalregression
verification ladder T0 review T1 audit T2 compute T3 formal

The pith

A machine-rendered reading of the paper's core claim, the machinery that carries it, and where it could break.

The reading

The paper introduces the Lasso distribution, whose density is proportional to $\exp(-\tfrac{1}{2} a x^2 + b x - c|x|)$, and argues that this is exactly the full conditional distribution of each regression coefficient in a Bayesian lasso Gibbs sampler. It derives closed-form expressions for the normalizing constant, moments, moment-generating function, and cumulative distribution function, and shows the distribution belongs to the exponential family. The paper also provides a numerically stable inverse-CDF sampler, built on a rational approximation of Mill's ratio, so that each coefficient can be drawn directly instead of as a weighted mixture of two truncated normals. The authors report that the resulting sampler mixes well and runs faster than their comparison implementations on benchmark regressions.

What carries the argument

The load-bearing object is the unnormalized kernel $\exp(-\tfrac{1}{2} a x^2 + b x - c|x|)$ with parameters $a\geq 0$, $b\in\mathbb{R}$, $c\geq 0$. The paper normalizes it by expressing the integral as a sum of two Mill's ratios $m(x)=\bar{\Phi}(x)/\phi(x)$, evaluates $m$ with a Remez-optimized degree-(8,9) rational approximation, and represents the density as a mixture of two truncated normals (one positive, one negative). That mixture representation gives the mixing weight $w$, the CDF, the inverse CDF, and closed-form moments; the inverse CDF is then the basis of the direct sampling algorithm used inside the Gibbs update.

What would settle it

For a concrete counter-check, take $(a,b,c)=(2,1,3)$, for which the closed-form mixing weight is $w\approx 0.374$; generate $10^6$ draws from the package's sampler and estimate $w$ as the fraction of draws below the mode. If the estimated weight disagrees with $0.374$ at three decimal places, the inverse-CDF branch rule is not what the paper's formulas imply.

Watch

Extended reading notes

Core claim

The central claim is that the kernel $p(\beta_j|\cdot)\propto \exp(-\tfrac{1}{2} a\beta_j^2 + b\beta_j - c|\beta_j|)$ that arises in the Bayesian lasso is not an awkward special case to be sampled by combining two truncated normals, but a member of a well-behaved parametric family: the Lasso distribution $\mathrm{Lasso}(a,b,c)$. The paper derives the normalizing constant in terms of Mill's ratio, writes the density as a mixture of a positively and a negatively truncated normal with mixing weight $w$, and obtains closed-form moments via the MGF $M(t)=Z(a,b+t,c)/Z(a,b,c)$. It then gives a four-case inverse-CDF sampler with careful branch selection to avoid overflow and underflow, and embeds this sampler in a modified Gibbs algorithm for the Bayesian lasso. If correct, every coefficient update in the Gibbs sampler becomes one direct draw from a known distribution, replacing a more fragile two-component truncated-normal step.

Load-bearing premise

The sampler's correctness rests on the branch-selection rule that chooses which algebraic form of the mixing weight $w$ (or $1-w$) to use based on the sign of $b$, together with the accuracy of the rational approximation to Mill's ratio over the whole parameter range.

Editorial extensions

If this is right

  • Each coefficient's full conditional in the Bayesian lasso can be sampled in one inverse-CDF draw, removing the need to simulate two truncated normals and combine them by hand.
  • Because the Lasso distribution is in the exponential family with sufficient statistic $(X^2, X, |X|)$, posterior analysis for this kernel inherits standard exponential-family structure.
  • The sampler runs in $O(N p \min(n,p))$ time when the design-matrix products are precomputed, which makes high-dimensional settings (large $p$, or $n$ much smaller than $p$) more tractable than samplers that require matrix square roots.
  • The same distributional machinery applies anywhere the kernel $\exp(-a x^2/2 + b x - c|x|)$ appears, beyond lasso regression.

Reading between the lines

Editorial extensions of the paper, not claims the author makes directly.

  • The mixture-of-truncated-normals representation suggests a direct two-stage sampler (choose component with probability $w$, then draw from a truncated normal), which could serve as a cross-check against the inverse-CDF code and might be more efficient for moderate parameters.
  • Closed-form moments make method-of-moments estimation of $(a,b,c)$ from sample data a natural next step, potentially giving a fast estimator for the shrinkage and asymmetry parameters without MCMC.
  • The exponential-family form invites generalized linear model extensions where the same kernel appears as a conditional prior or penalty, e.g., penalized regressions with asymmetric $\ell_1$ penalties.
  • A direct test of the branch-selection rule in the inverse-CDF sampler against the closed-form mixing weight on a grid of parameter values would determine whether the implemented R package matches the paper's stated formulas.
Share X Bluesky LinkedIn Reddit HN

Editorial analysis

A structured set of objections, weighed in public.

Desk editor's note, referee report, and a circularity audit.

Referee Report

3 major / 5 minor

Summary. The paper introduces a univariate probability distribution, the Lasso distribution, with density proportional to exp(-a x^2/2 + b x - c|x|), and derives its normalizing constant, CDF, inverse CDF, moments, MGF, and exponential-family representation. It then proposes an inverse-CDF sampling algorithm based on a Remez rational approximation to Mill's ratio, and presents two Gibbs samplers for Bayesian lasso regression: a modified Hans sampler that draws each regression coefficient from the Lasso distribution, and a modified Park--Casella sampler. The methods are implemented in an R package (BayesianLasso) and compared with existing samplers on three benchmark datasets.

Significance. The distributional contribution is genuinely useful: the full conditional for each coefficient in the Hans (2009) Gibbs sampler is exactly of this form, and naming it, deriving its properties, and providing a stable sampler is a practical advance for Bayesian lasso computation. The derivations of the PDF, CDF, inverse CDF, moments, and MGF are self-contained and, on inspection, algebraically consistent. The paper also ships a concrete, CRAN-available implementation, which is a real strength: the reader can reproduce the distributional computations and check the samplers. The performance comparison is less rigorous than the distributional theory: it is based on single runs with no reported uncertainty, and the algorithmic description contains a sign error that, if implemented literally, would break the sampler. The central theoretical claims are sound, but the practical sampling claim needs correction and verification.

major comments (3)
  1. [§3.1] The rule for choosing between w1 and w2 is reversed. In §3, w is defined as the weight of the negative-support truncated-normal component B ~ TN^-(μ2, σ²) in the mixture representation E(X^r) = (1-w)E(A^r) + wE(B^r). With v1 = (c-|b|)/√a and v2 = (c+|b|)/√a, the negative-component weight is m(v2)/(m(v1)+m(v2)) when b>0, which equals w1, not w2, as defined in §3.1. The text states w = w1 if b ≤ 0 and w = w2 if b > 0, and then uses the branch u ≤ w for the negative-support component. For the paper's own example (a,b,c)=(2,1,3), the true negative-component weight is w1 ≈ 0.365, but the printed rule selects w2 ≈ 0.635. A literal implementation would therefore draw from the negative branch about 63.5% of the time instead of 36.5%, producing samples inconsistent with the Lasso density. This is load-bearing because the paper's central algorithmic claim is that the inverse-CDF sampler is correct and numerically stable. The R package may or may not follow the printed rule; either way, the paper's description is wrong and must be corrected. Please also specify explicitly the four cases used when computing P^{-1}(u) in terms of the sign of b and the comparison of u with w.
  2. [§3.1] The accuracy claim for the Remez approximation of Mill's ratio is asserted without supporting evidence. The text says a degree (8,9) rational polynomial achieves 'up to 12 significant figures' on [0,600] and preserves 11 significant figures up to x ≈ 2000, but no error bound, reference to a supplementary file, or numerical verification is provided. Since this approximation is the basis for computing the normalizing constant Z and hence the inverse-CDF sampler, the claim of a 'numerically stable' algorithm depends on this accuracy. Please provide a verified error bound (e.g., a plot or table of relative error over [0,600] and beyond, with the reference interval used by the Remez optimization) or a published reference that establishes the stated accuracy.
  3. [§6, Table 1] The empirical performance comparison is based on a single run per method per dataset, with no standard errors, no number of chains, no seeds, and no statement of the hyperparameters (a1, b1, u1, v1) used in the priors. The efficiency numbers in Table 1 are therefore not statistically meaningful, and the conclusion that the modified Hans sampler is 'the most efficient' or 'second most efficient' is not supported beyond the specific runs shown. This is load-bearing for the claim of practical improvement in Bayesian lasso computation. Please report multiple independent runs (or at least a measure of variability across chains), the hyperparameter settings, and the random-number seeds, or recast the section as an illustrative comparison rather than a performance benchmark.
minor comments (5)
  1. [§3] The parameter domain is stated as a ≥ 0, but the normalizing constant, moments, and inverse CDF all use σ = 1/√a, which is undefined at a = 0. If a = 0 is intended to be included (e.g., as a Laplace limiting case), please state the limiting formulas or explicitly restrict the definition to a > 0.
  2. [Algorithm 1] Line 11 of Algorithm 1, 'RSS <- ||y||_2^2 - 2y^TXβ^(i)_j + (β^(i)_j)^T X^T X β^(i)_j', uses β^(i)_j as both a scalar (in β^(i)_j^T) and in matrix products; this is dimensionally inconsistent. The intended expression presumably involves the full vector β^(i). Please correct the notation. Similarly, line 15 writes 'y <- Xβ^(i)', but y is the response vector and should not be overwritten; the left-hand side should be the predicted vector.
  3. [§3.1] The text says 'The computation of P^{-1}(u) is carried out in four distinct cases, determined by whether u ≤ w and whether b > 0,' but the four cases are never enumerated. The subsequent sentence only describes safeguards. Please list the four cases explicitly with the corresponding formulas, especially after correcting the w1/w2 rule.
  4. [§4.3] In the caption of Figure 2, the empirical density is described as a dashed line, but the legend in the figure code and the plot description call it red; please make the caption and legend consistent.
  5. [§3] The phrase 'cumulative density function' is used in the first paragraph of Section 3; the standard term is 'cumulative distribution function'.

Circularity Check

0 steps flagged · score 0.0 of 10

No significant circularity: the distributional properties are derived from the stated density kernel, and the Gibbs application is a definitional equivalence rather than a circular prediction.

full rationale

The paper's central derivations are self-contained. The Lasso density is defined as p(X,a,b,c) = Z^-1 exp(-1/2 a X^2 + b X - c|X|), and the moments, MGF, CDF, and inverse CDF are then derived from this density in the supplement by explicit integration. These derivations do not assume the results they conclude: Z is computed by splitting the integral at zero; E(X^r) is expressed as a weighted sum of truncated-normal moments with w explicitly derived from Z; and M(t) is derived as Z(a,b+t,c)/Z(a,b,c). The Gibbs application in Section 5 is the closest thing to a definitional step: the full conditional kernel in equation (3) is exactly the Lasso density kernel, so saying the full conditional is Lasso(a,b,c) is an equivalence by construction. This is not circular reasoning; it is naming and then exploiting a kernel that arises in the model, and the sampler, exponential-family representation, and moment formulas carry independent content. Self-citations such as Ormerod et al. (2025) for the R package are code pointers rather than load-bearing evidence for the mathematical claims, and the benchmark comparisons against external packages (monomvn, bayeslm, rstan, bayesreg) provide outside validation. The possible sign-rule issue in Section 3.1 regarding the choice of w1 versus w2 when b > 0 is a correctness bug, not a circularity: it does not make any output equal to an input by construction. No fitted parameter is relabeled as a prediction, and no uniqueness theorem is imported from the authors' prior work.

Assumptions & free parameters 1 free parameters · 5 assumptions · 1 invented entities

The central properties are derived from first principles with standard calculus; no parameters are fitted to data. The Remez coefficients for Mill's ratio are numerical approximations introduced for the sampler. The model assumptions are stated in Section 2, and the identity between the full conditional and the Lasso distribution is established in the supplement.

free parameters (1)
  • Remez approximation coefficients for Mill's ratio = p0..p8 and q0..q9 as listed in Section 3.1
    Chosen by numerical optimization (Remez algorithm) to approximate Mill's ratio over [0,600]; the claimed 12-significant-figure accuracy is asserted without independent verification.
assumptions (5)
  • domain assumption Linear regression model y ~ N(Xβ, σ2 I_n) with centered responses and standardized predictors.
    Standard model setup stated in Section 2.
  • domain assumption Hierarchical prior: β_j | σ2, a_j ~ N(0, σ2/(a_j λ2)), a_j ~ IG(1, 1/2).
    Alternative parameterization in Section 2 that motivates the Lasso kernel; not the standard Park-Casella prior.
  • domain assumption The full-conditional kernel in (3) is exactly the Lasso density.
    The identity is derived in the Supplementary Material by integrating out the auxiliary variables a_j; it is the bridge between the model and the new distribution.
  • ad hoc to paper Remez rational approximation for Mill's ratio is accurate to 12 significant figures on [0,600].
    Asserted in Section 3.1 without an error bound or machine-checked proof; the sampler's tail accuracy depends on it.
  • standard math Asymptotic tail formula for Gaussian quantiles (Mächler, 2022) is valid in the extreme underflow cases.
    Used in Section 3.1 for the inverse CDF when arguments underflow; accepted from cited literature.
invented entities (1)
  • Lasso distribution
    purpose: To name and organize the full-conditional kernel for regression coefficients in Bayesian lasso and to enable efficient direct sampling.
    The density is equivalent to the kernel in Hans (2009); the novelty is the naming and derived properties, not an empirical entity with an external falsifiable handle.

how reviews work

0 comments
Cite this review

Pith. "Pith review of The Lasso Distribution: Properties, Sampling Methods, and Applications in Bayesian Lasso Regression." pith.science (2026). https://pith.science/paper/6JSSN5RE

@misc{pith2026250607394,
  author       = {Pith},
  title        = {Pith review of: The Lasso Distribution: Properties, Sampling Methods, and Applications in Bayesian Lasso Regression},
  year         = {2026},
  howpublished = {\url{https://pith.science/paper/6JSSN5RE}},
  note         = {Machine review of arXiv:2506.07394}
}
read the original abstract

In this paper, we introduce a new probability distribution, the Lasso distribution. We derive several fundamental properties of the distribution, including closed-form expressions for its moments and moment-generating function. Additionally, we present an efficient and numerically stable algorithm for generating random samples from the distribution, facilitating its use in both theoretical and applied settings. We establish that the Lasso distribution belongs to the exponential family. A direct application of the Lasso distribution arises in the context of an existing Gibbs sampler, where the full conditional distribution of each regression coefficient follows this distribution. This leads to a more computationally efficient and theoretically grounded sampling scheme. To facilitate the adoption of our methodology, we provide an R package, BayesianLasso, available on CRAN, implementing the proposed methods. Our findings offer new insights into the probabilistic structure underlying the Lasso penalty and provide practical improvements in Bayesian inference for high-dimensional regression problems.

Figures

Figures reproduced from arXiv: 2506.07394 by the authors.

Figure 1
Figure 1. The Lasso density function is depicted for different parameter values. In the top left panel, parameter a is varied while the other parameters are fixed. In the top right and bottom left panels, parameters b and c are varied, respectively, while the other parameters are fixed. as a mixture of two truncated normal distributions, the moments can be calculated as E(X r ) = (1 − w) · E(A r ) + w · E(B r ), where the mix… view at source ↗
Figure 2
Figure 2. Theoretical and empirical densities for Lasso(2, 1, 3). The theoretical density plot of Lasso(2, 1, 3) using dlasso(x, a, b, c, logarithm = FALSE) is shown as a solid black line, and the empirical density is shown as a dashed line. 4.4 Quantile function The qlasso() function computes quantiles of the Lasso distribution. Below, we compute the quantiles for probability values p = {0.1, 0.3, 0.6}: p_values <- c(0.1, 0.… view at source ↗

Discussion (0). Sign in to comment.

Reference graph

Works this paper leans on

29 extracted references · 19 canonical work pages

  1. [1]

    Bürkner, J

    P.-C. Bürkner, J. Gabry, M. Kay, and A. Vehtari. posterior: Tools for working with posterior distributions, 2023. URL https://mc-stan.org/posterior/. R package version 1.4.1

    work page 2023

  2. [2]

    Croissant and S

    Y. Croissant and S. Graves. Ecdat: Data Sets for Econometrics, 2022. URL https://CRAN.R-project.org/package=Ecdat. R package version 0.4-2

    work page 2022

  3. [3]

    Eddelbuettel

    D. Eddelbuettel. Seamless R and C++ Integration with Rcpp . Springer, New York, 2013. doi:10.1007/978-1-4614-6868-4. ISBN 978-1-4614-6867-7

    doi:10.1007/978-1-4614-6868-4 2013

  4. [4]

    Eddelbuettel and J

    D. Eddelbuettel and J. J. Balamuta. Extending R with C++ : A Brief Introduction to Rcpp . The American Statistician, 72 0 (1): 0 28--36, 2018. doi:10.1080/00031305.2017.1375990

    arXiv 2018

  5. [5]

    Eddelbuettel and R

    D. Eddelbuettel and R. Fran c ois. Rcpp : Seamless R and C++ integration. Journal of Statistical Software, 40 0 (8): 0 1--18, 2011. doi:10.18637/jss.v040.i08

    doi:10.18637/jss.v040.i08 2011

  6. [6]

    Eddelbuettel and C

    D. Eddelbuettel and C. Sanderson. Rcpparmadillo: Accelerating R with high-performance C ++ linear algebra. Computational Statistics and Data Analysis, 71: 0 1054--1063, March 2014. doi:10.1016/j.csda.2013.02.005

    doi:10.1016/j.csda.2013.02.005 2014

  7. [7]

    Eddelbuettel, R

    D. Eddelbuettel, R. Francois, J. Allaire, K. Ushey, Q. Kou, N. Russell, I. Ucar, D. Bates, and J. Chambers. Rcpp: Seamless R and C++ Integration , 2024 a . URL https://CRAN.R-project.org/package=Rcpp. R package version 1.0.13-1

    work page 2024

  8. [8]

    Eddelbuettel, R

    D. Eddelbuettel, R. Francois, D. Bates, B. Ni, and C. Sanderson. RcppArmadillo: 'Rcpp' Integration for the 'Armadillo' Templated Linear Algebra Library, 2024 b . URL https://CRAN.R-project.org/package=RcppArmadillo. R package version 14.2.2-1

    work page 2024

Show all 29 references
  1. [9]

    Efron, T

    B. Efron, T. Hastie, I. Johnstone, and R. Tibshirani. Least angle regression. The Annals of Statistics, 32 0 (2): 0 407 -- 499, 2004. doi:10.1214/009053604000000067. URL https://doi.org/10.1214/009053604000000067

    2004 doi

  2. [10]

    Gelman and D

    A. Gelman and D. B. Rubin. Inference from iterative simulation using multiple sequences. Statistical science, 7 0 (4): 0 457--472, 1992. URL https://doi.org/10.1214/ss/1177011136

    1992

  3. [11]

    G. H. Golub and C. F. van Loan. Matrix Computations. The Johns Hopkins University Press, fourth edition, 2013. ISBN 9781421407944

    2013

  4. [12]

    R. B. Gramacy. monomvn: Estimation for MVN and Student-t Data with Monotone Missingness, 2024. URL https://CRAN.R-project.org/package=monomvn. R package version 1.9-21

    2024

  5. [13]

    C. Hans. Bayesian lasso regression. Biometrika, 96 0 (4): 0 835--845, 2009. URL https://doi.org/10.1093/biomet/asp047

    2009 doi

  6. [14]

    Hastie and B

    T. Hastie and B. Efron. lars: Least Angle Regression, Lasso and Forward Stagewise, 2022. URL https://CRAN.R-project.org/package=lars. R package version 1.3

    2022

  7. [15]

    Hastie, R

    T. Hastie, R. Tibshirani, and M. Wainwright. Statistical Learning with Sparsity: The Lasso and Generalizations. Chapman & Hall/CRC, 2015. ISBN 1498712169

    2015

  8. [16]

    J. He, P. R. Hahn, H. Lopes, and A. Herren. bayeslm: Efficient Sampling for Gaussian Linear Regression with Arbitrary Priors, 2022. URL https://CRAN.R-project.org/package=bayeslm. R package version 1.0.1

    2022

  9. [17]

    Makalic and D

    E. Makalic and D. F. Schmidt. High-dimensional Bayesian regularised regression with the bayesreg package. arXiv preprint arXiv:1611.06649, 2016

    2016 arXiv

  10. [18]

    Marsaglia

    G. Marsaglia. Evaluating the normal distribution. Journal of Statistical Software, 11: 0 1--11, 2004

    2004

  11. [19]

    J. P. Mills. Table of the ratio: Area to bounding ordinate, for any protion of normal curve. Biometrika, 18 0 (3-4): 0 395--400, 1926

    1926

  12. [20]

    M. Mächler. Asymptotic Tail Formulas For Gaussian Quantiles , 2022. URL https://cran.r-project.org/package=DPQ/vignettes/qnorm-asymp.pdf

    2022

  13. [21]

    R. M. Neal. Slice sampling. The Annals of Statistics, 31 0 (3): 0 705--767, 2003

    2003

  14. [22]

    Ormerod, M

    J. Ormerod, M. J. Davoudabadi, G. Tarr, and S. Mueller. BayesianLasso : Bayesian Lasso Regression and Tools for the Lasso Distribution , 2025. URL https://github.com/garthtarr/BayesianLasso. R package on GitHub

    2025

  15. [23]

    Park and G

    T. Park and G. Casella. The B ayesian L asso. Journal of the American Statistical Association, 103 0 (482): 0 681--686, 2008

    2008

  16. [24]

    M. Redmond. Communities and Crime . UCI Machine Learning Repository, 2002. URL https://doi.org/10.24432/C53W3X

    2002 doi

  17. [25]

    Reemtsen

    R. Reemtsen. Modifications of the first R emez algorithm. SIAM Journal on Numerical Analysis, 27 0 (2): 0 507--518, 1990

    1990

  18. [26]

    RStan : the R interface to Stan , 2020

    Stan Development Team . RStan : the R interface to Stan , 2020. URL http://mc-stan.org/

    2020

  19. [27]

    LaplacesDemon: Complete Environment for Bayesian Inference , 2021

    Statisticat and LLC. LaplacesDemon: Complete Environment for Bayesian Inference , 2021. URL https://web.archive.org/web/20150206004624/http://www.bayesian-inference.com/software. R package version 16.1.6

    2021

  20. [28]

    Tibshirani

    R. Tibshirani. Regression shrinkage and selection via the lasso. Journal of the Royal Statistical Society Series B: Statistical Methodology, 58 0 (1): 0 267--288, 1996

    1996

  21. [29]

    Vehtari, A

    A. Vehtari, A. Gelman, D. Simpson, B. Carpenter, and P.-C. Bürkner. Rank-normalization, folding, and localization: An improved rhat for assessing convergence of MCMC (with discussion). Bayesian Analysis, 2021

    2021

Pith tools

Reviewed August 7, 2026 · model on record in the stance chip above.