sethaxen
commented
2 years ago I don't know if this will have any negative impact on sampling.
After further thought, yes, this is detrimental to sampling with NUTS, because it makes it very unlikely that a u-turn will be encountered, so the maximum tree depth of 10 is usually saturated:
julia> mean(==(10), chns[:tree_depth])
0.961So 96% of transitions saturated the tree depth (took 1024 leapfrog steps). Compare with the same model with Stan:
julia> using CmdStan, MCMCChains
julia> code = """
data {
int<lower=1> n;
real<lower=0> eta;
}
parameters {
corr_matrix[n] R;
}
model {
R ~ lkj_corr(eta);
}
""";
julia> stanmodel = Stanmodel(name="lkj", model=code, output_format=:mcmcchains);
...
julia> _, chns, _ = stan(stanmodel, Dict("n" => 3, "eta" => 2), "./");
...
Inference for Stan model: lkj_model
4 chains: each with iter=(1000,1000,1000,1000); warmup=(0,0,0,0); thin=(1,1,1,1); 4000 iterations saved.
Warmup took (0.013, 0.014, 0.014, 0.014) seconds, 0.055 seconds total
Sampling took (0.027, 0.026, 0.029, 0.026) seconds, 0.11 seconds total
Mean MCSE StdDev 5% 50% 95% N_Eff N_Eff/s R_hat
lp__ -1.7e+00 3.4e-02 1.4e+00 -4.5 -1.3e+00 -0.20 1.7e+03 1.6e+04 1.0e+00
accept_stat__ 0.91 1.8e-03 0.12 0.64 0.95 1.0 4.7e+03 4.4e+04 1.0e+00
stepsize__ 0.80 5.6e-02 0.080 0.67 0.85 0.88 2.0e+00 1.9e+01 2.8e+13
treedepth__ 2.2 9.1e-02 0.53 1.0 2.0 3.0 3.4e+01 3.1e+02 1.0e+00
n_leapfrog__ 4.5 3.5e-01 2.0 3.0 3.0 7.0 3.4e+01 3.1e+02 1.0e+00
divergent__ 0.00 nan 0.00 0.00 0.00 0.00 nan nan nan
energy__ 3.2 5.0e-02 1.9 0.87 2.8 6.8 1.5e+03 1.3e+04 1.0e+00
R[1,1] 1.0e+00 nan 6.7e-16 1.0 1.0e+00 1.0 nan nan nan
R[1,2] 4.0e-03 7.0e-03 4.1e-01 -0.67 -3.2e-03 0.68 3.5e+03 3.2e+04 1.0e+00
R[1,3] 2.9e-03 6.6e-03 4.1e-01 -0.68 2.2e-03 0.67 3.9e+03 3.6e+04 1.0e+00
R[2,1] 4.0e-03 7.0e-03 4.1e-01 -0.67 -3.2e-03 0.68 3.5e+03 3.2e+04 1.0e+00
R[2,2] 1.0e+00 nan 6.7e-16 1.0 1.0e+00 1.0 nan nan nan
R[2,3] -1.7e-02 7.1e-03 4.1e-01 -0.69 -1.5e-02 0.66 3.3e+03 3.1e+04 1.0e+00
R[3,1] 2.9e-03 6.6e-03 4.1e-01 -0.68 2.2e-03 0.67 3.9e+03 3.6e+04 1.0e+00
R[3,2] -1.7e-02 7.1e-03 4.1e-01 -0.69 -1.5e-02 0.66 3.3e+03 3.1e+04 1.0e+00
R[3,3] 1.0e+00 nan 6.7e-16 1.0 1.0e+00 1.0 nan nan nan
Samples were drawn using hmc with nuts.
For each parameter, N_Eff is a crude measure of effective sample size,
and R_hat is the potential scale reduction factor on split chains (at
convergence, R_hat=1).
julia> mean(==(10), chns[:treedepth__])
0.0
julia> mean(chns[:treedepth__])
2.203No transitions hit the max tree depth, and on average each transition took only 2 tree doublings (7 leapfrog steps).
torfjelde
An $n \times n$ correlation matrix has ${n \choose 2} = \frac{n (n-1)}{2}$ degrees of freedom. This is the same as the number of elements in a strict upper triangular $n \times n$ matrix. The
CorrBijectorworks by mapping from the correlation matrix first to its unique upper Cholesky factor and then to a strictly upper triangular matrix of unconstrained entries.The trouble is that in the unconstrained space, we now have $n \times n$ parameters, of which ${n+1 \choose 2} = \frac{n(n+1)}{2}$ have no impact on the log density. These extra parameters have an implicit improper uniform prior on the reals, which makes the posterior distribution in unconstrained space improper. Because these parameters have infinite variance, during adaptation, HMC will learn this, and they will explode in value. I don't know if this will have any negative impact on sampling.
In this demo, we're sampling the uniform distribution on the correlation matrices.
Note the number of parameters. We should have
3*(3-1)/2 = 3DOFs, but instead we have3*3. And note that3*(3+1)/2=6of the degrees of freedom have adapted variances of ~1e20.There are several ways to solve this, neither of which seem to be possible in Bijectors right now:
logabsdetjacto contain this prior term. However, when I tried this, it seemed to have no effect, since I guesslogabsdetjac(b::CorrBijector, X::AbstractArray{<:AbstractMatrix{<:Real}})is being called instead oflogabsdetjac(::Inverse{CorrBijector}, y::AbstractMatrix{<:Real}). When mapping fromXtoy, these extra parameters are all set to 0, so we have no way of setting this prior.