# Compare two Turing models with LOO: SE is NaN

**URL:** <https://discourse.julialang.org/t/compare-two-turing-models-with-loo-se-is-nan/102269>\
**Category:** Probabilistic Programming\
**Tags:** turing, bayesian-inference\
**Created:** [July 30, 2023, 2:30pm UTC](https://discourse.julialang.org/t/compare-two-turing-models-with-loo-se-is-nan/102269 "2023-07-30T14:30:57Z")\
**Posts on this page:** 8\
**Page:** 1

<div class="post-metadata">

**Author:** ![DominiqueMakowski](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dominiquemakowski/32/51410_2.png) [@DominiqueMakowski](https://discourse.julialang.org/u/DominiqueMakowski)\
**Post date:** [July 30, 2023, 2:30pm UTC](https://discourse.julialang.org/t/compare-two-turing-models-with-loo-se-is-nan/102269/1 "2023-07-30T14:30:57Z")

</div>

I am trying to compare two [SequentialSamplingModels](https://github.com/itsdfish/SequentialSamplingModels.jl) with loo using ParetoSmooth.jl, but I think I am doing something wrong:

```julia
using Turing
using SequentialSamplingModels
using Random
using Distributions
using DataFrames
using StatsPlots
using StatsModels
using StatsBase
using ParetoSmooth

# Generate data (1000 obs)
Random.seed!(6)

dist = LBA(ν=[3.0, 2.0], A=0.8, k=0.2, τ=0.3)
data = rand(dist, 1000)

# ---------------------
# Models
@model function model_lba(data; min_rt=minimum(data.rt))
    ν ~ filldist(Normal(0, 1), 2)
    A ~ truncated(Normal(0.8, 0.4), 0.0, Inf)
    k ~ truncated(Normal(0.2, 0.2), 0.0, Inf)
    τ ~ Uniform(0.0, min_rt)

    data ~ LBA(; ν, A, k, τ)
end
chain_lba = sample(model_lba(data), NUTS(), 1000)

@model function model_lnr(data; min_rt=minimum(data.rt))
    ν ~ filldist(Normal(0, 1), 2)
    σ ~ truncated(Normal(0, 1), 0.0, Inf)
    τ ~ Uniform(0.0, min_rt)

    data ~ LNR(; ν, σ, τ)
end
chain_lnr = sample(model_lnr(data), NUTS(), 1000)

```

When I run the following:

```julia
rez1 = psis_loo(model_lba(data), chain_lba)
rez2 = psis_loo(model_lnr(data), chain_lnr)
loo_compare((lba=rez1, lnr=rez2))

```

While the comparison works, the `psis_loo()` function returns

```julia
┌ Warning: Some Pareto k values are extremely high (>1). PSIS will not produce consistent estimates.
└ @ ParetoSmooth C:\Users\domma\.julia\packages\ParetoSmooth\Ml7Gb\src\InternalHelpers.jl:47
Results of PSIS-LOO-CV with 1000 Monte Carlo samples and 1 data points. Total Monte Carlo SE of NaN.
┌───────────┬────────┬──────────┬────────┬─────────┐
│ │ total │ se_total │ mean │ se_mean │
├───────────┼────────┼──────────┼────────┼─────────┤
│ cv_elpd │ 452.85 │ NaN │ 452.85 │ NaN │
│ naive_lpd │ 456.28 │ NaN │ 456.28 │ NaN │
│ p_eff │ 3.43 │ NaN │ 3.43 │ NaN │

```

I am not sure where it gets the “1 data points” from, as there are 1000 observations 🤔

Moreover, in [this post](https://discourse.mc-stan.org/t/interpreting-elpd-diff-loo-package/1628/2), Aki mentions that it is useful to compare the ELPD relative to their SEs to have an idea of the magnitude of the difference. Hence I am wondering if this info (or some standardized difference) is available or can be computed? Thanks for any tips for model comparison!

---

<div class="post-metadata">

**Author:** ![Christopher\_Fisher](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/christopher_fisher/32/26132_2.png) [@Christopher\_Fisher](https://discourse.julialang.org/u/Christopher_Fisher)\
**Post date:** [July 30, 2023, 2:59pm UTC](https://discourse.julialang.org/t/compare-two-turing-models-with-loo-se-is-nan/102269/2 "2023-07-30T14:59:16Z")

</div>

> [@DominiqueMakowski](#):
>
> `(model_lba(data), chain_lba)`

I think the problem is that `pointwise_log_likelihoods` does not compute pointwise correctly for this type of model. The following should be 1000X1000X1 I believe

`pointwise_log_likelihoods(model_lba(data), chain_lba)`

I’ll continue digging to see where the error occurs.

---

<div class="post-metadata">

**Author:** ![Christopher\_Fisher](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/christopher_fisher/32/26132_2.png) [@Christopher\_Fisher](https://discourse.julialang.org/u/Christopher_Fisher)\
**Post date:** [July 30, 2023, 4:20pm UTC](https://discourse.julialang.org/t/compare-two-turing-models-with-loo-se-is-nan/102269/3 "2023-07-30T16:20:37Z")

</div>

I thought the issue was that [nsamples](https://github.com/JuliaStats/Distributions.jl/blob/c8d3e4b52ea2c04f655510af85f93fa3876f258f/src/common.jl#L112) was not defined. Adding that definition didn’t fix the problem. I’m not sure where in Turing the data size is computed.

---

<div class="post-metadata">

**Author:** ![DominiqueMakowski](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dominiquemakowski/32/51410_2.png) [@DominiqueMakowski](https://discourse.julialang.org/u/DominiqueMakowski)\
**Post date:** [July 30, 2023, 5:01pm UTC](https://discourse.julialang.org/t/compare-two-turing-models-with-loo-se-is-nan/102269/4 "2023-07-30T17:01:49Z")

</div>

It works if we use the list of tuples specification:

```julia
@model function model_lba(data; min_rt=0.2)
    # Priors
    ν ~ filldist(Normal(0, 1), 2)
    A ~ truncated(Normal(0.8, 0.4), 0.0, Inf)
    k ~ truncated(Normal(0.2, 0.2), 0.0, Inf)
    τ ~ Uniform(0.0, min_rt)

    # Likelihood
    for i in 1:length(data)
        data[i] ~ LBA(; ν, A, k, τ)
    end
end

dat = [(choice=data.choice[i], rt=data.rt[i]) for i in 1:length(data.rt)]
chain_lba = sample(model_lba(dat, min_rt=minimum(data.rt)), NUTS(), 1000)

rez1 = psis_loo(model_lba(dat, min_rt=minimum(data.rt)), chain_lba)

```

```julia
[ Info: No source provided for samples; variables are assumed to be from a Markov Chain. If the samples are independent, specify this with keyword argument `source=:other`.
Results of PSIS-LOO-CV with 1000 Monte Carlo samples and 1000 data points. Total Monte Carlo SE of 0.084.
┌───────────┬────────┬──────────┬───────┬─────────┐
│ │ total │ se_total │ mean │ se_mean │
├───────────┼────────┼──────────┼───────┼─────────┤
│ cv_elpd │ 453.32 │ 27.43 │ 0.45 │ 0.03 │
│ naive_lpd │ 457.57 │ 27.26 │ 0.46 │ 0.03 │
│ p_eff │ 4.25 │ 0.30 │ 0.00 │ 0.00 │

```

---

<div class="post-metadata">

**Author:** ![Christopher\_Fisher](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/christopher_fisher/32/26132_2.png) [@Christopher\_Fisher](https://discourse.julialang.org/u/Christopher_Fisher)\
**Post date:** [July 30, 2023, 5:42pm UTC](https://discourse.julialang.org/t/compare-two-turing-models-with-loo-se-is-nan/102269/5 "2023-07-30T17:42:52Z")

</div>

That makes sense. Thanks for reporting. My guess is that `length` or something similar is called to extract the length of the vector.

---

<div class="post-metadata">

**Author:** ![Christopher\_Fisher](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/christopher_fisher/32/26132_2.png) [@Christopher\_Fisher](https://discourse.julialang.org/u/Christopher_Fisher)\
**Post date:** [July 30, 2023, 7:57pm UTC](https://discourse.julialang.org/t/compare-two-turing-models-with-loo-se-is-nan/102269/6 "2023-07-30T19:57:53Z")

</div>

I thought maybe the problem is that SequentialSamplingModels was not complying with the interface, but the problem occurs with other models. Consider the following:

```julia
using Turing
using ParetoSmooth
using Distributions

@model function model(data)
    μ ~ Normal()
    data ~ Normal(μ, 1)
end

data = rand(Normal(0, 1), 100)

chain = sample(model(data), NUTS(), 1000)
rez1 = psis_loo(model(data), chain)

```

Replacing the vectorized form with a for loop fixes the problem. I wonder whether the Turing should destructure the arrays to compute pointwise correctly. Of course, that may lead to other problems.

---

<div class="post-metadata">

**Author:** ![Christopher\_Fisher](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/christopher_fisher/32/26132_2.png) [@Christopher\_Fisher](https://discourse.julialang.org/u/Christopher_Fisher)\
**Post date:** [August 2, 2023, 10:02am UTC](https://discourse.julialang.org/t/compare-two-turing-models-with-loo-se-is-nan/102269/7 "2023-08-02T10:02:25Z")

</div>

Using a for loop is currently the proper approach for using LOO with ParetoSmooth.jl. I will make a PR to explain that in the documentation. There might be a plan to improve the interface at some point in the future. Another alternative is Arviz.jl.

---

<div class="post-metadata">

**Author:** ![penelopeysm](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/penelopeysm/32/213172_2.png) [@penelopeysm](https://discourse.julialang.org/u/penelopeysm)\
**Post date:** [April 19, 2026, 1:17am UTC](https://discourse.julialang.org/t/compare-two-turing-models-with-loo-se-is-nan/102269/8 "2026-04-19T01:17:20Z")

</div>

> [@Christopher\_Fisher](#):
>
> Replacing the vectorized form with a for loop fixes the problem. I wonder whether the Turing should destructure the arrays to compute pointwise correctly.

I’m aware I’m necroing an old post, but since DynamicPPL 0.40 (specifically [this PR](https://github.com/TuringLang/DynamicPPL.jl/pull/1279)) this model will return a vector of log-likelihoods for `data` so will work correctly with LOO.

```julia
using Turing
@model function f(data)
    μ ~ Normal()
    data ~ Normal(μ, 1)
end

```

It is likely that even more forms of models can be supported in the future (e.g. `data ~ MvNormal(...)`, see [Implement factorised pointwise probabilities by penelopeysm · Pull Request #1331 · TuringLang/DynamicPPL.jl · GitHub](https://github.com/TuringLang/DynamicPPL.jl/pull/1331)).

Now the catch is that ParetoSmooth is not really maintained and so won’t work with DynamicPPL 0.40 as of the time of writing. (Technically, ParetoSmooth doesn’t have to talk directly to DynamicPPL: they can communicate ‘via’ the chain of log-likelihoods, but you would have to keep them in separate environments and serialise the chain, which is a bit silly.)

However, you can use FlexiChains.jl + PosteriorStats.jl, which will Just Work:

```julia
using FlexiChains, PosteriorStats
m = f(randn(10))
chn = sample(m, NUTS(), 100; chain_type=VNChain)
result = PosteriorStats.loo(m, chn)

```