# LKJ prior throwing positivity exceptions

**URL:** <https://discourse.mc-stan.org/t/lkj-prior-throwing-positivity-exceptions/41421>\
**Category:** Modeling\
**Tags:** fitting-issues\
**Created:** [July 3, 2026, 11:44pm UTC](https://discourse.mc-stan.org/t/lkj-prior-throwing-positivity-exceptions/41421 "2026-07-03T23:44:27Z")\
**Posts on this page:** 11\
**Page:** 1

<div class="post-metadata">

**Author:** ![kyle.hofmann](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/kyle.hofmann/32/19952_2.png) [@kyle.hofmann](https://discourse.mc-stan.org/u/kyle.hofmann)\
**Post date:** [July 3, 2026, 11:44pm UTC](https://discourse.mc-stan.org/t/lkj-prior-throwing-positivity-exceptions/41421/1 "2026-07-03T23:44:27Z")

</div>

I am trying to use an LKJ prior. I have simplified my model down to the following. It’s so simple that I feel like I must be doing something wrong, yet it does not work for me:

```no-highlight
parameters {
  cholesky_factor_corr[2] L;
}
model {
  L ~ lkj_corr_cholesky(1.0);
}

```

I’m calling this `lkj_test.stan`. I run this through CmdStanPy. Again, after simplifying things as much as possible, here is exactly what I run:

```no-highlight
import cmdstanpy

model = cmdstanpy.CmdStanModel(stan_file="lkj_test.stan")
fit = model.sample({}, chains=1)
print(fit.summary())

```

The result is:

```no-highlight
16:09:46 - cmdstanpy - INFO - CmdStan start processing
chain 1: 100%|█████████████████████████████████████████████████████████████████████| 2000/2000 [00:00<00:00, 134301.53it/s, (Sampling completed)]
                                                                                                                                                 
16:09:46 - cmdstanpy - INFO - CmdStan done processing.
16:09:46 - cmdstanpy - WARNING - Non-fatal error during sampling:
Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in 'lkj_test.stan', line 5, column 2 to column 29)
        Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in 'lkj_test.stan', line 5, column 2 to column 29)
        Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in 'lkj_test.stan', line 5, column 2 to column 29)
        Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in 'lkj_test.stan', line 5, column 2 to column 29)
Consider re-running with show_console=True if the above output is unclear!
            Mean MCSE StdDev MAD 5% 50% 95% ESS_bulk ESS_tail ESS_bulk/s R_hat
lp__ -0.609607 0.047758 0.855861 0.365846 -2.327080 -0.260325 -0.001084 394.434 410.966 98608.5 1.00091
L[1,1] 1.000000 NaN 0.000000 0.000000 1.000000 1.000000 1.000000 NaN NaN NaN NaN
L[1,2] 0.000000 NaN 0.000000 0.000000 0.000000 0.000000 0.000000 NaN NaN NaN NaN
L[2,1] -0.022324 0.031549 0.572107 0.707185 -0.912745 -0.009924 0.890539 301.995 331.944 75498.8 1.00651
L[2,2] 0.788447 0.011914 0.225666 0.169100 0.312378 0.877953 0.999458 394.434 410.966 98608.5 1.00091

```

While I don’t see anything wrong with the final `L`, the exceptions are worrying. I looked at the samples, and they looked like valid Cholesky factorizations of correlation matrices; but anything invalid should have been rejected, so I don’t think that’s any evidence of anything. When I run my actual model, I get the same exception many dozens of times.

As an experiment, I also tried

```no-highlight
parameters {
  corr_matrix[2] M;
}
model {
  M ~ lkj_corr(1.0);
}

```

Calling this `lkj_test2.stan` and running it precisely as before results in:

```no-highlight
16:21:25 - cmdstanpy - INFO - CmdStan start processing
16:24:32 - cmdstanpy - INFO - CmdStan start processing
chain 1: 100%|██████████████████████████████████████████████████████████████████████| 2000/2000 [00:00<00:00, 78904.82it/s, (Sampling completed)]
                                                                                                                                                 
16:24:32 - cmdstanpy - INFO - CmdStan done processing.
16:24:32 - cmdstanpy - WARNING - Non-fatal error during sampling:
Exception: lkj_corr_lpdf: Correlation matrix is not positive definite. (in 'lkj_test2.stan', line 5, column 2 to column 20)
        Exception: lkj_corr_lpdf: Correlation matrix is not positive definite. (in 'lkj_test2.stan', line 5, column 2 to column 20)
        Exception: lkj_corr_lpdf: Correlation matrix is not positive definite. (in 'lkj_test2.stan', line 5, column 2 to column 20)
        Exception: lkj_corr_lpdf: Correlation matrix is not positive definite. (in 'lkj_test2.stan', line 5, column 2 to column 20)
Consider re-running with show_console=True if the above output is unclear!
            Mean MCSE StdDev MAD 5% 50% 95% ESS_bulk ESS_tail ESS_bulk/s R_hat
lp__ -0.663346 0.082399 1.006740 0.414771 -2.299060 -0.307206 -0.003113 425.331 212.841 53166.3 1.00199
M[1,1] 1.000000 NaN 0.000000 0.000000 1.000000 1.000000 1.000000 NaN NaN NaN NaN
M[1,2] -0.005443 0.031301 0.587459 0.766037 -0.930505 -0.008561 0.880960 265.173 212.546 33146.6 1.00977
M[2,1] -0.005443 0.031301 0.587459 0.766037 -0.930505 -0.008561 0.880960 265.173 212.546 33146.6 1.00977
M[2,2] 1.000000 NaN 0.000000 0.000000 1.000000 1.000000 1.000000 NaN NaN NaN NaN

```

Again, the final `M` and the individual samples look okay to me, but the exceptions are worrying. Also again, my actual model produces the same error, but many more times.

My question is: Am I doing something wrong? Or is this just the expected behavior when sampling from the LKJ distribution? If it is the expected behavior, is there some way I can suppress these exceptions (but not others)?

I am using Python 3.14.4 and CmdStanPy 1.3.0.

---

<div class="post-metadata">

**Author:** ![hjhornbeck](https://avatars.discourse-cdn.com/v4/letter/h/b77776/32.png) [@hjhornbeck](https://discourse.mc-stan.org/u/hjhornbeck)\
**Post date:** [July 5, 2026, 2:43am UTC](https://discourse.mc-stan.org/t/lkj-prior-throwing-positivity-exceptions/41421/2 "2026-07-05T02:43:22Z")

</div>

I don’t have a lot of experience with Cholesky matrices, but I have seen similar messages pop up when I invoke the Gamma distribution. There, Stan insists on y \> 0 for \Gamma( y | \alpha, \beta ) , so I usually define my inputs with something like `real<lower=0> y;`. Placing hard constraints on a variable can cause havok with MCMC algorithms, due to tonnes of rejected samples; Stan tries to solve this by internally sampling an unconstrained value then transforming it into a constrained one, adjusting the gradient appropriately. For positive-only values, [that transform](https://mc-stan.org/docs/reference-manual/transforms.html#lower-bound-inverse-transform) is y = \exp( \dot y ). Problem solved, right?

```stan
transformed data {
   print( "exp( -747 ) = ", exp( -747 ), "." );
}
// STDOUT: exp( -747 ) = 0.

```

Alas, floating point values only have a certain amount of precision associated with them, so in reality the hard constraint has merely shifted location. That’s still a huge improvement, most of the time. Stan’s default initialization routine picks random unconstrained values between -2 and 2, the gradient steers it away from those problematic negative values, and you never see an exception.

So if you _are_ seeing an exception, NUTS must be drifting into the Danger Zone. Unfortunately, there are quite a few potential reasons for this: you may have a broken model that favours values near zero for certain constrained variables, you may have a poorly-constrained model that doesn’t prevent the sampler from drifting, you may have not have enough data to steer the sampler away, you may be running across a corner case within Stan’s internal routines, or it may be some combination of all of the above. Even if any of the above apply, that doesn’t necessarily mean your samples are unrepresentative. Most MCMC algorithms reject far more than a dozen samples during their runtime, and yet we still consider the sample to be representative of the true distribution.

There’s really no simple solution here, you’ve got to roll up your sleeves and figure out the reason(s) yourself. One thing that may help: if I add the following to your toy example …

```stan
transformed parameters {
    if( L[1,1] == 0. )
        reject( "ERROR: invalid value in L[1,1]! Unconstrained value = ", lower_bound_unconstrain( L[1,1], 0. ), "." );
    if( L[2,2] == 0. )
        reject( "ERROR: invalid value in L[2,2]! Unconstrained value = ", lower_bound_unconstrain( L[2,2], 0. ), "." );
    if( L[2,1] == 0. )
        reject( "ERROR: invalid value in L[2,1]! Unconstrained value = ", lower_bound_unconstrain( L[2,1], 0. ), "." );
}

```

… I get different error messages.

```no-highlight
Exception: ERROR: invalid value in L[2,2]! Unconstrained value = -inf.
Exception: ERROR: invalid value in L[2,2]! Unconstrained value = -inf.
Exception: ERROR: invalid value in L[2,2]! Unconstrained value = -inf.
Exception: ERROR: invalid value in L[2,2]! Unconstrained value = -inf.

```

The ERRORs for`L[1,1]` and `L[2,1]`never pop up. I figure this could be Stan constraining the value of`L[2,2]`based on the contents of other elements in the matrix, and with no data or prior to steer it that occasionally leads to problematic values. The `-inf` is almost certainly misleading here, the underlying code is probably just \log\left( \exp\left( \dot y \right) \right) without compensation for floating point rounding. Still, it was worth showing off [`lower_bound_unconstrain()`](https://mc-stan.org/docs/functions-reference/transform_functions.html#transforms-for-scalars).

---

<div class="post-metadata">

**Author:** ![kyle.hofmann](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/kyle.hofmann/32/19952_2.png) [@kyle.hofmann](https://discourse.mc-stan.org/u/kyle.hofmann)\
**Post date:** [July 5, 2026, 4:10am UTC](https://discourse.mc-stan.org/t/lkj-prior-throwing-positivity-exceptions/41421/3 "2026-07-05T04:10:17Z")

</div>

I’m not sure that explains the problem; but I tried digging through the source, and the behavior I’m seeing looks impossible, so I must be misunderstanding something.

I found `‎lkj_corr_cholesky_rng` in [math/stan/math/prim/prob/lkj\_corr\_cholesky\_rng.hpp](https://github.com/stan-dev/math/blob/master/stan/math/prim/prob/lkj_corr_cholesky_rng.hpp#L13). When `K = 2`, this creates an array `CPCs` with only a single entry. It sets `alpha = eta + 0.5` and then runs a pair of nested loops. Both loops execute only once. In the only loop iteration, `alpha` gets updated to `eta` and `CPCs[0]` is set to `2.0 * beta_rng(alpha, alpha, rng) - 1.0`. For my choice of \eta = 1, this is a \operatorname{Beta}(1, 1)-distributed random variable shifted and scaled to [-1, 1], i.e., it’s uniform on [-1, 1].

As the last step of `‎lkj_corr_cholesky_rng`, `CPCs` is sent to `read_corr_L` in [stan/math/prim/fun/read\_corr\_L.hpp](https://github.com/stan-dev/math/blob/master/stan/math/prim/fun/read_corr_L.hpp#L38). This creates an array `acc` with a single entry and a 2\times 2 matrix `L`. It first sets `L` to zero. Then, if I’m reading this right, it sets (using C++ indexing) `L[0,0] = 1.0`, `L[1,0] = CPCs[0]`, `acc[0] = 1.0 - CPCs[0].square()`, and `L[1,1] = sqrt(acc[0])`. That is, if `CPCs[0]` is c, it produces the matrix L = \begin{pmatrix} 1 & 0 \\ c & \sqrt{1 - c^2} \end{pmatrix}. This is a valid Cholesky factorization of a correlation matrix; indeed, it’s easy to calculate LL^T = \begin{pmatrix} 1 & c \\ c & 1 \end{pmatrix}.

The error messages your suggestion yields seem informative with regards to what variable is causing the exceptions, but not with regards to why. The `-inf`s appearing in your error messages look erroneous to me. Ultimately, `lower_bound_unconstrain` calls `lb_free` in [stan/math/prim/constraint/lb\_free.hpp](https://github.com/stan-dev/math/blob/master/stan/math/prim/constraint/lb_free.hpp#L30). This function assumes that it’s been passed exponentiated values, so it takes a log. The fact that it’s reporting `-inf` means that it’s actually complaining about `L[2,2]` being zero.

The exception I’m seeing is coming from `lkj_corr_cholesky_lpdf` in [stan/math/prim/prob/lkj\_corr\_cholesky\_lpdf.hpp](https://github.com/stan-dev/math/blob/master/stan/math/prim/prob/lkj_corr_cholesky_lpdf.hpp). Early on, this calls `check_cholesky_factor` in [stan/math/prim/err/check\_cholesky\_factor.hpp](https://github.com/stan-dev/math/blob/master/stan/math/prim/err/check_cholesky_factor.hpp). The check that fails is the last one in the function, `check_positive` applied to the diagonal. And that function is in [stan/math/prim/err/check\_positive.hpp](https://github.com/stan-dev/math/blob/master/stan/math/prim/err/check_positive.hpp), and all it does is check that what it’s passed is (strictly) greater than zero.

`L[2,2]` is a square root of a non-negative number. Mathematically, it’s almost surely positive. Numerically, there’s some tiny but positive probability that `c` is `1.0` and that we get zero exactly, but that’s extraordinarily unlikely unless \eta \ll 1. Which it’s not. That leaves me still confused: I can’t see how this is happening, and yet it is.

---

<div class="post-metadata">

**Author:** ![spinkney](https://avatars.discourse-cdn.com/v4/letter/s/dec6dc/32.png) [@spinkney](https://discourse.mc-stan.org/u/spinkney)\
**Post date:** [July 6, 2026, 1:26pm UTC](https://discourse.mc-stan.org/t/lkj-prior-throwing-positivity-exceptions/41421/4 "2026-07-06T13:26:54Z")

</div>

This is almost surely an overflow issue where the value is overflowing to 1 in which causes the diagonal to be 0 and triggers the warning. As @kyle.hofmann noted the mathematics say this is impossible but with finite-precision computing this happens.

In the original model, the dimension of the correlation matrix is only 2. I almost never use the built-in Cholesky factor of correlation matrices for this. Since you can easily parameterize it with 1 parameter yourself as

```stan
real<lower=-1, upper=1> phi;

```

We, the Stan devs, need to update the parameterization to something like [Updated cholesky corr parameterization testing](https://discourse.mc-stan.org/t/updated-cholesky-corr-parameterization-testing/38827) which shouldn’t have this issue.

---

<div class="post-metadata">

**Author:** ![hjhornbeck](https://avatars.discourse-cdn.com/v4/letter/h/b77776/32.png) [@hjhornbeck](https://discourse.mc-stan.org/u/hjhornbeck)\
**Post date:** [July 6, 2026, 10:01pm UTC](https://discourse.mc-stan.org/t/lkj-prior-throwing-positivity-exceptions/41421/5 "2026-07-06T22:01:33Z")

</div>

> I found `‎lkj_corr_cholesky_rng` in [math/stan/math/prim/prob/lkj\_corr\_cholesky\_rng.hpp](https://github.com/stan-dev/math/blob/master/stan/math/prim/prob/lkj_corr_cholesky_rng.hpp#L13).

I’m not sure why you’re looking there. The PRNG sampling routines might be relevant during initialization, but if `check_positive()` on line 15 was the trigger then adding a `transformed parameters` section wouldn’t have changed the exception message. Line 22 in [lkj\_corr\_cholesky\_lpdf.hpp](https://github.com/stan-dev/math/blob/master/stan/math/prim/prob/lkj_corr_cholesky_lpdf.hpp) is more likely.

> The `-inf`s appearing in your error messages look erroneous to me. Ultimately, `lower_bound_unconstrain` calls `lb_free` in [stan/math/prim/constraint/lb\_free.hpp](https://github.com/stan-dev/math/blob/master/stan/math/prim/constraint/lb_free.hpp#L30). This function assumes that it’s been passed exponentiated values, so it takes a log. The fact that it’s reporting `-inf` means that it’s actually complaining about `L[2,2]` being zero.

Yep, I said as much here:

> The `-inf` is almost certainly misleading here, the underlying code is probably just log⁡(exp⁡(˙𝑦)) without compensation for floating point rounding.

While that might have been a bit too subtle, the conditional was `( L[2,2] == 0. )` .

> That leaves me still confused: I can’t see how this is happening, and yet it is.

Ah, my inexperience with Cholesky correlation matrices is showing. As you and @spinkney pointed out, `L[2,1]` must be between -1 and 1. Thus Stan uses [this transformation](https://mc-stan.org/docs/reference-manual/transforms.html#log-odds-and-the-logistic-sigmoid) to go from unconstrained to constrained values …

y = \left( 1 + \exp(-\dot y) \right)^{-1}

… which has precision issues on _two_ sides. To make what @spinkney said a bit more concrete:

```stan
    print( "inv_logit( -747 ) = ", inv_logit( -747 ), "." ); // STDOUT: inv_logit( -747 ) = 0.
    print( "inv_logit( 747 ) = ", inv_logit( 747 ), "." ); // STDOUT: inv_logit( 747 ) = 1.

```

Again, whether that’s an issue or not depends on your model. Does it require more precision than an IEEE double-precision float gives you? If so, then you’ve either got to get creative or ditch Stan.

---

<div class="post-metadata">

**Author:** ![kyle.hofmann](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/kyle.hofmann/32/19952_2.png) [@kyle.hofmann](https://discourse.mc-stan.org/u/kyle.hofmann)\
**Post date:** [July 7, 2026, 2:04am UTC](https://discourse.mc-stan.org/t/lkj-prior-throwing-positivity-exceptions/41421/6 "2026-07-07T02:04:49Z")

</div>

> [@hjhornbeck](#):
>
> I’m not sure why you’re looking there.

I assumed that it would be relevant, but I didn’t put the pieces together to see that it wasn’t related to where the exception is actually being thrown. Sorry about that; this is my first time using Stan.

The transformation between unconstrained and constrained seems to happen at [stan/math/prim/constraint/cholesky\_corr\_constrain.hpp#L17](https://github.com/stan-dev/math/blob/master/stan/math/prim/constraint/cholesky_corr_constrain.hpp#L17). This calls `corr_constrain` from [stan/math/prim/constraint/corr\_constrain.hpp#L27](https://github.com/stan-dev/math/blob/master/stan/math/prim/constraint/corr_constrain.hpp#L27), and `corr_constrain` uses \tanh as its transformation. So not quite the logit and its inverse (it can’t be those because they map between [0, 1] and \mathbf{R}, while this function needs to map between [-1, 1] and \mathbf{R}), but different only by some shifting and scaling. And in standard double precision arithmetic, \tanh(19) = 1.0, so we don’t have to wander very far before we’re in trouble.

I’m still rather concerned that I get these exceptions every time even in the trivial case of my original post. My actual data has non-trivial correlation (which is _a priori_ expected; it’s one of those “when people like something they do it more” phenomena), so the sampler is going to gravitate towards larger values of the correlation coefficient, and I guess that makes numerical issues more likely. My plan, after I figured this out, was to ask about translating the LKJ distribution so that I could put its peak somewhere other than the identity. But I’m starting to think it might be safer to switch to a Wishart or maybe inverse Wishart distribution. I don’t have a strong preference among them, and for this application, I suspect they’ll all work fine in the end.

> [@spinkney](#):
>
> I almost never use the built-in Cholesky factor of correlation matrices for this.

Okay, that’s good to know. I didn’t know whether parameterizing by hand like that would be liable to cause issues elsewhere. Thanks!

---

<div class="post-metadata">

**Author:** ![Bob\_Carpenter](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bob_carpenter/32/9230_2.png) [@Bob\_Carpenter](https://discourse.mc-stan.org/u/Bob_Carpenter)\
**Post date:** [July 8, 2026, 7:03pm UTC](https://discourse.mc-stan.org/t/lkj-prior-throwing-positivity-exceptions/41421/7 "2026-07-08T19:03:24Z")

</div>

> [@kyle.hofmann](#):
>
> this is my first time using Stan.

Welcome! This is a pretty ambitious place to start.

> [@kyle.hofmann](#):
>
> so the sampler is going to gravitate towards larger values of the correlation coefficient, and I guess that makes numerical issues more likely.

This shouldn’t be a problem unless the correlations get close to 1 or the variances get close to 0. More information will tend to keep the results away from the boundaries as the posteriors concentrate away from uniform.

> [@kyle.hofmann](#):
>
> Am I doing something wrong? Or is this just the expected behavior when sampling from the LKJ distribution?

No. Yes.

> [@kyle.hofmann](#):
>
> anything invalid should have been rejected, so I don’t think that’s any evidence of anything

The rejections will cause duplicated states in the sampler, but that’s not the only cause.

We intentionally label these as “Non-fatal error” because it’s recoverable. What’s usually going on is that early in sampling you will get these exceptiosn or you will get them occassionally during sampling. They can indicate there might be bias if they’re all due to the same cause and show up as divergences in the same region of the posterior.

> [@hjhornbeck](#):
>
> the gradient steers it away from those problematic negative values,

Negative values themselves are necessary in order to get values less than 1. The problem starts when the absolute value of unconstrained variable exceeds 350 or so and starts underflowing 0 (negative) or overflowing to +infinity (positive).

> [@hjhornbeck](#):
>
> you may have a broken model that favours values near zero for certain constrained variables

This isn’t necessarily a broken model so much as our sampler not being up to the job. @andrewgelman likes to recommend things like half normal priors on scale parameters in hierarchical models. If the data is consistent with complete pooling, this can run into numerical difficulties compared to a “zero avoiding” prior like lognormal. So while Bayes wants a really small value in these cases, the `exp()` transform can be a problem.

In another topic here, @aseyboldt posted a [more stable altnerative he calls “bisoftplus”](https://discourse.mc-stan.org/t/models-where-stan-outperforms-nutpie-walnuts/41095/33). Like we did for sum-to-zero and simplex, we may wind up changing our default built-in transform as a result.

> [@hjhornbeck](#):
>
> The ERRORs for`L[1,1]` and `L[2,1]` never pop up

If `L[2,2]` is bad, you wont see messages for `L[2,1]` because rejections stop where they’re at.

> [@spinkney](#):
>
> I almost never use the built-in Cholesky factor of correlation matrices for this.

Right, you take the correlation `phi`, and the scales `sigma[1]` and `sigma[2]` and create

```stan
parameters {
  vector<lower=0>[2] sigma;
  real<lower=-1, upper=1> phi;
}
transformed parameters {
  matrix[2, 2] Sigma;
  {
    vector[2] var = sigma^2;
    real cov = sigma[1] * sigma[2] * phi;
    Sigma = [[var[1], cov], [cov, var[2]];
  }
}

```

---

<div class="post-metadata">

**Author:** ![kyle.hofmann](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/kyle.hofmann/32/19952_2.png) [@kyle.hofmann](https://discourse.mc-stan.org/u/kyle.hofmann)\
**Post date:** [July 10, 2026, 6:16am UTC](https://discourse.mc-stan.org/t/lkj-prior-throwing-positivity-exceptions/41421/8 "2026-07-10T06:16:24Z")

</div>

> [@Bob\_Carpenter](#):
>
> This is a pretty ambitious place to start.

It’s the ambitious stuff that gives me enough of a kick to learn something new. But the learning curve has been more gentle than I expected!

> [@Bob\_Carpenter](#):
>
> No. Yes.

Okay, great. Thanks, that helps.

> [@Bob\_Carpenter](#):
>
> a [more stable alternative he calls “bisoftplus”](https://discourse.mc-stan.org/t/models-where-stan-outperforms-nutpie-walnuts/41095/33)

Oh, that’s quite interesting. I like the fact that it’s basically linear for positive x and doesn’t go to zero too quickly for negative x. I looked for a better parameterization, but I didn’t find one. (My attempt was: The minimal polynomial of y = x + \sqrt{x^2 + 1} is y^2 - 2xy - 1. Over \mathbf{C} it defines a Riemann surface which double covers the sphere and is branched over \pm i. It’s a rational algebraic curve, so it has a rational parameterization. For example, fix a point on the curve like (0, 1) and consider the line of slope t through (0, 1), i.e., y = tx + 1; we can solve this pair of equations to express any one of x, y, or t in terms of the other two. But when we express y in terms of t, we find that y(t) is a parabola, and that’s not a helpful way to parameterize the positive reals for this application. And when we express x in terms of t, half the slopes correspond to the wrong branch, i.e,. to y = x - \sqrt{x^2 + 1}, so we end up with half the slopes being disallowed. That’s kinda like where we started with half the real numbers being disallowed, so it isn’t an improvement either. Consequently this idea seems to make things no better.)

In any case, I am busily sampling as we speak. A test run on a subset of my data produced good-looking results, so I’m feeling pretty good right now. Thanks again for your help!

---

<div class="post-metadata">

**Author:** ![Bob\_Carpenter](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bob_carpenter/32/9230_2.png) [@Bob\_Carpenter](https://discourse.mc-stan.org/u/Bob_Carpenter)\
**Post date:** [July 10, 2026, 3:04pm UTC](https://discourse.mc-stan.org/t/lkj-prior-throwing-positivity-exceptions/41421/9 "2026-07-10T15:04:50Z")

</div>

I’m afraid I don’t know what minimal polynomials are or where t came from.

From a change-of-variables perspective, this should be fine. Let me be super pedantic about what we’re proposing.

In case anyone wants to implement, here are the details in the style of the _Reference Manual_ chapter on transforms. Suppose we have a variable Y \in (0, \infty) and define

\qquad X = f(Y) = (Y^2 - 1) / (2 \cdot Y).

The inverse is

\qquad Y = f^{-1}(X) + \sqrt{X^2 + 1}.

Given the density of the constrained Y, p\_Y(y), the unconstrained density of X is

\qquad p\_X(x) = p\_Y(f^{-1}(x)) \cdot \left| \nabla f^{-1}(x) \right|,

which reduces to

\qquad p\_X(x) = p\_Y\!\left(x + \sqrt{x^2 + 1}\right) \cdot \left(1 + \dfrac{x}{\sqrt{x^2 + 1}}\right).

---

<div class="post-metadata">

**Author:** ![kyle.hofmann](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/kyle.hofmann/32/19952_2.png) [@kyle.hofmann](https://discourse.mc-stan.org/u/kyle.hofmann)\
**Post date:** [July 11, 2026, 2:16am UTC](https://discourse.mc-stan.org/t/lkj-prior-throwing-positivity-exceptions/41421/10 "2026-07-11T02:16:02Z")

</div>

Ah, well, you can safely ignore what I said! I was hoping to find something that behaved similarly but had a different formula. There is nothing wrong with the given one, but I thought I had an idea, so I worked it out. The idea failed—the result didn’t solve the original problem of mapping between (0, \infty) and \mathbf{R}. So while I found it to be a fun exercise, it’s not something anyone else ought to spend serious time with!

(But in case you’re curious: The minimal polynomial of x + \sqrt{x^2 + 1} is the smallest degree polynomial with leading coefficient 1 such that plugging in x + \sqrt{x^2 + 1} gives zero. To write down that polynomial, we introduce a new variable y. The polynomial y^2 - 2xy - 1 is quadratic in y, has leading coefficient 1, and if we substitute y = x + \sqrt{x^2 + 1}, we get zero. There’s no degree 1 polynomial with this property because of the square root, so y^2 - 2xy - 1 is the minimal polynomial of x + \sqrt{x^2 + 1}. If you’re familiar with field theory and Galois theory, this is the minimal polynomial of x + \sqrt{x^2 + 1} over the field \mathbf{R}(x) of rational functions with real coefficients. This polynomial turns up in what you wrote, because solving y^2 - 2xy - 1 = 0 for x gives x = (y^2 - 1)/(2y) just like you had; minimality says that this expression for x in terms of y is, in a certain limited sense, as simple as possible. A degree 2 polynomial is small enough that this is quite special, and the theory of Riemann surfaces made me hope that I could do something fun and clever. That was where t came in; but alas, it was not meant to be!)

Oh, and a last follow-up: So far, my results look great! Thanks for developing Stan; without it, what I wanted to do would have been insurmountable, and instead it has been easy.

---

<div class="post-metadata">

**Author:** ![Bob\_Carpenter](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bob_carpenter/32/9230_2.png) [@Bob\_Carpenter](https://discourse.mc-stan.org/u/Bob_Carpenter)\
**Post date:** [July 14, 2026, 7:00pm UTC](https://discourse.mc-stan.org/t/lkj-prior-throwing-positivity-exceptions/41421/11 "2026-07-14T19:00:14Z")

</div>

> [@kyle.hofmann](#):
>
> If you’re familiar with field theory and Galois theory

Funny that you should mention that. Way back in 1982 as a pure math undergrad, I was going to sign up for the matrices class, but my advisor (a combinatorialist) told me that matrices were for engineers and I could skip it. So I learned Galois theory in my abstract algebra class. Why not learn about morphisms of rings and general linear operators rather than concrete ones? It was a bit of an uphill when I came to stats without knowing what a determinant was.

> [@kyle.hofmann](#):
>
> Thanks for developing Stan; without it, what I wanted to do would have been insurmountable, and instead it has been easy.

You’re welcome, and thanks for the kind words. It’s what keeps us going. That and we use it ourselves.