{"id":"47191979-4c89-4742-baee-7fbe0027fb80","arxiv_id":"2506.07394","paper_version":3,"verdict":"CONDITIONAL","confidence":"MODERATE","novelty_score":4.0,"correctness_risk":"medium","formal_verification":"none","parameter_count":1,"one_line_summary":"The authors define the Lasso distribution, derive its properties including an inverse-CDF sampler, and use it to build a faster Bayesian lasso Gibbs sampler in the BayesianLasso R package.","lead":"The paper introduces a named 'Lasso distribution' for the regression coefficient full conditional in Bayesian lasso, with derived moments, a moment generating function, and a numerically stable sampler. It packages this as an R package and claims faster Bayesian lasso Gibbs sampling on benchmark datasets.","discovery_kind":"extension","skeptic_critique":{"model":"deepseek-v4-flash","headline":"The Section 3.1 sign rule for w1/w2 is backwards: for b>0 the negative-component weight is w1, not w2, so the inverse-CDF sampler as printed draws the wrong mixture branch.","rationale":"The reader's weakest assumption identifies the same load-bearing flaw: the sign-dependent rule for w1 and w2 in Section 3.1 is reversed. My independent calculation confirms the reversal: for b > 0 the negative-support mixture weight equals w1, not w2, and the paper's own example parameters make the error quantitatively large (0.635 vs 0.365). This is a genuine internal inconsistency, not a matter of outside consensus. Because the inverse-CDF sampler is a central claimed contribution and is used in the Gibbs application, the manuscript cannot be accepted as written. The fix is small—either swap the labels or swap the branch condition—so the correct disposition is conditional acceptance pending correction and verification against the package. I do not see a second concern of comparable weight: the density, normalizing constant, MGF, and moment derivations in the supplementary material are essentially standard calculations and appear internally consistent; the Remez/Mill's-ratio approximation is asserted without a published error bound, but that is a reproducibility detail rather than a structural flaw. The performance comparison is a single run without error bars, but that weakens the empirical claims without undermining the central algorithmic claim. Therefore the reader's CONDITIONAL verdict is appropriate and my assessment does not move it.","tokens_in":13952,"tokens_out":5884,"duration_ms":66961,"concrete_test":"Download the BayesianLasso source from CRAN/GitHub and inspect the implementation of rlasso()/qlasso(): locate where w is assigned from w1 and w2 and whether the branch condition uses u <= w. Independently compute, for (a,b,c) = (2,1,3), w1 = 1/(1 + m(v1)/m(v2)) and w2 = 1/(1 + m(v2)/m(v1)) with m(x) = pnorm(-x)/dnorm(x), and compare with the true negative-component weight. Then generate 10^6 draws with rlasso(1e6, 2, 1, 3) and compare the empirical CDF against plasso() on a fine grid; if the printed sign rule is implemented, a KS or chi-square test will reject dramatically, whereas the corrected rule should give agreement within Monte Carlo error.","verdict_should_be":"UNCHANGED","load_bearing_attack":"Section 3.1 of the manuscript states w = w1 if b <= 0 and w = w2 if b > 0, then uses the branch 'u <= w' for the negative-support component. This is reversed. From Section 3, w is defined as the coefficient of E(B^r) in the moment formula, i.e. the weight of the negative-support truncated-normal component B ~ TN-(mu2, sigma^2). In terms of v1 = (c-|b|)/sqrt(a) and v2 = (c+|b|)/sqrt(a), that component's weight is m(v2)/(m(v1)+m(v2)) when b > 0, which is exactly w1 as defined in Section 3.1, not w2. For the paper's own example (a,b,c) = (2,1,3), v1 = sqrt(2), v2 = 2*sqrt(2), m(v1) ~ 0.545, m(v2) ~ 0.314, so the true negative-component weight is w1 ~ 0.365. The printed rule sets w = w2 ~ 0.635, so a sampler implemented literally would draw from the negative branch about 63.5% of the time instead of 36.5%, producing samples inconsistent with the Lasso density. The closed-form density, moments, and MGF derivations do not depend on this sampling rule and appear consistent, but the central algorithmic claim—an efficient numerically stable inverse-CDF sampler—is not correct as written. The R package may or may not follow the text; either way the paper's description is wrong.","agreement_with_reader":"agree"},"referee_report":{"model":"deepseek-v4-flash","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.","tokens_in":14314,"tokens_out":3526,"duration_ms":44755,"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":[{"comment":"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.","section":"§3.1"},{"comment":"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.","section":"§3.1"},{"comment":"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.","section":"§6, Table 1"}],"minor_comments":[{"comment":"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.","section":"§3"},{"comment":"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.","section":"Algorithm 1"},{"comment":"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.","section":"§3.1"},{"comment":"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.","section":"§4.3"},{"comment":"The phrase 'cumulative density function' is used in the first paragraph of Section 3; the standard term is 'cumulative distribution function'.","section":"§3"}],"recommendation":"major_revision","confidential_remarks":"The sign error in Section 3.1 is easily verified numerically, and the authors should also check the R package implementation against the corrected rule, ideally adding a stochastic test (e.g., compare empirical quantiles from rlasso() against the theoretical quantile function) to the package tests. The performance comparison in Section 6 is too thin for the strength of the claims made; I would suggest the authors either add a proper multi-run comparison or soften the conclusions accordingly."},"author_rebuttal":null,"desk_editor":{"model":"deepseek-v4-flash","letter":"Here's my read. The paper does real work: it takes the full-conditional kernel that Hans (2009) used for Bayesian lasso and treats it as a named distribution, deriving the normalizing constant, CDF, inverse CDF, moments, MGF, and exponential-family representation. I checked the derivations in the supplementary material; they're consistent. The Mill's ratio approximation for the normalizing constant is a practical touch, and the R package is a genuine deliverable. The modified Hans Gibbs sampler is a sensible way to use the distribution, and the complexity analysis is honest.\n\nThe soft spot is the sampling rule in Section 3.1. The paper defines w = w1 for b ≤ 0 and w2 for b > 0, and sets 1-w accordingly. That's reversed. From the moment formula, w is the weight on the negative-support truncated normal. For b > 0, the negative component's weight is w1 = m(v2)/(m(v1)+m(v2)), not w2. The paper's own example (a,b,c)=(2,1,3) gives w1 ≈ 0.365, but the printed rule sends you to w2 ≈ 0.635. An inverse-CDF sampler implemented literally from the text would therefore draw from the wrong mixture branch. This may be a typo—the package code might follow the correct rule—but as printed, the central algorithmic claim is wrong. The derivation of w itself in Section 3 is fine; only the sign-based renaming in 3.1 is flipped.\n\nAlso, the benchmark table is a single run per dataset, no error bars, and no hyperparameter details for the rival packages. That's enough for a pilot comparison but not for strong efficiency claims. The Remez approximation's 12-significant-figure accuracy is asserted without an error bound; minor.\n\nOn novelty: calling the Lasso distribution 'new' is generous—it's the Hans kernel renamed—but the closed-form properties and the numerically stable sampler are not in Hans, and they're useful.\n\nBottom line: this is a serious paper that deserves refereeing. It needs a major revision to fix the sign rule (and confirm the package matches), and the benchmark reproducibility should be tightened. If the package is correct, the fix is small; if the package follows the text, the sampler is biased and the paper is not usable as written. Either way, the referee will need to see the corrected algorithm and a reproducibility statement. I'd send it out.","headline":"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.","tokens_in":14774,"tokens_out":4566,"would_cite":false,"duration_ms":46754,"reading_group":"maybe","serious_thinker":"yes","would_accept_peer_review":true},"rs_alignment":null,"lean_confirmation":null,"pith_extraction":{"msc":["62E15","62F15","62J07","65C10"],"pacs":[],"model":"deepseek-v4-flash","headline":"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.","keywords":["Lasso distribution","Bayesian lasso","Gibbs sampler","truncated normal mixture","inverse-CDF sampling","Mills ratio","exponential family","high-dimensional regression"],"falsifier":"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.","tokens_in":13789,"feed_emoji":"📊","tokens_out":5930,"duration_ms":56796,"temperature":0.7,"pith_summary":"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.","feed_headline":"Bayesian lasso gains a named distribution and a faster sampler","feed_subtitle":"Closed-form moments and an inverse-CDF sampler turn each coefficient update into one direct draw.","key_machinery":"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.","core_discovery":"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.","pith_inferences":["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."],"forward_implications":["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."],"supporting_citations":[{"why":"Supplies the Bayesian lasso Gibbs sampler whose coefficient full conditional is the kernel the Lasso distribution normalizes; the paper modifies this sampler.","marker":"Hans (2009)"},{"why":"Provides the scale-mixture Bayesian lasso formulation whose alternative parameterization motivates the Lasso distribution.","marker":"Park and Casella (2008)"},{"why":"Defines lasso regression, the target model this work builds on.","marker":"Tibshirani (1996)"},{"why":"Defines the Mills ratio in which the normalizing constant of the Lasso distribution is expressed.","marker":"Mills (1926)"},{"why":"Gives the Remez algorithm used to construct the rational approximation of the Mills ratio that underpins the sampling scheme.","marker":"Reemtsen (1990)"},{"why":"Provides the asymptotic tail formula used to evaluate the inverse normal CDF in extreme tails of the Lasso distribution.","marker":"Mächler (2022)"},{"why":"Supplies the slice sampler used for the variance and penalty parameters in the proposed Gibbs algorithms.","marker":"Neal (2003)"},{"why":"Documents the precision limits of evaluating the normal density in double precision, motivating the numerical safeguards in the implementation.","marker":"Marsaglia (2004)"}],"fun_headline_variants":["Lasso distribution: closed-form moments and a direct sampler","Bayesian lasso coefficients drawn directly from new Lasso distribution","New Lasso distribution speeds up Bayesian lasso Gibbs sampler","One direct draw per coefficient: the Lasso distribution in action","Lasso distribution: exponential family, closed forms, and R package"],"cache_read_input_tokens":3200,"weakest_assumption_plain":"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.","fun_headline_variants_meta":{"raw":{"variants":["Lasso distribution: closed-form moments and a direct sampler","Bayesian lasso coefficients drawn directly from new Lasso distribution","New Lasso distribution speeds up Bayesian lasso Gibbs sampler","One direct draw per coefficient: the Lasso distribution in action","Lasso distribution: exponential family, closed forms, and R package"]},"model":"deepseek-v4-flash","effort":"low","cost_usd":0.00068,"raw_usage":{"total_tokens":3067,"prompt_tokens":903,"completion_tokens":2164,"prompt_tokens_details":{"cached_tokens":384},"prompt_cache_hit_tokens":384,"prompt_cache_miss_tokens":519,"completion_tokens_details":{"reasoning_tokens":2079}},"tokens_in":519,"tokens_out":2164,"duration_ms":17638,"temperature":1.0,"reasoning_tokens":2079,"cache_read_input_tokens":384,"cache_creation_input_tokens":0},"cache_creation_input_tokens":0},"created_at":"2026-08-07T05:36:09.196825+00:00","model_set":{"reader":"deepseek-v4-flash"},"falsifier":"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.","supporting_citations":[],"review_version":1}