From 9333a541fb3e8dcd987b206f7e82a820f60ae65c Mon Sep 17 00:00:00 2001 From: Seth Axen Date: Tue, 29 Sep 2026 00:18:21 +0200 Subject: [PATCH 01/34] refactor: Call StatsFuns directly for distribution functions --- src/discretediag.jl | 8 ++++---- src/gelmandiag.jl | 2 +- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/src/discretediag.jl b/src/discretediag.jl index 0c0071dc..97d39d2d 100644 --- a/src/discretediag.jl +++ b/src/discretediag.jl @@ -71,7 +71,7 @@ function weiss(X::AbstractMatrix) stat = (n / ca) * sum(chi_stat) pval = NaN if ((m_tot - 1) * (d - 1)) >= 1 - pval = Distributions.ccdf(Distributions.Chisq((m_tot - 1) * (d - 1)), stat) + pval = StatsFuns.chisqccdf((m_tot - 1) * (d - 1), stat) end return (stat, m_tot, pval, ca) @@ -303,14 +303,14 @@ function diag_all( stat = t * sum(chi_stat) df0 = (m - 1) * (d - 1) if m > 1 && !isnan(stat) - pval = Distributions.ccdf(Distributions.Chisq(df0), stat) + pval = StatsFuns.chisqccdf(df0, stat) end elseif method == :weiss stat = (t / ca) * sum(chi_stat) df0 = (m - 1) * (d - 1) pval = NaN if m > 1 && !isnan(stat) - pval = Distributions.ccdf(Distributions.Chisq(df0), stat) + pval = StatsFuns.chisqccdf(df0, stat) end elseif method == :DARBOOT stat = t * sum(chi_stat) @@ -339,7 +339,7 @@ function diag_all( stat = hot_stat df0 = df if df > 0 && !isnan(hot_stat) - pval = Distributions.ccdf(Distributions.Chisq(df), hot_stat) + pval = StatsFuns.chisqccdf(df, hot_stat) end elseif method == :billingsleyBOOT stat = hot_stat diff --git a/src/gelmandiag.jl b/src/gelmandiag.jl index dd7b3ded..73db1906 100644 --- a/src/gelmandiag.jl +++ b/src/gelmandiag.jl @@ -44,7 +44,7 @@ function _gelmandiag(psi::AbstractArray{<:Real,3}; alpha::Real=0.05) estimates[i] = sqrt(correction * (rfixed + rrandom)) if !isnan(rrandom) - rrandom *= Distributions.quantile(Distributions.FDist(B_df, W_df[i]), q) + rrandom *= StatsFuns.fdistinvcdf(B_df, W_df[i], q) end upperlimits[i] = sqrt(correction * (rfixed + rrandom)) end From 4c6e068db48c84e5f17ad324f55f36cc9ea50dbe Mon Sep 17 00:00:00 2001 From: Seth Axen Date: Tue, 29 Sep 2026 00:20:27 +0200 Subject: [PATCH 02/34] feat: Use our own categorical sampler instead of Distributions --- src/discretediag.jl | 23 ++++++++++++----------- src/utils.jl | 17 +++++++++++++++++ 2 files changed, 29 insertions(+), 11 deletions(-) diff --git a/src/discretediag.jl b/src/discretediag.jl index 97d39d2d..01feccda 100644 --- a/src/discretediag.jl +++ b/src/discretediag.jl @@ -185,7 +185,7 @@ function bd_inner(Y::AbstractMatrix, m::Int) end @doc raw""" - simulate_DAR1!(X::Matrix{Int}, phi::Float64, sampler) + simulate_DAR1!(X::Matrix{Int}, phi::Float64, prob::AbstractVector{<:Real}) Simulate a DAR(1) model independently in each column of `X`. @@ -196,30 +196,31 @@ X_t = \alpha_t X_{t-1} + (1 - \alpha_t) \epsilon_{t-1}, where ```math \begin{aligned} -X_1 \sim \text{sampler}, \\ +X_1 \sim \operatorname{Categorical}(\mathrm{prob}), \\ \alpha_t \sim \operatorname{Bernoulli}(\phi), \\ -\epsilon_{t-1} \sim \text{sampler}, +\epsilon_{t-1} \sim \operatorname{Categorical}(\mathrm{prob}), \end{aligned} ``` are independent random variables. """ -function simulate_DAR1!(X::Matrix{Int}, phi::Float64, sampler) +function simulate_DAR1!(X::Matrix{Int}, phi::Float64, prob::AbstractVector{<:Real}) n = size(X, 1) n > 0 || error("output matrix must be non-empty") + rng = Random.default_rng() # for all simulations @inbounds for j in axes(X, 2) # sample first value from categorical distribution with probabilities `prob` - X[1, j] = rand(sampler) + X[1, j] = _wsample(rng, prob) for t in 2:n # compute next value - X[t, j] = if rand() <= phi + X[t, j] = if rand(rng) <= phi # copy previous value with probability `phi` X[t - 1, j] else # sample value with probability `1-phi` - rand(sampler) + _wsample(rng, prob) end end end @@ -230,9 +231,10 @@ end function simulate_MC(N::Int, P::Matrix{Float64}) X = zeros(Int, N) n, m = size(P) - X[1] = StatsBase.sample(1:n) + rng = Random.default_rng() + X[1] = rand(rng, 1:n) for i in 2:N - X[i] = StatsBase.wsample(1:n, vec(P[X[i - 1], :])) + X[i] = _wsample(rng, vec(P[X[i - 1], :])) end return X end @@ -314,11 +316,10 @@ function diag_all( end elseif method == :DARBOOT stat = t * sum(chi_stat) - sampler_phat = Distributions.sampler(Distributions.Categorical(phat)) bstats = zeros(nsim) Y = Matrix{Int}(undef, t, d) for b in 1:nsim - simulate_DAR1!(Y, phia, sampler_phat) + simulate_DAR1!(Y, phia, phat) s = hangartner_inner(Y, m)[1] @inbounds bstats[b] = s end diff --git a/src/utils.jl b/src/utils.jl index 0fb3f721..481a545f 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -183,6 +183,23 @@ function _rank_normalize!(values, x) return values end +# sample an index `i` of `weights` with probability proportional to `weights[i]`. +# for `Float64` weights, equivalent to `StatsBase.wsample(rng, eachindex(weights), weights)` +function _wsample(rng::Random.AbstractRNG, weights::AbstractVector{T}) where {T<:Real} + t = rand(rng, float(T)) * sum(weights) + i = firstindex(weights) + ilast = lastindex(weights) + cw = weights[i] + while true + cw ≥ t && return i + # may happen with floating point weights due to numerical inaccuracies. + # never return an index with zero weight. + i == ilast && return findlast(!iszero, weights, i) + i += 1 + cw += weights[i] + end +end + # transform the ranks to quantiles of a standard normal distribution applying the # "α-β correction" recommended in Eq 6.10.3 of # Blom. Statistical Estimates and Transformed Beta-Variables. Wiley; New York, 1958 From c5320b2dc093a4272a7063323f8138c53064edcd Mon Sep 17 00:00:00 2001 From: Seth Axen Date: Tue, 29 Sep 2026 00:37:52 +0200 Subject: [PATCH 03/34] feat!: Move rstar Distributions integration to extension --- Project.toml | 6 ++++++ ext/MCMCDiagnosticToolsDistributionsExt.jl | 16 +++++++++++++++ src/MCMCDiagnosticTools.jl | 7 +++++++ src/rstar.jl | 24 +++++++++++++++------- 4 files changed, 46 insertions(+), 7 deletions(-) create mode 100644 ext/MCMCDiagnosticToolsDistributionsExt.jl diff --git a/Project.toml b/Project.toml index 98a82a85..7ee7f6cc 100644 --- a/Project.toml +++ b/Project.toml @@ -20,6 +20,12 @@ StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91" StatsFuns = "4c63d2b9-4356-54db-8cca-17b64c39e42c" Tables = "bd369af6-aec1-5ad0-b16a-f7cc5008161c" +[weakdeps] +Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" + +[extensions] +MCMCDiagnosticToolsDistributionsExt = "Distributions" + [compat] AbstractFFTs = "0.5, 1" DataAPI = "1.12" diff --git a/ext/MCMCDiagnosticToolsDistributionsExt.jl b/ext/MCMCDiagnosticToolsDistributionsExt.jl new file mode 100644 index 00000000..5e38cf20 --- /dev/null +++ b/ext/MCMCDiagnosticToolsDistributionsExt.jl @@ -0,0 +1,16 @@ +module MCMCDiagnosticToolsDistributionsExt + +using MCMCDiagnosticTools: MCMCDiagnosticTools +using Distributions: Distributions + +function MCMCDiagnosticTools._rstar_distribution(predictions, ytest, nclasses::Int) + # create Poisson binomial distribution with support `0:length(predictions)` + distribution = Distributions.PoissonBinomial(map(Distributions.pdf, predictions, ytest)) + + # scale distribution to support in `[0, nclasses]` + scaled_distribution = (nclasses//length(predictions)) * distribution + + return scaled_distribution +end + +end diff --git a/src/MCMCDiagnosticTools.jl b/src/MCMCDiagnosticTools.jl index ecd049b7..beb754d9 100644 --- a/src/MCMCDiagnosticTools.jl +++ b/src/MCMCDiagnosticTools.jl @@ -35,4 +35,11 @@ include("heideldiag.jl") include("mcse.jl") include("rafterydiag.jl") include("rstar.jl") + +function __init__() + if isdefined(Base.Experimental, :register_error_hint) + Base.Experimental.register_error_hint(_rstar_error_hint, MethodError) + end +end + end diff --git a/src/rstar.jl b/src/rstar.jl index f4f2e6ff..2a1f46f1 100644 --- a/src/rstar.jl +++ b/src/rstar.jl @@ -147,7 +147,8 @@ within-chain convergence. The training of the classifier can be inspected by adj If the classifier is deterministic, i.e., if it predicts a class, the value of the ``R^*`` statistic is returned (algorithm 1). If the classifier is probabilistic, i.e., if it outputs probabilities of classes, the scaled Poisson-binomial distribution of the ``R^*`` statistic -is returned (algorithm 2). +is returned (algorithm 2). Probabilistic classifiers require +[Distributions.jl](https://github.com/JuliaStats/Distributions.jl) to be loaded. !!! note The correctness of the statistic depends on the convergence of the `classifier` used @@ -164,6 +165,7 @@ julia> samples = fill(4.0, 100, 3, 2); One can compute the distribution of the ``R^*`` statistic (algorithm 2) with a probabilistic classifier. For instance, we can use a gradient-boosted trees model with `nrounds = 100` sequentially stacked trees and learning rate `eta = 0.05`: +Because MLJBase already loads Distributions, we don't need to do so ourselves. ```jldoctest rstar julia> model = EvoTreeClassifier(; nrounds=100, eta=0.05); @@ -254,14 +256,22 @@ function _rstar( length(predictions) == length(ytest) || error("numbers of predictions and targets must be equal") - # create Poisson binomial distribution with support `0:length(predictions)` - distribution = Distributions.PoissonBinomial(map(Distributions.pdf, predictions, ytest)) - - # scale distribution to support in `[0, nclasses]` nclasses = length(MMI.classes(ytest)) - scaled_distribution = (nclasses//length(predictions)) * distribution + return _rstar_distribution(predictions, ytest, nclasses) +end + +# Implemented in MCMCDiagnosticToolsDistributionsExt +function _rstar_distribution end - return scaled_distribution +function _rstar_error_hint(io, exc, argtypes, kwargs) + if exc.f === _rstar_distribution + print( + io, + "\n`rstar` with a probabilistic classifier requires Distributions to be loaded, " * + "e.g. with `using Distributions`.", + ) + end + return nothing end # unsupported types of predictions and targets From 0a2092ffc35b906454ada46070721cda7f47fd6c Mon Sep 17 00:00:00 2001 From: Seth Axen Date: Tue, 29 Sep 2026 00:38:46 +0200 Subject: [PATCH 04/34] chore: Remove Distributions as dependency --- Project.toml | 1 - src/MCMCDiagnosticTools.jl | 1 - 2 files changed, 2 deletions(-) diff --git a/Project.toml b/Project.toml index 7ee7f6cc..782eccbe 100644 --- a/Project.toml +++ b/Project.toml @@ -10,7 +10,6 @@ projects = ["docs", "test"] AbstractFFTs = "621f4979-c628-5d54-868e-fcf4e3e8185c" DataAPI = "9a962f9c-6df0-11e9-0e5d-c546b8b5ee8a" DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8" -Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" MLJModelInterface = "e80e1ace-859a-464e-9ed9-23947d8ae3ea" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" diff --git a/src/MCMCDiagnosticTools.jl b/src/MCMCDiagnosticTools.jl index beb754d9..a4f49021 100644 --- a/src/MCMCDiagnosticTools.jl +++ b/src/MCMCDiagnosticTools.jl @@ -3,7 +3,6 @@ module MCMCDiagnosticTools using AbstractFFTs: AbstractFFTs using DataAPI: DataAPI using DataStructures: DataStructures -using Distributions: Distributions using MLJModelInterface: MLJModelInterface as MMI using SpecialFunctions: SpecialFunctions using StatsBase: StatsBase From 1601eae865456cc6ceac42b93538cc6295a3ddb3 Mon Sep 17 00:00:00 2001 From: Seth Axen Date: Tue, 29 Sep 2026 00:46:40 +0200 Subject: [PATCH 05/34] feat: Remove DataStructures as dependency --- Project.toml | 2 -- src/MCMCDiagnosticTools.jl | 1 - src/utils.jl | 6 +++--- 3 files changed, 3 insertions(+), 6 deletions(-) diff --git a/Project.toml b/Project.toml index 782eccbe..d6a567db 100644 --- a/Project.toml +++ b/Project.toml @@ -9,7 +9,6 @@ projects = ["docs", "test"] [deps] AbstractFFTs = "621f4979-c628-5d54-868e-fcf4e3e8185c" DataAPI = "9a962f9c-6df0-11e9-0e5d-c546b8b5ee8a" -DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" MLJModelInterface = "e80e1ace-859a-464e-9ed9-23947d8ae3ea" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" @@ -28,7 +27,6 @@ MCMCDiagnosticToolsDistributionsExt = "Distributions" [compat] AbstractFFTs = "0.5, 1" DataAPI = "1.12" -DataStructures = "0.18.13, 0.19" Distributions = "0.25.57" LinearAlgebra = "1.6" MLJModelInterface = "1.11" diff --git a/src/MCMCDiagnosticTools.jl b/src/MCMCDiagnosticTools.jl index a4f49021..45653a97 100644 --- a/src/MCMCDiagnosticTools.jl +++ b/src/MCMCDiagnosticTools.jl @@ -2,7 +2,6 @@ module MCMCDiagnosticTools using AbstractFFTs: AbstractFFTs using DataAPI: DataAPI -using DataStructures: DataStructures using MLJModelInterface: MLJModelInterface as MMI using SpecialFunctions: SpecialFunctions using StatsBase: StatsBase diff --git a/src/utils.jl b/src/utils.jl index 481a545f..459b3f10 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -50,7 +50,7 @@ found. function unique_indices(x) inds = eachindex(x) T = eltype(inds) - ind_map = DataStructures.SortedDict{eltype(x),Vector{T}}() + ind_map = Dict{eltype(x),Vector{T}}() for i in inds xi = x[i] inds_xi = get!(ind_map, xi) do @@ -58,8 +58,8 @@ function unique_indices(x) end push!(inds_xi, i) end - unique = collect(keys(ind_map)) - indices = collect(values(ind_map)) + unique = sort!(collect(keys(ind_map))) + indices = map(Base.Fix1(getindex, ind_map), unique) return unique, indices end From 0d3df583d435945526d49c0b2bceabfea3365725 Mon Sep 17 00:00:00 2001 From: Seth Axen Date: Tue, 29 Sep 2026 00:53:01 +0200 Subject: [PATCH 06/34] fix: fix bug in _wsample --- src/utils.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/utils.jl b/src/utils.jl index 459b3f10..0e6e4bbc 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -194,7 +194,7 @@ function _wsample(rng::Random.AbstractRNG, weights::AbstractVector{T}) where {T< cw ≥ t && return i # may happen with floating point weights due to numerical inaccuracies. # never return an index with zero weight. - i == ilast && return findlast(!iszero, weights, i) + i == ilast && return something(findlast(!iszero, weights), i) i += 1 cw += weights[i] end From a3d6ee963a30645474d21f4f2b4144eafccae870 Mon Sep 17 00:00:00 2001 From: Seth Axen Date: Tue, 29 Sep 2026 00:55:05 +0200 Subject: [PATCH 07/34] test: Add _wsample tests --- test/utils.jl | 74 +++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 74 insertions(+) diff --git a/test/utils.jl b/test/utils.jl index db853627..96828c7f 100644 --- a/test/utils.jl +++ b/test/utils.jl @@ -3,6 +3,7 @@ using Test using OffsetArrays using Random using Statistics +using StatsBase: StatsBase @testset "unique_indices" begin @testset "indices=$(eachindex(inds))" for inds in [ @@ -106,6 +107,79 @@ end end end +# RNG that always generates the same `Float64`, used to control `_wsample` +struct FixedRNG <: Random.AbstractRNG + x::Float64 +end +function Random.rand(rng::FixedRNG, ::Random.SamplerTrivial{Random.CloseOpen01{Float64}}) + return rng.x +end + +@testset "_wsample" begin + @testset "basic" begin + w = [0.2, 0.0, 0.5, 0.3] + # the first index whose cumulative weight is at least `rand() * sum(w)` is sampled + @test MCMCDiagnosticTools._wsample(FixedRNG(0.0), w) == 1 + @test MCMCDiagnosticTools._wsample(FixedRNG(0.1), w) == 1 + @test MCMCDiagnosticTools._wsample(FixedRNG(0.3), w) == 3 + @test MCMCDiagnosticTools._wsample(FixedRNG(0.9), w) == 4 + @test @inferred(MCMCDiagnosticTools._wsample(Random.default_rng(), w)) isa Int + end + + @testset "all-zero weights return the first index" begin + @test MCMCDiagnosticTools._wsample(FixedRNG(0.5), zeros(3)) == 1 + end + + @testset "never returns an index with zero weight" begin + # emulate floating-point error causing `rand() * sum(w)` to exceed the cumulative + # sum of the weights + w = [0.2, 0.5, 0.3, 0.0, 0.0] + rng = FixedRNG(nextfloat(1.0)) + @test MCMCDiagnosticTools._wsample(rng, w) == 3 + @test MCMCDiagnosticTools._wsample(rng, w) == + StatsBase.wsample(rng, eachindex(w), w) + end + + @testset "offset indices" begin + w = OffsetArray([0.2, 0.0, 0.5, 0.3], -2) + @test MCMCDiagnosticTools._wsample(FixedRNG(0.0), w) == -1 + @test MCMCDiagnosticTools._wsample(FixedRNG(0.3), w) == 1 + end + + @testset "consistent with StatsBase.wsample" begin + p = rand(10) + p ./= sum(p) + weights = ( + [1.0], + rand(5), + [rand(5); zeros(3)], + [0.0; rand(3); 0.0; rand(2)], + p, + rand(0:3, 8), + ) + for w in weights + seed = rand(UInt) + rng1, rng2 = Xoshiro(seed), Xoshiro(seed) + @test all(1:100) do _ + return MCMCDiagnosticTools._wsample(rng1, w) == + StatsBase.wsample(rng2, eachindex(w), w) + end + end + end + + @testset "sampling frequencies" begin + w = [1.0, 0.0, 3.0, 6.0] + rng = Xoshiro(42) + ndraws = 100_000 + counts = zeros(Int, length(w)) + for _ in 1:ndraws + counts[MCMCDiagnosticTools._wsample(rng, w)] += 1 + end + @test counts[2] == 0 + @test counts ./ ndraws ≈ w ./ sum(w) atol = 0.01 + end +end + @testset "_fold_around_median" begin @testset for sz in ((1000,), (1000, 4), (1000, 4, 8), (1000, 4, 8, 2)) x = rand(sz...) From f0a1550ff834ad22e56f3e16ec45832c90d1defb Mon Sep 17 00:00:00 2001 From: Seth Axen Date: Tue, 29 Sep 2026 00:55:47 +0200 Subject: [PATCH 08/34] feat: Roll our own _tiedrank! --- src/utils.jl | 29 +++++++++++++++++++++++++++-- test/utils.jl | 46 ++++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 73 insertions(+), 2 deletions(-) diff --git a/src/utils.jl b/src/utils.jl index 0e6e4bbc..4a8e6997 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -177,12 +177,37 @@ function _rank_normalize!(values, x) fill!(values, missing) return values end - rank = StatsBase.tiedrank(x) - _normal_quantiles_from_ranks!(values, rank) + _tiedrank!(values, x) + _normal_quantiles_from_ranks!(values, values) map!(StatsFuns.norminvcdf, values, values) return values end +# rank `x` in ascending order, assigning tied values the mean of their ranks. +# equivalent to `StatsBase.tiedrank(x)` +function _tiedrank!(ranks::AbstractArray, x::AbstractArray) + _tiedrank!(vec(ranks), vec(x)) + return ranks +end +function _tiedrank!(ranks::AbstractVector, x::AbstractVector) + perm = sortperm(x) + n = lastindex(perm) + start = firstindex(perm) + offset = 1 - start + while start ≤ n + xstart = x[perm[start]] + # find the last index `stop` of the run of values tied with xstart + stop = start + while stop < n && x[perm[stop + 1]] == xstart + stop += 1 + end + rank = (start + stop) / 2 + offset + @views ranks[perm[start:stop]] .= rank + start = stop + 1 + end + return ranks +end + # sample an index `i` of `weights` with probability proportional to `weights[i]`. # for `Float64` weights, equivalent to `StatsBase.wsample(rng, eachindex(weights), weights)` function _wsample(rng::Random.AbstractRNG, weights::AbstractVector{T}) where {T<:Real} diff --git a/test/utils.jl b/test/utils.jl index 96828c7f..025dc01d 100644 --- a/test/utils.jl +++ b/test/utils.jl @@ -107,6 +107,52 @@ end end end +@testset "_tiedrank!" begin + @testset "basic" begin + x = [3.0, 1.0, 2.0, 1.0, 5.0, 1.0] + ranks = similar(x) + @test @inferred(MCMCDiagnosticTools._tiedrank!(ranks, x)) === ranks + @test ranks == [5, 2, 4, 2, 6, 2] + @test isempty(MCMCDiagnosticTools._tiedrank!(Float64[], Float64[])) + end + + @testset "special values" begin + # as in StatsBase, -0.0 and 0.0 are tied, while NaNs are ranked last and not tied + x = [NaN, 0.0, Inf, -0.0, NaN, -Inf, 1.0] + ranks = MCMCDiagnosticTools._tiedrank!(similar(x), x) + @test ranks == [6, 2.5, 5, 2.5, 7, 1, 4] + @test ranks == StatsBase.tiedrank(x) + end + + @testset "consistent with StatsBase.tiedrank" begin + @testset for T in (Float64, Float32, Int), sz in ((20,), (20, 3)) + x = T <: Integer ? rand(T(-3):T(3), sz) : round.(randn(T, sz); digits=1) + expected = StatsBase.tiedrank(x) + @test MCMCDiagnosticTools._tiedrank!(similar(x, Float64), x) == expected + # ranks with the same float type as `x`, as in `_rank_normalize!` + @test MCMCDiagnosticTools._tiedrank!(similar(x, float(T)), x) == expected + end + end + + @testset "array types" begin + x = round.(randn(20, 3, 2); digits=1) + # parameter slices, as passed by `_rank_normalize` + xslice = view(x, :, :, 2) + @test MCMCDiagnosticTools._tiedrank!(similar(xslice), xslice) == + StatsBase.tiedrank(xslice) + + xoff = OffsetArray(x[:, :, 1], -3, 5) + ranks = MCMCDiagnosticTools._tiedrank!(similar(xoff), xoff) + @test axes(ranks) == axes(xoff) + @test parent(ranks) == StatsBase.tiedrank(parent(xoff)) + + xvoff = OffsetArray(x[:, 1, 1], -5) + ranks = MCMCDiagnosticTools._tiedrank!(similar(xvoff), xvoff) + @test axes(ranks) == axes(xvoff) + @test parent(ranks) == StatsBase.tiedrank(parent(xvoff)) + end +end + # RNG that always generates the same `Float64`, used to control `_wsample` struct FixedRNG <: Random.AbstractRNG x::Float64 From b54e546ad08f35af6cc3d93bd3a660b6eaa8d7bd Mon Sep 17 00:00:00 2001 From: Seth Axen Date: Tue, 29 Sep 2026 00:56:18 +0200 Subject: [PATCH 09/34] feat: Roll our own _counts function --- src/rafterydiag.jl | 4 ++-- src/utils.jl | 13 +++++++++++++ test/utils.jl | 15 +++++++++++++++ 3 files changed, 30 insertions(+), 2 deletions(-) diff --git a/src/rafterydiag.jl b/src/rafterydiag.jl index a6960838..e1a651a7 100644 --- a/src/rafterydiag.jl +++ b/src/rafterydiag.jl @@ -44,7 +44,7 @@ function rafterydiag( test = dichot[1:kthin:nx] ntest = length(test) temp = test[1:(ntest - 2)] + 2 * test[2:(ntest - 1)] + 4 * test[3:ntest] - trantest = reshape(StatsBase.counts(temp, 0:7), 2, 2, 2) + trantest = reshape(_counts(temp, 0:7), 2, 2, 2) g2 = 0.0 for i1 in 1:2, i2 in 1:2, i3 in 1:2 tt = trantest[i1, i2, i3] @@ -58,7 +58,7 @@ function rafterydiag( bic = g2 - 2.0 * log(ntest - 2.0) end - tranfinal = StatsBase.counts(test[1:(ntest - 1)] + 2 * test[2:ntest], 0:3) + tranfinal = _counts(test[1:(ntest - 1)] + 2 * test[2:ntest], 0:3) alpha = tranfinal[3] / (tranfinal[1] + tranfinal[3]) beta = tranfinal[2] / (tranfinal[2] + tranfinal[4]) kthin *= step(range) diff --git a/src/utils.jl b/src/utils.jl index 4a8e6997..2c6ac9a8 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -225,6 +225,19 @@ function _wsample(rng::Random.AbstractRNG, weights::AbstractVector{T}) where {T< end end +# count the number of occurrences of each value of `levels` in `x`. +# equivalent to `StatsBase.counts(x, levels)` +function _counts(x::AbstractArray{<:Integer}, levels::UnitRange{<:Integer}) + counts = zeros(Int, length(levels)) + offset = 1 - first(levels) + for xi in x + if xi in levels + counts[xi + offset] += 1 + end + end + return counts +end + # transform the ranks to quantiles of a standard normal distribution applying the # "α-β correction" recommended in Eq 6.10.3 of # Blom. Statistical Estimates and Transformed Beta-Variables. Wiley; New York, 1958 diff --git a/test/utils.jl b/test/utils.jl index 025dc01d..186a12f1 100644 --- a/test/utils.jl +++ b/test/utils.jl @@ -226,6 +226,21 @@ end end end +@testset "_counts" begin + x = [0, 3, 3, 7, 8, -1, 3, 0] + @test @inferred(MCMCDiagnosticTools._counts(x, 0:7)) == [2, 0, 0, 3, 0, 0, 0, 1] + @test MCMCDiagnosticTools._counts(x, 3:3) == [3] + @test MCMCDiagnosticTools._counts(reshape(x, 2, 4), 0:7) == [2, 0, 0, 3, 0, 0, 0, 1] + @test MCMCDiagnosticTools._counts(Int[], 0:3) == zeros(Int, 4) + + @testset "consistent with StatsBase.counts" begin + @testset for levels in (0:3, 0:7, -3:5) + x = rand(-5:9, 50) + @test MCMCDiagnosticTools._counts(x, levels) == StatsBase.counts(x, levels) + end + end +end + @testset "_fold_around_median" begin @testset for sz in ((1000,), (1000, 4), (1000, 4, 8), (1000, 4, 8, 2)) x = rand(sz...) From eeb5d456e6fca97075b162d7784ed22029e0e06e Mon Sep 17 00:00:00 2001 From: Seth Axen Date: Tue, 29 Sep 2026 00:56:37 +0200 Subject: [PATCH 10/34] fix: Call Statistics.quantile directly --- src/rafterydiag.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/rafterydiag.jl b/src/rafterydiag.jl index e1a651a7..639a3000 100644 --- a/src/rafterydiag.jl +++ b/src/rafterydiag.jl @@ -35,7 +35,7 @@ function rafterydiag( kthin = -1 burnin = total = NaN else - dichot = Int[(x .<= StatsBase.quantile(x, q))...] + dichot = Int[(x .<= Statistics.quantile(x, q))...] kthin = 0 bic = 1.0 local test, ntest From f11c8727e640cbbfa508b6099f984acf2cc1e15e Mon Sep 17 00:00:00 2001 From: Seth Axen Date: Tue, 29 Sep 2026 00:58:13 +0200 Subject: [PATCH 11/34] refactor: Move expectand proxy for mad to extension --- Project.toml | 2 ++ ext/MCMCDiagnosticToolsStatsBaseExt.jl | 12 ++++++++++++ src/ess_rhat.jl | 4 ---- 3 files changed, 14 insertions(+), 4 deletions(-) create mode 100644 ext/MCMCDiagnosticToolsStatsBaseExt.jl diff --git a/Project.toml b/Project.toml index d6a567db..73c9a714 100644 --- a/Project.toml +++ b/Project.toml @@ -20,9 +20,11 @@ Tables = "bd369af6-aec1-5ad0-b16a-f7cc5008161c" [weakdeps] Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" +StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91" [extensions] MCMCDiagnosticToolsDistributionsExt = "Distributions" +MCMCDiagnosticToolsStatsBaseExt = "StatsBase" [compat] AbstractFFTs = "0.5, 1" diff --git a/ext/MCMCDiagnosticToolsStatsBaseExt.jl b/ext/MCMCDiagnosticToolsStatsBaseExt.jl new file mode 100644 index 00000000..b2342571 --- /dev/null +++ b/ext/MCMCDiagnosticToolsStatsBaseExt.jl @@ -0,0 +1,12 @@ +module MCMCDiagnosticToolsStatsBaseExt + +using MCMCDiagnosticTools: MCMCDiagnosticTools +using Statistics: Statistics +using StatsBase: StatsBase + +function MCMCDiagnosticTools._expectand_proxy(::typeof(StatsBase.mad), x) + x_folded = MCMCDiagnosticTools._fold_around_median(x) + return MCMCDiagnosticTools._expectand_proxy(Statistics.median, x_folded) +end + +end diff --git a/src/ess_rhat.jl b/src/ess_rhat.jl index f31eb251..2c26d2a9 100644 --- a/src/ess_rhat.jl +++ b/src/ess_rhat.jl @@ -640,10 +640,6 @@ end function _expectand_proxy(::typeof(Statistics.std), x) return (x .- Statistics.mean(x; dims=_sample_dims(x))) .^ 2 end -function _expectand_proxy(::typeof(StatsBase.mad), x) - x_folded = _fold_around_median(x) - return _expectand_proxy(Statistics.median, x_folded) -end function _expectand_proxy(f::Base.Fix2{typeof(Statistics.quantile),<:Real}, x) y = similar(x) # currently quantile does not support a dims keyword argument From 755914982b706ef0a42e19b9356598cfb868a5bc Mon Sep 17 00:00:00 2001 From: Seth Axen Date: Tue, 29 Sep 2026 00:58:28 +0200 Subject: [PATCH 12/34] chore: Remove StatsBase as package extension --- Project.toml | 1 - src/MCMCDiagnosticTools.jl | 1 - 2 files changed, 2 deletions(-) diff --git a/Project.toml b/Project.toml index 73c9a714..44135c8e 100644 --- a/Project.toml +++ b/Project.toml @@ -14,7 +14,6 @@ MLJModelInterface = "e80e1ace-859a-464e-9ed9-23947d8ae3ea" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b" Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" -StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91" StatsFuns = "4c63d2b9-4356-54db-8cca-17b64c39e42c" Tables = "bd369af6-aec1-5ad0-b16a-f7cc5008161c" diff --git a/src/MCMCDiagnosticTools.jl b/src/MCMCDiagnosticTools.jl index 45653a97..c1cea28e 100644 --- a/src/MCMCDiagnosticTools.jl +++ b/src/MCMCDiagnosticTools.jl @@ -4,7 +4,6 @@ using AbstractFFTs: AbstractFFTs using DataAPI: DataAPI using MLJModelInterface: MLJModelInterface as MMI using SpecialFunctions: SpecialFunctions -using StatsBase: StatsBase using StatsFuns: StatsFuns, sqrt2 using Tables: Tables From 4b8103dbc8ea19ac4e648b2582a54201ffda34d3 Mon Sep 17 00:00:00 2001 From: Seth Axen Date: Tue, 29 Sep 2026 00:59:10 +0200 Subject: [PATCH 13/34] chore: Remove unused DataAPI dependency --- Project.toml | 2 -- src/MCMCDiagnosticTools.jl | 1 - 2 files changed, 3 deletions(-) diff --git a/Project.toml b/Project.toml index 44135c8e..dcc7c6b6 100644 --- a/Project.toml +++ b/Project.toml @@ -8,7 +8,6 @@ projects = ["docs", "test"] [deps] AbstractFFTs = "621f4979-c628-5d54-868e-fcf4e3e8185c" -DataAPI = "9a962f9c-6df0-11e9-0e5d-c546b8b5ee8a" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" MLJModelInterface = "e80e1ace-859a-464e-9ed9-23947d8ae3ea" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" @@ -27,7 +26,6 @@ MCMCDiagnosticToolsStatsBaseExt = "StatsBase" [compat] AbstractFFTs = "0.5, 1" -DataAPI = "1.12" Distributions = "0.25.57" LinearAlgebra = "1.6" MLJModelInterface = "1.11" diff --git a/src/MCMCDiagnosticTools.jl b/src/MCMCDiagnosticTools.jl index c1cea28e..987124a8 100644 --- a/src/MCMCDiagnosticTools.jl +++ b/src/MCMCDiagnosticTools.jl @@ -1,7 +1,6 @@ module MCMCDiagnosticTools using AbstractFFTs: AbstractFFTs -using DataAPI: DataAPI using MLJModelInterface: MLJModelInterface as MMI using SpecialFunctions: SpecialFunctions using StatsFuns: StatsFuns, sqrt2 From 8d469ef2fca37c96659c2b3d22e93a6ae03ac703 Mon Sep 17 00:00:00 2001 From: Seth Axen Date: Tue, 29 Sep 2026 01:37:31 +0200 Subject: [PATCH 14/34] test: Bump EvoTrees test dep lower bound --- test/Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/Project.toml b/test/Project.toml index c3a66f98..8eb8f026 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -26,7 +26,7 @@ MCMCDiagnosticTools = {path = ".."} Aqua = "0.8.14" Distributions = "0.25.57" DynamicHMC = "3.2" -EvoTrees = "0.16.2, 0.18, 0.19" +EvoTrees = "0.16.3, 0.18, 0.19" FFTW = "1.3" LogDensityProblems = "0.12, 1, 2" LogExpFunctions = "0.3.26, 1" From 59b41fa2ba355fce91c8187d665337ea0a47fa7a Mon Sep 17 00:00:00 2001 From: Seth Axen Date: Tue, 29 Sep 2026 01:38:25 +0200 Subject: [PATCH 15/34] test: Test trailing zeros walkback only on StatsBase versions that have it --- test/utils.jl | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/test/utils.jl b/test/utils.jl index 186a12f1..6b382dda 100644 --- a/test/utils.jl +++ b/test/utils.jl @@ -182,8 +182,11 @@ end w = [0.2, 0.5, 0.3, 0.0, 0.0] rng = FixedRNG(nextfloat(1.0)) @test MCMCDiagnosticTools._wsample(rng, w) == 3 - @test MCMCDiagnosticTools._wsample(rng, w) == - StatsBase.wsample(rng, eachindex(w), w) + # StatsBase avoids returning an index with zero weight since v0.34.13 + if pkgversion(StatsBase) ≥ v"0.34.13" + @test MCMCDiagnosticTools._wsample(rng, w) == + StatsBase.wsample(rng, eachindex(w), w) + end end @testset "offset indices" begin From 68b1a3f1b4ccdac7e6174445c972dcba04164002 Mon Sep 17 00:00:00 2001 From: Seth Axen Date: Tue, 29 Sep 2026 01:39:14 +0200 Subject: [PATCH 16/34] chore: Increment minor version number --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index dcc7c6b6..8d8f225d 100644 --- a/Project.toml +++ b/Project.toml @@ -1,6 +1,6 @@ name = "MCMCDiagnosticTools" uuid = "be115224-59cd-429b-ad48-344e309966f0" -version = "0.3.19" +version = "0.4.0" authors = ["David Widmann", "Seth Axen"] [workspace] From 89f09e78cfe460d528cad8f535dafee791c337f6 Mon Sep 17 00:00:00 2001 From: Seth Axen Date: Tue, 29 Sep 2026 02:01:48 +0200 Subject: [PATCH 17/34] Revert "chore: Increment minor version number" This reverts commit 68b1a3f1b4ccdac7e6174445c972dcba04164002. --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 8d8f225d..dcc7c6b6 100644 --- a/Project.toml +++ b/Project.toml @@ -1,6 +1,6 @@ name = "MCMCDiagnosticTools" uuid = "be115224-59cd-429b-ad48-344e309966f0" -version = "0.4.0" +version = "0.3.19" authors = ["David Widmann", "Seth Axen"] [workspace] From 8291310e5fe2942efbeed1622bb8a3397c2dcf09 Mon Sep 17 00:00:00 2001 From: Seth Axen Date: Tue, 29 Sep 2026 02:15:45 +0200 Subject: [PATCH 18/34] chore: Bump test and doc compat to match package --- docs/Project.toml | 2 +- test/Project.toml | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/docs/Project.toml b/docs/Project.toml index 893ee69d..e5ec787a 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -16,7 +16,7 @@ MCMCDiagnosticTools = {path = ".."} Documenter = "1" DocumenterInterLinks = "1" EvoTrees = "0.18.2, 0.19" -MCMCDiagnosticTools = "0.3.15" +MCMCDiagnosticTools = "0.4" MLJBase = "1.10" MLJIteration = "0.6" StatisticalMeasures = "0.3" diff --git a/test/Project.toml b/test/Project.toml index 8eb8f026..d02f13e9 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -30,7 +30,7 @@ EvoTrees = "0.16.3, 0.18, 0.19" FFTW = "1.3" LogDensityProblems = "0.12, 1, 2" LogExpFunctions = "0.3.26, 1" -MCMCDiagnosticTools = "0.3.15" +MCMCDiagnosticTools = "0.4" MLJBase = "0.20.18, 0.21, 1" MLJDecisionTreeInterface = "0.3, 0.4, 0.5" MLJLIBSVMInterface = "0.2" From 227972779534236b6f56bdb5d9ea7e5be511cfbe Mon Sep 17 00:00:00 2001 From: Seth Axen Date: Tue, 29 Sep 2026 02:18:02 +0200 Subject: [PATCH 19/34] Reapply "chore: Increment minor version number" This reverts commit 89f09e78cfe460d528cad8f535dafee791c337f6. --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index dcc7c6b6..8d8f225d 100644 --- a/Project.toml +++ b/Project.toml @@ -1,6 +1,6 @@ name = "MCMCDiagnosticTools" uuid = "be115224-59cd-429b-ad48-344e309966f0" -version = "0.3.19" +version = "0.4.0" authors = ["David Widmann", "Seth Axen"] [workspace] From ab44c02fc7250c8548b72065028dd03883b93c1f Mon Sep 17 00:00:00 2001 From: Seth Axen Date: Tue, 29 Sep 2026 02:19:07 +0200 Subject: [PATCH 20/34] test: Add test for rstar error hint --- test/rstar.jl | 36 ++++++++++++++++++++++++++++++++++++ 1 file changed, 36 insertions(+) diff --git a/test/rstar.jl b/test/rstar.jl index 763dfe31..9b9a12a5 100644 --- a/test/rstar.jl +++ b/test/rstar.jl @@ -189,6 +189,42 @@ end @test_throws ArgumentError MCMCDiagnosticTools._rstar(1.0, rand(2), rand(2)) end + @testset "probabilistic predictions without Distributions" begin + # Distributions is loaded here, so check in a new process that only loads + # MCMCDiagnosticTools + code = """ + using MCMCDiagnosticTools + const MMI = MCMCDiagnosticTools.MMI + # `MMI.classes` otherwise requires MLJBase, which loads Distributions + MMI.classes(y::Vector{Int}) = sort!(unique(y)) + T = AbstractVector{MMI.Density{MMI.Finite}} + err = try + MCMCDiagnosticTools._rstar(T, [0.3, 0.8], [1, 2]) + nothing + catch e + e + end + println( + Base.get_extension(MCMCDiagnosticTools, :MCMCDiagnosticToolsDistributionsExt) === + nothing, + ) + println(err isa MethodError && err.f === MCMCDiagnosticTools._rstar_distribution) + err === nothing || showerror(stdout, err) + """ + cmd = `$(Base.julia_cmd()) --startup-file=no --project=$(Base.active_project()) -e $code` + out = read(cmd, String) + extension_not_loaded, method_error_raised, msg = split(out, '\n'; limit=3) + @test extension_not_loaded == "true" + @test method_error_raised == "true" + @test occursin("requires Distributions to be loaded", msg) + + # other `MethodError`s get no hint + @test !occursin( + "requires Distributions to be loaded", + sprint(showerror, MethodError(rstar, (1,))), + ) + end + @testset "single chain: method ambiguity issue" begin samples = rand(1:5, N) rng = MersenneTwister(42) From 3305b2fa763c28317b5f0ac4c8c1e1e0829d00c6 Mon Sep 17 00:00:00 2001 From: Seth Axen Date: Tue, 29 Sep 2026 13:42:48 +0200 Subject: [PATCH 21/34] perf: Use alias table for DAR(1) sampling in discretediag --- Project.toml | 2 ++ src/MCMCDiagnosticTools.jl | 1 + src/discretediag.jl | 15 ++++++++------- 3 files changed, 11 insertions(+), 7 deletions(-) diff --git a/Project.toml b/Project.toml index 8d8f225d..8d1d70f5 100644 --- a/Project.toml +++ b/Project.toml @@ -8,6 +8,7 @@ projects = ["docs", "test"] [deps] AbstractFFTs = "621f4979-c628-5d54-868e-fcf4e3e8185c" +AliasTables = "66dad0bd-aa9a-41b7-9441-69ab47430ed8" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" MLJModelInterface = "e80e1ace-859a-464e-9ed9-23947d8ae3ea" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" @@ -26,6 +27,7 @@ MCMCDiagnosticToolsStatsBaseExt = "StatsBase" [compat] AbstractFFTs = "0.5, 1" +AliasTables = "1" Distributions = "0.25.57" LinearAlgebra = "1.6" MLJModelInterface = "1.11" diff --git a/src/MCMCDiagnosticTools.jl b/src/MCMCDiagnosticTools.jl index 987124a8..d2d91673 100644 --- a/src/MCMCDiagnosticTools.jl +++ b/src/MCMCDiagnosticTools.jl @@ -1,6 +1,7 @@ module MCMCDiagnosticTools using AbstractFFTs: AbstractFFTs +using AliasTables: AliasTables using MLJModelInterface: MLJModelInterface as MMI using SpecialFunctions: SpecialFunctions using StatsFuns: StatsFuns, sqrt2 diff --git a/src/discretediag.jl b/src/discretediag.jl index 01feccda..e3056914 100644 --- a/src/discretediag.jl +++ b/src/discretediag.jl @@ -185,7 +185,7 @@ function bd_inner(Y::AbstractMatrix, m::Int) end @doc raw""" - simulate_DAR1!(X::Matrix{Int}, phi::Float64, prob::AbstractVector{<:Real}) + simulate_DAR1!(X::Matrix{Int}, phi::Float64, sampler) Simulate a DAR(1) model independently in each column of `X`. @@ -196,14 +196,14 @@ X_t = \alpha_t X_{t-1} + (1 - \alpha_t) \epsilon_{t-1}, where ```math \begin{aligned} -X_1 \sim \operatorname{Categorical}(\mathrm{prob}), \\ +X_1 \sim \text{sampler}, \\ \alpha_t \sim \operatorname{Bernoulli}(\phi), \\ -\epsilon_{t-1} \sim \operatorname{Categorical}(\mathrm{prob}), +\epsilon_{t-1} \sim \text{sampler}, \end{aligned} ``` are independent random variables. """ -function simulate_DAR1!(X::Matrix{Int}, phi::Float64, prob::AbstractVector{<:Real}) +function simulate_DAR1!(X::Matrix{Int}, phi::Float64, sampler) n = size(X, 1) n > 0 || error("output matrix must be non-empty") rng = Random.default_rng() @@ -211,7 +211,7 @@ function simulate_DAR1!(X::Matrix{Int}, phi::Float64, prob::AbstractVector{<:Rea # for all simulations @inbounds for j in axes(X, 2) # sample first value from categorical distribution with probabilities `prob` - X[1, j] = _wsample(rng, prob) + X[1, j] = rand(rng, sampler) for t in 2:n # compute next value @@ -220,7 +220,7 @@ function simulate_DAR1!(X::Matrix{Int}, phi::Float64, prob::AbstractVector{<:Rea X[t - 1, j] else # sample value with probability `1-phi` - _wsample(rng, prob) + rand(rng, sampler) end end end @@ -316,10 +316,11 @@ function diag_all( end elseif method == :DARBOOT stat = t * sum(chi_stat) + sampler_phat = AliasTables.AliasTable(phat) bstats = zeros(nsim) Y = Matrix{Int}(undef, t, d) for b in 1:nsim - simulate_DAR1!(Y, phia, phat) + simulate_DAR1!(Y, phia, sampler_phat) s = hangartner_inner(Y, m)[1] @inbounds bstats[b] = s end From 8c04c722168a0cdbe817a15007ff44b70e144d20 Mon Sep 17 00:00:00 2001 From: Seth Axen Date: Tue, 29 Sep 2026 13:42:55 +0200 Subject: [PATCH 22/34] perf: Use alias tables for Markov chain sampling in discretediag --- src/discretediag.jl | 23 ++++++++++---- src/utils.jl | 17 ---------- test/discretediag.jl | 17 ++++++++++ test/utils.jl | 76 -------------------------------------------- 4 files changed, 34 insertions(+), 99 deletions(-) diff --git a/src/discretediag.jl b/src/discretediag.jl index e3056914..034927ae 100644 --- a/src/discretediag.jl +++ b/src/discretediag.jl @@ -228,13 +228,22 @@ function simulate_DAR1!(X::Matrix{Int}, phi::Float64, sampler) return X end -function simulate_MC(N::Int, P::Matrix{Float64}) +# Construct a sampler of the next state of a Markov chain for each current state from the +# rows of the transition matrix `P`. States without transition probabilities (i.e. that +# were never observed to transition) always transition to the first state. +function transition_samplers(P::Matrix{Float64}) + return map(eachrow(P)) do p + all(iszero, p) && return AliasTables.AliasTable(Float64.(eachindex(p) .== 1)) + return AliasTables.AliasTable(p) + end +end + +function simulate_MC(N::Int, samplers::AbstractVector{<:AliasTables.AliasTable}) X = zeros(Int, N) - n, m = size(P) rng = Random.default_rng() - X[1] = rand(rng, 1:n) + X[1] = rand(rng, 1:length(samplers)) for i in 2:N - X[i] = _wsample(rng, vec(P[X[i - 1], :])) + X[i] = rand(rng, samplers[X[i - 1]]) end return X end @@ -328,9 +337,10 @@ function diag_all( df0 = Statistics.mean(non_nan_bstats) pval = Statistics.mean(stat <= x for x in non_nan_bstats) elseif method == :MCBOOT + samplers_mP = transition_samplers(mP) bstats = zeros(Float64, nsim) for b in 1:nsim - Y = reduce(hcat, [simulate_MC(t, mP) for j in 1:d]) + Y = reduce(hcat, [simulate_MC(t, samplers_mP) for j in 1:d]) s = hangartner_inner(Y, m)[1] bstats[b] = s end @@ -345,9 +355,10 @@ function diag_all( end elseif method == :billingsleyBOOT stat = hot_stat + samplers_mP = transition_samplers(mP) bstats = zeros(Float64, nsim) for b in 1:nsim - Y = reduce(hcat, [simulate_MC(t, mP) for j in 1:d]) + Y = reduce(hcat, [simulate_MC(t, samplers_mP) for j in 1:d]) (s, sd) = bd_inner(Y, m)[1:2] bstats[b] = s / sd end diff --git a/src/utils.jl b/src/utils.jl index 2c6ac9a8..e513b5eb 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -208,23 +208,6 @@ function _tiedrank!(ranks::AbstractVector, x::AbstractVector) return ranks end -# sample an index `i` of `weights` with probability proportional to `weights[i]`. -# for `Float64` weights, equivalent to `StatsBase.wsample(rng, eachindex(weights), weights)` -function _wsample(rng::Random.AbstractRNG, weights::AbstractVector{T}) where {T<:Real} - t = rand(rng, float(T)) * sum(weights) - i = firstindex(weights) - ilast = lastindex(weights) - cw = weights[i] - while true - cw ≥ t && return i - # may happen with floating point weights due to numerical inaccuracies. - # never return an index with zero weight. - i == ilast && return something(findlast(!iszero, weights), i) - i += 1 - cw += weights[i] - end -end - # count the number of occurrences of each value of `levels` in `x`. # equivalent to `StatsBase.counts(x, levels)` function _counts(x::AbstractArray{<:Integer}, levels::UnitRange{<:Integer}) diff --git a/test/discretediag.jl b/test/discretediag.jl index 2796aa96..65f5ba36 100644 --- a/test/discretediag.jl +++ b/test/discretediag.jl @@ -31,4 +31,21 @@ @test_throws ArgumentError discretediag(samples; frac=x) end end + + @testset "transition_samplers" begin + # state 2 was never observed to transition + P = [0.2 0.8 0.0; 0.0 0.0 0.0; 0.5 0.0 0.5] + samplers = MCMCDiagnosticTools.transition_samplers(P) + @test length(samplers) == 3 + rng = Xoshiro(42) + ndraws = 100_000 + for i in (1, 3) + counts = zeros(Int, 3) + for _ in 1:ndraws + counts[rand(rng, samplers[i])] += 1 + end + @test counts ./ ndraws ≈ P[i, :] atol = 0.01 + end + @test all(rand(rng, samplers[2]) == 1 for _ in 1:100) + end end diff --git a/test/utils.jl b/test/utils.jl index 6b382dda..961d67ee 100644 --- a/test/utils.jl +++ b/test/utils.jl @@ -153,82 +153,6 @@ end end end -# RNG that always generates the same `Float64`, used to control `_wsample` -struct FixedRNG <: Random.AbstractRNG - x::Float64 -end -function Random.rand(rng::FixedRNG, ::Random.SamplerTrivial{Random.CloseOpen01{Float64}}) - return rng.x -end - -@testset "_wsample" begin - @testset "basic" begin - w = [0.2, 0.0, 0.5, 0.3] - # the first index whose cumulative weight is at least `rand() * sum(w)` is sampled - @test MCMCDiagnosticTools._wsample(FixedRNG(0.0), w) == 1 - @test MCMCDiagnosticTools._wsample(FixedRNG(0.1), w) == 1 - @test MCMCDiagnosticTools._wsample(FixedRNG(0.3), w) == 3 - @test MCMCDiagnosticTools._wsample(FixedRNG(0.9), w) == 4 - @test @inferred(MCMCDiagnosticTools._wsample(Random.default_rng(), w)) isa Int - end - - @testset "all-zero weights return the first index" begin - @test MCMCDiagnosticTools._wsample(FixedRNG(0.5), zeros(3)) == 1 - end - - @testset "never returns an index with zero weight" begin - # emulate floating-point error causing `rand() * sum(w)` to exceed the cumulative - # sum of the weights - w = [0.2, 0.5, 0.3, 0.0, 0.0] - rng = FixedRNG(nextfloat(1.0)) - @test MCMCDiagnosticTools._wsample(rng, w) == 3 - # StatsBase avoids returning an index with zero weight since v0.34.13 - if pkgversion(StatsBase) ≥ v"0.34.13" - @test MCMCDiagnosticTools._wsample(rng, w) == - StatsBase.wsample(rng, eachindex(w), w) - end - end - - @testset "offset indices" begin - w = OffsetArray([0.2, 0.0, 0.5, 0.3], -2) - @test MCMCDiagnosticTools._wsample(FixedRNG(0.0), w) == -1 - @test MCMCDiagnosticTools._wsample(FixedRNG(0.3), w) == 1 - end - - @testset "consistent with StatsBase.wsample" begin - p = rand(10) - p ./= sum(p) - weights = ( - [1.0], - rand(5), - [rand(5); zeros(3)], - [0.0; rand(3); 0.0; rand(2)], - p, - rand(0:3, 8), - ) - for w in weights - seed = rand(UInt) - rng1, rng2 = Xoshiro(seed), Xoshiro(seed) - @test all(1:100) do _ - return MCMCDiagnosticTools._wsample(rng1, w) == - StatsBase.wsample(rng2, eachindex(w), w) - end - end - end - - @testset "sampling frequencies" begin - w = [1.0, 0.0, 3.0, 6.0] - rng = Xoshiro(42) - ndraws = 100_000 - counts = zeros(Int, length(w)) - for _ in 1:ndraws - counts[MCMCDiagnosticTools._wsample(rng, w)] += 1 - end - @test counts[2] == 0 - @test counts ./ ndraws ≈ w ./ sum(w) atol = 0.01 - end -end - @testset "_counts" begin x = [0, 3, 3, 7, 8, -1, 3, 0] @test @inferred(MCMCDiagnosticTools._counts(x, 0:7)) == [2, 0, 0, 3, 0, 0, 0, 1] From 8f502f699ab2bdf3b3737861eb552b24e6882cb7 Mon Sep 17 00:00:00 2001 From: Seth Axen Date: Tue, 29 Sep 2026 14:57:35 +0200 Subject: [PATCH 23/34] perf: Replace broadcast over views with for loop --- src/utils.jl | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/utils.jl b/src/utils.jl index e513b5eb..c48cd27e 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -202,7 +202,9 @@ function _tiedrank!(ranks::AbstractVector, x::AbstractVector) stop += 1 end rank = (start + stop) / 2 + offset - @views ranks[perm[start:stop]] .= rank + for i in start:stop + ranks[perm[i]] = rank + end start = stop + 1 end return ranks From 6eadc190ac004f404cbcfcb5681b7157ad5b2f8e Mon Sep 17 00:00:00 2001 From: Seth Axen Date: Tue, 29 Sep 2026 16:44:35 +0200 Subject: [PATCH 24/34] refactor: Move rstar implementation to MLJBase extension --- Project.toml | 8 +- ext/MCMCDiagnosticToolsDistributionsExt.jl | 16 --- ext/MCMCDiagnosticToolsMLJBaseExt.jl | 160 +++++++++++++++++++++ src/MCMCDiagnosticTools.jl | 2 - src/rstar.jl | 150 ++----------------- test/rstar.jl | 25 ++-- 6 files changed, 187 insertions(+), 174 deletions(-) delete mode 100644 ext/MCMCDiagnosticToolsDistributionsExt.jl create mode 100644 ext/MCMCDiagnosticToolsMLJBaseExt.jl diff --git a/Project.toml b/Project.toml index 8d1d70f5..a3255395 100644 --- a/Project.toml +++ b/Project.toml @@ -10,19 +10,20 @@ projects = ["docs", "test"] AbstractFFTs = "621f4979-c628-5d54-868e-fcf4e3e8185c" AliasTables = "66dad0bd-aa9a-41b7-9441-69ab47430ed8" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" -MLJModelInterface = "e80e1ace-859a-464e-9ed9-23947d8ae3ea" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b" Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" StatsFuns = "4c63d2b9-4356-54db-8cca-17b64c39e42c" -Tables = "bd369af6-aec1-5ad0-b16a-f7cc5008161c" [weakdeps] Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" +MLJBase = "a7f614a8-145f-11e9-1d2a-a57a1082229d" +MLJModelInterface = "e80e1ace-859a-464e-9ed9-23947d8ae3ea" StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91" +Tables = "bd369af6-aec1-5ad0-b16a-f7cc5008161c" [extensions] -MCMCDiagnosticToolsDistributionsExt = "Distributions" +MCMCDiagnosticToolsMLJBaseExt = ["Distributions", "MLJBase", "MLJModelInterface", "Tables"] MCMCDiagnosticToolsStatsBaseExt = "StatsBase" [compat] @@ -30,6 +31,7 @@ AbstractFFTs = "0.5, 1" AliasTables = "1" Distributions = "0.25.57" LinearAlgebra = "1.6" +MLJBase = "0.20.18, 0.21, 1" MLJModelInterface = "1.11" Random = "1.6" SpecialFunctions = "1.8.4, 2" diff --git a/ext/MCMCDiagnosticToolsDistributionsExt.jl b/ext/MCMCDiagnosticToolsDistributionsExt.jl deleted file mode 100644 index 5e38cf20..00000000 --- a/ext/MCMCDiagnosticToolsDistributionsExt.jl +++ /dev/null @@ -1,16 +0,0 @@ -module MCMCDiagnosticToolsDistributionsExt - -using MCMCDiagnosticTools: MCMCDiagnosticTools -using Distributions: Distributions - -function MCMCDiagnosticTools._rstar_distribution(predictions, ytest, nclasses::Int) - # create Poisson binomial distribution with support `0:length(predictions)` - distribution = Distributions.PoissonBinomial(map(Distributions.pdf, predictions, ytest)) - - # scale distribution to support in `[0, nclasses]` - scaled_distribution = (nclasses//length(predictions)) * distribution - - return scaled_distribution -end - -end diff --git a/ext/MCMCDiagnosticToolsMLJBaseExt.jl b/ext/MCMCDiagnosticToolsMLJBaseExt.jl new file mode 100644 index 00000000..89a8e98f --- /dev/null +++ b/ext/MCMCDiagnosticToolsMLJBaseExt.jl @@ -0,0 +1,160 @@ +module MCMCDiagnosticToolsMLJBaseExt + +using Distributions: Distributions +using MCMCDiagnosticTools: + MCMCDiagnosticTools, shuffle_split_stratified, split_chain_indices +using MLJModelInterface: MLJModelInterface as MMI +using Random: Random +using Statistics: Statistics +using Tables: Tables + +import MCMCDiagnosticTools: + _check_model_supports_continuous_inputs, + _check_model_supports_multiclass_predictions, + _check_model_supports_multiclass_targets, + _rstar + +function _rstar( + rng::Random.AbstractRNG, + classifier, + x, + y::AbstractVector{Int}; + subset::Real, + split_chains::Int, + verbosity::Int, +) + # check the arguments + _check_model_supports_continuous_inputs(classifier) + _check_model_supports_multiclass_targets(classifier) + _check_model_supports_multiclass_predictions(classifier) + MMI.nrows(x) != length(y) && throw(DimensionMismatch()) + 0 < subset < 1 || throw(ArgumentError("`subset` must be a number in (0, 1)")) + + # randomly sub-select training and testing set + ysplit = split_chain_indices(y, split_chains) + train_ids, test_ids = shuffle_split_stratified(rng, ysplit, subset) + 0 < length(train_ids) < length(y) || + throw(ArgumentError("training and test data subsets must not be empty")) + + xtable = _astable(x) + ycategorical = MMI.categorical(ysplit) + + # train classifier on training data + data = MMI.reformat(classifier, xtable, ycategorical) + train_data = MMI.selectrows(classifier, train_ids, data...) + fitresult, _ = MMI.fit(classifier, verbosity, train_data...) + + # compute predictions on test data + # we exploit that MLJ demands that + # reformat(model, args...)[1] = reformat(model, args[1]) + # (https://alan-turing-institute.github.io/MLJ.jl/dev/adding_models_for_general_use/#Implementing-a-data-front-end) + test_data = MMI.selectrows(classifier, test_ids, data[1]) + predictions = _predict(classifier, fitresult, test_data...) + + # compute statistic + ytest = ycategorical[test_ids] + result = _rstar(MMI.scitype(predictions), predictions, ytest) + + return result +end + +# check that the model supports the inputs and targets, and has predictions of the desired form +function _check_model_supports_continuous_inputs(classifier) + # ideally we would not allow MMI.Unknown but some models do not implement the traits + input_scitype_classifier = MMI.input_scitype(classifier) + if input_scitype_classifier !== MMI.Unknown && + !(MMI.Table(MMI.Continuous) <: input_scitype_classifier) + throw( + ArgumentError( + "classifier does not support tables of continuous values as inputs" + ), + ) + end + return nothing +end +function _check_model_supports_multiclass_targets(classifier) + target_scitype_classifier = MMI.target_scitype(classifier) + if target_scitype_classifier !== MMI.Unknown && + !(AbstractVector{<:MMI.Finite} <: target_scitype_classifier) + throw( + ArgumentError( + "classifier does not support vectors of multi-class labels as targets" + ), + ) + end + return nothing +end +function _check_model_supports_multiclass_predictions(classifier) + if !( + MMI.predict_scitype(classifier) <: Union{ + MMI.Unknown, + AbstractVector{<:MMI.Finite}, + AbstractVector{<:MMI.Density{<:MMI.Finite}}, + } + ) + throw( + ArgumentError( + "classifier does not support vectors of multi-class labels or their densities as predictions", + ), + ) + end + return nothing +end + +_astable(x::AbstractVecOrMat) = Tables.table(x) +_astable(x) = Tables.istable(x) ? x : throw(ArgumentError("Argument is not a valid table")) + +# Workaround for https://github.com/JuliaAI/MLJBase.jl/issues/863 +# `MLJModelInterface.predict` sometimes returns predictions and sometimes predictions + additional information +# TODO: Remove once the upstream issue is fixed +function _predict(model::MMI.Model, fitresult, x) + y = MMI.predict(model, fitresult, x) + return if :predict in MMI.reporting_operations(model) + first(y) + else + y + end +end + +# R⋆ for deterministic predictions (algorithm 1) +function _rstar( + ::Type{<:AbstractVector{<:MMI.Finite}}, + predictions::AbstractVector, + ytest::AbstractVector, +) + length(predictions) == length(ytest) || + error("numbers of predictions and targets must be equal") + mean_accuracy = Statistics.mean(p == y for (p, y) in zip(predictions, ytest)) + nclasses = length(MMI.classes(ytest)) + return nclasses * mean_accuracy +end + +# R⋆ for probabilistic predictions (algorithm 2) +function _rstar( + ::Type{<:AbstractVector{<:MMI.Density{<:MMI.Finite}}}, + predictions::AbstractVector, + ytest::AbstractVector, +) + length(predictions) == length(ytest) || + error("numbers of predictions and targets must be equal") + + # create Poisson binomial distribution with support `0:length(predictions)` + distribution = Distributions.PoissonBinomial(map(Distributions.pdf, predictions, ytest)) + + # scale distribution to support in `[0, nclasses]` + nclasses = length(MMI.classes(ytest)) + scaled_distribution = (nclasses//length(predictions)) * distribution + + return scaled_distribution +end + +# unsupported types of predictions and targets +function _rstar(::Any, predictions, targets) + throw( + ArgumentError( + "unsupported types of predictions ($(typeof(predictions))) and targets ($(typeof(targets)))", + ), + ) +end + +end diff --git a/src/MCMCDiagnosticTools.jl b/src/MCMCDiagnosticTools.jl index d2d91673..ca8037bf 100644 --- a/src/MCMCDiagnosticTools.jl +++ b/src/MCMCDiagnosticTools.jl @@ -2,10 +2,8 @@ module MCMCDiagnosticTools using AbstractFFTs: AbstractFFTs using AliasTables: AliasTables -using MLJModelInterface: MLJModelInterface as MMI using SpecialFunctions: SpecialFunctions using StatsFuns: StatsFuns, sqrt2 -using Tables: Tables using LinearAlgebra: LinearAlgebra using Random: Random diff --git a/src/rstar.jl b/src/rstar.jl index 2a1f46f1..118d678e 100644 --- a/src/rstar.jl +++ b/src/rstar.jl @@ -28,97 +28,7 @@ function rstar( split_chains::Int=2, verbosity::Int=0, ) - # check the arguments - _check_model_supports_continuous_inputs(classifier) - _check_model_supports_multiclass_targets(classifier) - _check_model_supports_multiclass_predictions(classifier) - MMI.nrows(x) != length(y) && throw(DimensionMismatch()) - 0 < subset < 1 || throw(ArgumentError("`subset` must be a number in (0, 1)")) - - # randomly sub-select training and testing set - ysplit = split_chain_indices(y, split_chains) - train_ids, test_ids = shuffle_split_stratified(rng, ysplit, subset) - 0 < length(train_ids) < length(y) || - throw(ArgumentError("training and test data subsets must not be empty")) - - xtable = _astable(x) - ycategorical = MMI.categorical(ysplit) - - # train classifier on training data - data = MMI.reformat(classifier, xtable, ycategorical) - train_data = MMI.selectrows(classifier, train_ids, data...) - fitresult, _ = MMI.fit(classifier, verbosity, train_data...) - - # compute predictions on test data - # we exploit that MLJ demands that - # reformat(model, args...)[1] = reformat(model, args[1]) - # (https://alan-turing-institute.github.io/MLJ.jl/dev/adding_models_for_general_use/#Implementing-a-data-front-end) - test_data = MMI.selectrows(classifier, test_ids, data[1]) - predictions = _predict(classifier, fitresult, test_data...) - - # compute statistic - ytest = ycategorical[test_ids] - result = _rstar(MMI.scitype(predictions), predictions, ytest) - - return result -end - -# check that the model supports the inputs and targets, and has predictions of the desired form -function _check_model_supports_continuous_inputs(classifier) - # ideally we would not allow MMI.Unknown but some models do not implement the traits - input_scitype_classifier = MMI.input_scitype(classifier) - if input_scitype_classifier !== MMI.Unknown && - !(MMI.Table(MMI.Continuous) <: input_scitype_classifier) - throw( - ArgumentError( - "classifier does not support tables of continuous values as inputs" - ), - ) - end - return nothing -end -function _check_model_supports_multiclass_targets(classifier) - target_scitype_classifier = MMI.target_scitype(classifier) - if target_scitype_classifier !== MMI.Unknown && - !(AbstractVector{<:MMI.Finite} <: target_scitype_classifier) - throw( - ArgumentError( - "classifier does not support vectors of multi-class labels as targets" - ), - ) - end - return nothing -end -function _check_model_supports_multiclass_predictions(classifier) - if !( - MMI.predict_scitype(classifier) <: Union{ - MMI.Unknown, - AbstractVector{<:MMI.Finite}, - AbstractVector{<:MMI.Density{<:MMI.Finite}}, - } - ) - throw( - ArgumentError( - "classifier does not support vectors of multi-class labels or their densities as predictions", - ), - ) - end - return nothing -end - -_astable(x::AbstractVecOrMat) = Tables.table(x) -_astable(x) = Tables.istable(x) ? x : throw(ArgumentError("Argument is not a valid table")) - -# Workaround for https://github.com/JuliaAI/MLJBase.jl/issues/863 -# `MLJModelInterface.predict` sometimes returns predictions and sometimes predictions + additional information -# TODO: Remove once the upstream issue is fixed -function _predict(model::MMI.Model, fitresult, x) - y = MMI.predict(model, fitresult, x) - return if :predict in MMI.reporting_operations(model) - first(y) - else - y - end + return _rstar(rng, classifier, x, y; subset, split_chains, verbosity) end """ @@ -147,8 +57,11 @@ within-chain convergence. The training of the classifier can be inspected by adj If the classifier is deterministic, i.e., if it predicts a class, the value of the ``R^*`` statistic is returned (algorithm 1). If the classifier is probabilistic, i.e., if it outputs probabilities of classes, the scaled Poisson-binomial distribution of the ``R^*`` statistic -is returned (algorithm 2). Probabilistic classifiers require -[Distributions.jl](https://github.com/JuliaStats/Distributions.jl) to be loaded. +is returned (algorithm 2). + +!!! info + `rstar` requires [MLJBase.jl](https://github.com/JuliaAI/MLJBase.jl) to be loaded, e.g. + with `using MLJBase` or `using MLJ`. !!! note The correctness of the statistic depends on the convergence of the `classifier` used @@ -165,7 +78,6 @@ julia> samples = fill(4.0, 100, 3, 2); One can compute the distribution of the ``R^*`` statistic (algorithm 2) with a probabilistic classifier. For instance, we can use a gradient-boosted trees model with `nrounds = 100` sequentially stacked trees and learning rate `eta = 0.05`: -Because MLJBase already loads Distributions, we don't need to do so ourselves. ```jldoctest rstar julia> model = EvoTreeClassifier(; nrounds=100, eta=0.05); @@ -234,51 +146,15 @@ function rstar(classifier, x::AbstractArray; kwargs...) return rstar(Random.default_rng(), classifier, x; kwargs...) end -# R⋆ for deterministic predictions (algorithm 1) -function _rstar( - ::Type{<:AbstractVector{<:MMI.Finite}}, - predictions::AbstractVector, - ytest::AbstractVector, -) - length(predictions) == length(ytest) || - error("numbers of predictions and targets must be equal") - mean_accuracy = Statistics.mean(p == y for (p, y) in zip(predictions, ytest)) - nclasses = length(MMI.classes(ytest)) - return nclasses * mean_accuracy -end - -# R⋆ for probabilistic predictions (algorithm 2) -function _rstar( - ::Type{<:AbstractVector{<:MMI.Density{<:MMI.Finite}}}, - predictions::AbstractVector, - ytest::AbstractVector, -) - length(predictions) == length(ytest) || - error("numbers of predictions and targets must be equal") - - nclasses = length(MMI.classes(ytest)) - return _rstar_distribution(predictions, ytest, nclasses) -end - -# Implemented in MCMCDiagnosticToolsDistributionsExt -function _rstar_distribution end - function _rstar_error_hint(io, exc, argtypes, kwargs) - if exc.f === _rstar_distribution - print( - io, - "\n`rstar` with a probabilistic classifier requires Distributions to be loaded, " * - "e.g. with `using Distributions`.", - ) + if exc.f === _rstar + print(io, "\n`rstar` requires MLJBase to be loaded, e.g. with `using MLJBase`.") end return nothing end -# unsupported types of predictions and targets -function _rstar(::Any, predictions, targets) - throw( - ArgumentError( - "unsupported types of predictions ($(typeof(predictions))) and targets ($(typeof(targets)))", - ), - ) -end +# Implemented in MCMCDiagnosticToolsMLJBaseExt +function _rstar end +function _check_model_supports_continuous_inputs end +function _check_model_supports_multiclass_targets end +function _check_model_supports_multiclass_predictions end diff --git a/test/rstar.jl b/test/rstar.jl index 9b9a12a5..3a89d38b 100644 --- a/test/rstar.jl +++ b/test/rstar.jl @@ -189,26 +189,19 @@ end @test_throws ArgumentError MCMCDiagnosticTools._rstar(1.0, rand(2), rand(2)) end - @testset "probabilistic predictions without Distributions" begin - # Distributions is loaded here, so check in a new process that only loads + @testset "rstar without MLJBase" begin + # MLJBase is loaded here, so check in a new process that only loads # MCMCDiagnosticTools code = """ using MCMCDiagnosticTools - const MMI = MCMCDiagnosticTools.MMI - # `MMI.classes` otherwise requires MLJBase, which loads Distributions - MMI.classes(y::Vector{Int}) = sort!(unique(y)) - T = AbstractVector{MMI.Density{MMI.Finite}} err = try - MCMCDiagnosticTools._rstar(T, [0.3, 0.8], [1, 2]) + rstar(nothing, randn(10, 2)) nothing catch e e end - println( - Base.get_extension(MCMCDiagnosticTools, :MCMCDiagnosticToolsDistributionsExt) === - nothing, - ) - println(err isa MethodError && err.f === MCMCDiagnosticTools._rstar_distribution) + println(Base.get_extension(MCMCDiagnosticTools, :MCMCDiagnosticToolsMLJBaseExt) === nothing) + println(err isa MethodError) err === nothing || showerror(stdout, err) """ cmd = `$(Base.julia_cmd()) --startup-file=no --project=$(Base.active_project()) -e $code` @@ -216,12 +209,12 @@ end extension_not_loaded, method_error_raised, msg = split(out, '\n'; limit=3) @test extension_not_loaded == "true" @test method_error_raised == "true" - @test occursin("requires Distributions to be loaded", msg) + @test occursin("no method matching _rstar(", msg) + @test occursin("requires MLJBase to be loaded", msg) - # other `MethodError`s get no hint + # with MLJBase loaded, other `MethodError`s get no hint @test !occursin( - "requires Distributions to be loaded", - sprint(showerror, MethodError(rstar, (1,))), + "requires MLJBase to be loaded", sprint(showerror, MethodError(rstar, (1,))) ) end From 12a17795801454fc116eb05f2ae5464f792031f0 Mon Sep 17 00:00:00 2001 From: Seth Axen Date: Tue, 29 Sep 2026 16:49:46 +0200 Subject: [PATCH 25/34] chore: Use patch instead of minor version bump --- Project.toml | 2 +- docs/Project.toml | 2 +- test/Project.toml | 2 +- 3 files changed, 3 insertions(+), 3 deletions(-) diff --git a/Project.toml b/Project.toml index a3255395..0111c367 100644 --- a/Project.toml +++ b/Project.toml @@ -1,6 +1,6 @@ name = "MCMCDiagnosticTools" uuid = "be115224-59cd-429b-ad48-344e309966f0" -version = "0.4.0" +version = "0.3.20" authors = ["David Widmann", "Seth Axen"] [workspace] diff --git a/docs/Project.toml b/docs/Project.toml index e5ec787a..893ee69d 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -16,7 +16,7 @@ MCMCDiagnosticTools = {path = ".."} Documenter = "1" DocumenterInterLinks = "1" EvoTrees = "0.18.2, 0.19" -MCMCDiagnosticTools = "0.4" +MCMCDiagnosticTools = "0.3.15" MLJBase = "1.10" MLJIteration = "0.6" StatisticalMeasures = "0.3" diff --git a/test/Project.toml b/test/Project.toml index d02f13e9..8eb8f026 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -30,7 +30,7 @@ EvoTrees = "0.16.3, 0.18, 0.19" FFTW = "1.3" LogDensityProblems = "0.12, 1, 2" LogExpFunctions = "0.3.26, 1" -MCMCDiagnosticTools = "0.4" +MCMCDiagnosticTools = "0.3.15" MLJBase = "0.20.18, 0.21, 1" MLJDecisionTreeInterface = "0.3, 0.4, 0.5" MLJLIBSVMInterface = "0.2" From 2d8d4c0034c2742c4f3ab0852a8f43377651a09e Mon Sep 17 00:00:00 2001 From: Seth Axen Date: Tue, 29 Sep 2026 17:08:12 +0200 Subject: [PATCH 26/34] style: Apply formatting suggestions Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> --- ext/MCMCDiagnosticToolsMLJBaseExt.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ext/MCMCDiagnosticToolsMLJBaseExt.jl b/ext/MCMCDiagnosticToolsMLJBaseExt.jl index 89a8e98f..bbfb5281 100644 --- a/ext/MCMCDiagnosticToolsMLJBaseExt.jl +++ b/ext/MCMCDiagnosticToolsMLJBaseExt.jl @@ -150,7 +150,7 @@ end # unsupported types of predictions and targets function _rstar(::Any, predictions, targets) - throw( + return throw( ArgumentError( "unsupported types of predictions ($(typeof(predictions))) and targets ($(typeof(targets)))", ), From f8d7ee922f4599a9a6cc0f2a46e4696e3c28d7ea Mon Sep 17 00:00:00 2001 From: Seth Axen Date: Mon, 26 Oct 2026 12:24:00 +0100 Subject: [PATCH 27/34] refactor: Drop package name from extension module names --- Project.toml | 4 ++-- ext/{MCMCDiagnosticToolsMLJBaseExt.jl => MLJBaseExt.jl} | 2 +- ext/{MCMCDiagnosticToolsStatsBaseExt.jl => StatsBaseExt.jl} | 2 +- src/rstar.jl | 2 +- test/rstar.jl | 2 +- 5 files changed, 6 insertions(+), 6 deletions(-) rename ext/{MCMCDiagnosticToolsMLJBaseExt.jl => MLJBaseExt.jl} (99%) rename ext/{MCMCDiagnosticToolsStatsBaseExt.jl => StatsBaseExt.jl} (89%) diff --git a/Project.toml b/Project.toml index 0111c367..2269464c 100644 --- a/Project.toml +++ b/Project.toml @@ -23,8 +23,8 @@ StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91" Tables = "bd369af6-aec1-5ad0-b16a-f7cc5008161c" [extensions] -MCMCDiagnosticToolsMLJBaseExt = ["Distributions", "MLJBase", "MLJModelInterface", "Tables"] -MCMCDiagnosticToolsStatsBaseExt = "StatsBase" +MLJBaseExt = ["Distributions", "MLJBase", "MLJModelInterface", "Tables"] +StatsBaseExt = "StatsBase" [compat] AbstractFFTs = "0.5, 1" diff --git a/ext/MCMCDiagnosticToolsMLJBaseExt.jl b/ext/MLJBaseExt.jl similarity index 99% rename from ext/MCMCDiagnosticToolsMLJBaseExt.jl rename to ext/MLJBaseExt.jl index bbfb5281..696a99a3 100644 --- a/ext/MCMCDiagnosticToolsMLJBaseExt.jl +++ b/ext/MLJBaseExt.jl @@ -1,4 +1,4 @@ -module MCMCDiagnosticToolsMLJBaseExt +module MLJBaseExt using Distributions: Distributions using MCMCDiagnosticTools: diff --git a/ext/MCMCDiagnosticToolsStatsBaseExt.jl b/ext/StatsBaseExt.jl similarity index 89% rename from ext/MCMCDiagnosticToolsStatsBaseExt.jl rename to ext/StatsBaseExt.jl index b2342571..c89f0fba 100644 --- a/ext/MCMCDiagnosticToolsStatsBaseExt.jl +++ b/ext/StatsBaseExt.jl @@ -1,4 +1,4 @@ -module MCMCDiagnosticToolsStatsBaseExt +module StatsBaseExt using MCMCDiagnosticTools: MCMCDiagnosticTools using Statistics: Statistics diff --git a/src/rstar.jl b/src/rstar.jl index 118d678e..3e27c2c7 100644 --- a/src/rstar.jl +++ b/src/rstar.jl @@ -153,7 +153,7 @@ function _rstar_error_hint(io, exc, argtypes, kwargs) return nothing end -# Implemented in MCMCDiagnosticToolsMLJBaseExt +# Implemented in MLJBaseExt function _rstar end function _check_model_supports_continuous_inputs end function _check_model_supports_multiclass_targets end diff --git a/test/rstar.jl b/test/rstar.jl index 3a89d38b..44685de5 100644 --- a/test/rstar.jl +++ b/test/rstar.jl @@ -200,7 +200,7 @@ end catch e e end - println(Base.get_extension(MCMCDiagnosticTools, :MCMCDiagnosticToolsMLJBaseExt) === nothing) + println(Base.get_extension(MCMCDiagnosticTools, :MLJBaseExt) === nothing) println(err isa MethodError) err === nothing || showerror(stdout, err) """ From a03bcf0d9fc3e13e0073397b797d8d202e14c749 Mon Sep 17 00:00:00 2001 From: Seth Axen Date: Mon, 26 Oct 2026 12:29:09 +0100 Subject: [PATCH 28/34] refactor: Define model trait checks only in MLJBase extension --- ext/MLJBaseExt.jl | 6 +----- src/rstar.jl | 3 --- test/rstar.jl | 9 ++++----- 3 files changed, 5 insertions(+), 13 deletions(-) diff --git a/ext/MLJBaseExt.jl b/ext/MLJBaseExt.jl index 696a99a3..a22272f4 100644 --- a/ext/MLJBaseExt.jl +++ b/ext/MLJBaseExt.jl @@ -8,11 +8,7 @@ using Random: Random using Statistics: Statistics using Tables: Tables -import MCMCDiagnosticTools: - _check_model_supports_continuous_inputs, - _check_model_supports_multiclass_predictions, - _check_model_supports_multiclass_targets, - _rstar +import MCMCDiagnosticTools: _rstar function _rstar( rng::Random.AbstractRNG, diff --git a/src/rstar.jl b/src/rstar.jl index 3e27c2c7..c2e2c26e 100644 --- a/src/rstar.jl +++ b/src/rstar.jl @@ -155,6 +155,3 @@ end # Implemented in MLJBaseExt function _rstar end -function _check_model_supports_continuous_inputs end -function _check_model_supports_multiclass_targets end -function _check_model_supports_multiclass_predictions end diff --git a/test/rstar.jl b/test/rstar.jl index 44685de5..d3f42976 100644 --- a/test/rstar.jl +++ b/test/rstar.jl @@ -149,6 +149,7 @@ end end @testset "model traits requirements" begin + MLJBaseExt = Base.get_extension(MCMCDiagnosticTools, :MLJBaseExt) samples = randn(2, 3, 4) inputs_error = ArgumentError( @@ -156,9 +157,7 @@ end ) model = UnivariateDiscretizer() @test_throws inputs_error rstar(model, samples) - @test_throws inputs_error MCMCDiagnosticTools._check_model_supports_continuous_inputs( - model - ) + @test_throws inputs_error MLJBaseExt._check_model_supports_continuous_inputs(model) targets_error = ArgumentError( "classifier does not support vectors of multi-class labels as targets" @@ -173,10 +172,10 @@ end end for model in models @test_throws targets_error rstar(model, samples) - @test_throws targets_error MCMCDiagnosticTools._check_model_supports_multiclass_targets( + @test_throws targets_error MLJBaseExt._check_model_supports_multiclass_targets( model ) - @test_throws predictions_error MCMCDiagnosticTools._check_model_supports_multiclass_predictions( + @test_throws predictions_error MLJBaseExt._check_model_supports_multiclass_predictions( model ) end From 164bdc6924b829df84b793fbe81d395a1e53938e Mon Sep 17 00:00:00 2001 From: Seth Axen Date: Mon, 26 Oct 2026 12:29:55 +0100 Subject: [PATCH 29/34] Revert "perf: Replace broadcast over views with for loop" This reverts commit 8f502f699ab2bdf3b3737861eb552b24e6882cb7. --- src/utils.jl | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/src/utils.jl b/src/utils.jl index c48cd27e..e513b5eb 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -202,9 +202,7 @@ function _tiedrank!(ranks::AbstractVector, x::AbstractVector) stop += 1 end rank = (start + stop) / 2 + offset - for i in start:stop - ranks[perm[i]] = rank - end + @views ranks[perm[start:stop]] .= rank start = stop + 1 end return ranks From 364143c6963cd4aff3a85dde39ba7e640bef6c64 Mon Sep 17 00:00:00 2001 From: Seth Axen Date: Mon, 26 Oct 2026 12:29:55 +0100 Subject: [PATCH 30/34] Revert "chore: Remove StatsBase as package extension" This reverts commit 755914982b706ef0a42e19b9356598cfb868a5bc. --- Project.toml | 1 + src/MCMCDiagnosticTools.jl | 1 + 2 files changed, 2 insertions(+) diff --git a/Project.toml b/Project.toml index 2269464c..3f53b09e 100644 --- a/Project.toml +++ b/Project.toml @@ -13,6 +13,7 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b" Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" +StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91" StatsFuns = "4c63d2b9-4356-54db-8cca-17b64c39e42c" [weakdeps] diff --git a/src/MCMCDiagnosticTools.jl b/src/MCMCDiagnosticTools.jl index ca8037bf..75079ca8 100644 --- a/src/MCMCDiagnosticTools.jl +++ b/src/MCMCDiagnosticTools.jl @@ -3,6 +3,7 @@ module MCMCDiagnosticTools using AbstractFFTs: AbstractFFTs using AliasTables: AliasTables using SpecialFunctions: SpecialFunctions +using StatsBase: StatsBase using StatsFuns: StatsFuns, sqrt2 using LinearAlgebra: LinearAlgebra From c789ffe33f78026e2d14e0471174ba3633ca9344 Mon Sep 17 00:00:00 2001 From: Seth Axen Date: Mon, 26 Oct 2026 12:30:12 +0100 Subject: [PATCH 31/34] Revert "refactor: Move expectand proxy for mad to extension" This reverts commit f11c8727e640cbbfa508b6099f984acf2cc1e15e. --- Project.toml | 2 -- ext/StatsBaseExt.jl | 12 ------------ src/ess_rhat.jl | 4 ++++ 3 files changed, 4 insertions(+), 14 deletions(-) delete mode 100644 ext/StatsBaseExt.jl diff --git a/Project.toml b/Project.toml index 3f53b09e..820b5015 100644 --- a/Project.toml +++ b/Project.toml @@ -20,12 +20,10 @@ StatsFuns = "4c63d2b9-4356-54db-8cca-17b64c39e42c" Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" MLJBase = "a7f614a8-145f-11e9-1d2a-a57a1082229d" MLJModelInterface = "e80e1ace-859a-464e-9ed9-23947d8ae3ea" -StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91" Tables = "bd369af6-aec1-5ad0-b16a-f7cc5008161c" [extensions] MLJBaseExt = ["Distributions", "MLJBase", "MLJModelInterface", "Tables"] -StatsBaseExt = "StatsBase" [compat] AbstractFFTs = "0.5, 1" diff --git a/ext/StatsBaseExt.jl b/ext/StatsBaseExt.jl deleted file mode 100644 index c89f0fba..00000000 --- a/ext/StatsBaseExt.jl +++ /dev/null @@ -1,12 +0,0 @@ -module StatsBaseExt - -using MCMCDiagnosticTools: MCMCDiagnosticTools -using Statistics: Statistics -using StatsBase: StatsBase - -function MCMCDiagnosticTools._expectand_proxy(::typeof(StatsBase.mad), x) - x_folded = MCMCDiagnosticTools._fold_around_median(x) - return MCMCDiagnosticTools._expectand_proxy(Statistics.median, x_folded) -end - -end diff --git a/src/ess_rhat.jl b/src/ess_rhat.jl index 2c26d2a9..f31eb251 100644 --- a/src/ess_rhat.jl +++ b/src/ess_rhat.jl @@ -640,6 +640,10 @@ end function _expectand_proxy(::typeof(Statistics.std), x) return (x .- Statistics.mean(x; dims=_sample_dims(x))) .^ 2 end +function _expectand_proxy(::typeof(StatsBase.mad), x) + x_folded = _fold_around_median(x) + return _expectand_proxy(Statistics.median, x_folded) +end function _expectand_proxy(f::Base.Fix2{typeof(Statistics.quantile),<:Real}, x) y = similar(x) # currently quantile does not support a dims keyword argument From a0b79dee6e4e5c2aa34235455f5c799073291d1b Mon Sep 17 00:00:00 2001 From: Seth Axen Date: Mon, 26 Oct 2026 12:30:30 +0100 Subject: [PATCH 32/34] Revert "feat: Roll our own _counts function" This reverts commit b54e546ad08f35af6cc3d93bd3a660b6eaa8d7bd. --- src/rafterydiag.jl | 4 ++-- src/utils.jl | 13 ------------- test/utils.jl | 15 --------------- 3 files changed, 2 insertions(+), 30 deletions(-) diff --git a/src/rafterydiag.jl b/src/rafterydiag.jl index 639a3000..fe888bda 100644 --- a/src/rafterydiag.jl +++ b/src/rafterydiag.jl @@ -44,7 +44,7 @@ function rafterydiag( test = dichot[1:kthin:nx] ntest = length(test) temp = test[1:(ntest - 2)] + 2 * test[2:(ntest - 1)] + 4 * test[3:ntest] - trantest = reshape(_counts(temp, 0:7), 2, 2, 2) + trantest = reshape(StatsBase.counts(temp, 0:7), 2, 2, 2) g2 = 0.0 for i1 in 1:2, i2 in 1:2, i3 in 1:2 tt = trantest[i1, i2, i3] @@ -58,7 +58,7 @@ function rafterydiag( bic = g2 - 2.0 * log(ntest - 2.0) end - tranfinal = _counts(test[1:(ntest - 1)] + 2 * test[2:ntest], 0:3) + tranfinal = StatsBase.counts(test[1:(ntest - 1)] + 2 * test[2:ntest], 0:3) alpha = tranfinal[3] / (tranfinal[1] + tranfinal[3]) beta = tranfinal[2] / (tranfinal[2] + tranfinal[4]) kthin *= step(range) diff --git a/src/utils.jl b/src/utils.jl index e513b5eb..1af0dabb 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -208,19 +208,6 @@ function _tiedrank!(ranks::AbstractVector, x::AbstractVector) return ranks end -# count the number of occurrences of each value of `levels` in `x`. -# equivalent to `StatsBase.counts(x, levels)` -function _counts(x::AbstractArray{<:Integer}, levels::UnitRange{<:Integer}) - counts = zeros(Int, length(levels)) - offset = 1 - first(levels) - for xi in x - if xi in levels - counts[xi + offset] += 1 - end - end - return counts -end - # transform the ranks to quantiles of a standard normal distribution applying the # "α-β correction" recommended in Eq 6.10.3 of # Blom. Statistical Estimates and Transformed Beta-Variables. Wiley; New York, 1958 diff --git a/test/utils.jl b/test/utils.jl index 961d67ee..93e39e4a 100644 --- a/test/utils.jl +++ b/test/utils.jl @@ -153,21 +153,6 @@ end end end -@testset "_counts" begin - x = [0, 3, 3, 7, 8, -1, 3, 0] - @test @inferred(MCMCDiagnosticTools._counts(x, 0:7)) == [2, 0, 0, 3, 0, 0, 0, 1] - @test MCMCDiagnosticTools._counts(x, 3:3) == [3] - @test MCMCDiagnosticTools._counts(reshape(x, 2, 4), 0:7) == [2, 0, 0, 3, 0, 0, 0, 1] - @test MCMCDiagnosticTools._counts(Int[], 0:3) == zeros(Int, 4) - - @testset "consistent with StatsBase.counts" begin - @testset for levels in (0:3, 0:7, -3:5) - x = rand(-5:9, 50) - @test MCMCDiagnosticTools._counts(x, levels) == StatsBase.counts(x, levels) - end - end -end - @testset "_fold_around_median" begin @testset for sz in ((1000,), (1000, 4), (1000, 4, 8), (1000, 4, 8, 2)) x = rand(sz...) From 638a9a0fa2761f9f941eb2c207b054d1bc3f377a Mon Sep 17 00:00:00 2001 From: Seth Axen Date: Mon, 26 Oct 2026 12:30:51 +0100 Subject: [PATCH 33/34] Revert "feat: Roll our own _tiedrank!" This reverts commit f0a1550ff834ad22e56f3e16ec45832c90d1defb. --- src/utils.jl | 29 ++--------------------------- test/utils.jl | 46 ---------------------------------------------- 2 files changed, 2 insertions(+), 73 deletions(-) diff --git a/src/utils.jl b/src/utils.jl index 1af0dabb..c93b9b53 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -177,37 +177,12 @@ function _rank_normalize!(values, x) fill!(values, missing) return values end - _tiedrank!(values, x) - _normal_quantiles_from_ranks!(values, values) + rank = StatsBase.tiedrank(x) + _normal_quantiles_from_ranks!(values, rank) map!(StatsFuns.norminvcdf, values, values) return values end -# rank `x` in ascending order, assigning tied values the mean of their ranks. -# equivalent to `StatsBase.tiedrank(x)` -function _tiedrank!(ranks::AbstractArray, x::AbstractArray) - _tiedrank!(vec(ranks), vec(x)) - return ranks -end -function _tiedrank!(ranks::AbstractVector, x::AbstractVector) - perm = sortperm(x) - n = lastindex(perm) - start = firstindex(perm) - offset = 1 - start - while start ≤ n - xstart = x[perm[start]] - # find the last index `stop` of the run of values tied with xstart - stop = start - while stop < n && x[perm[stop + 1]] == xstart - stop += 1 - end - rank = (start + stop) / 2 + offset - @views ranks[perm[start:stop]] .= rank - start = stop + 1 - end - return ranks -end - # transform the ranks to quantiles of a standard normal distribution applying the # "α-β correction" recommended in Eq 6.10.3 of # Blom. Statistical Estimates and Transformed Beta-Variables. Wiley; New York, 1958 diff --git a/test/utils.jl b/test/utils.jl index 93e39e4a..8d74fdac 100644 --- a/test/utils.jl +++ b/test/utils.jl @@ -107,52 +107,6 @@ end end end -@testset "_tiedrank!" begin - @testset "basic" begin - x = [3.0, 1.0, 2.0, 1.0, 5.0, 1.0] - ranks = similar(x) - @test @inferred(MCMCDiagnosticTools._tiedrank!(ranks, x)) === ranks - @test ranks == [5, 2, 4, 2, 6, 2] - @test isempty(MCMCDiagnosticTools._tiedrank!(Float64[], Float64[])) - end - - @testset "special values" begin - # as in StatsBase, -0.0 and 0.0 are tied, while NaNs are ranked last and not tied - x = [NaN, 0.0, Inf, -0.0, NaN, -Inf, 1.0] - ranks = MCMCDiagnosticTools._tiedrank!(similar(x), x) - @test ranks == [6, 2.5, 5, 2.5, 7, 1, 4] - @test ranks == StatsBase.tiedrank(x) - end - - @testset "consistent with StatsBase.tiedrank" begin - @testset for T in (Float64, Float32, Int), sz in ((20,), (20, 3)) - x = T <: Integer ? rand(T(-3):T(3), sz) : round.(randn(T, sz); digits=1) - expected = StatsBase.tiedrank(x) - @test MCMCDiagnosticTools._tiedrank!(similar(x, Float64), x) == expected - # ranks with the same float type as `x`, as in `_rank_normalize!` - @test MCMCDiagnosticTools._tiedrank!(similar(x, float(T)), x) == expected - end - end - - @testset "array types" begin - x = round.(randn(20, 3, 2); digits=1) - # parameter slices, as passed by `_rank_normalize` - xslice = view(x, :, :, 2) - @test MCMCDiagnosticTools._tiedrank!(similar(xslice), xslice) == - StatsBase.tiedrank(xslice) - - xoff = OffsetArray(x[:, :, 1], -3, 5) - ranks = MCMCDiagnosticTools._tiedrank!(similar(xoff), xoff) - @test axes(ranks) == axes(xoff) - @test parent(ranks) == StatsBase.tiedrank(parent(xoff)) - - xvoff = OffsetArray(x[:, 1, 1], -5) - ranks = MCMCDiagnosticTools._tiedrank!(similar(xvoff), xvoff) - @test axes(ranks) == axes(xvoff) - @test parent(ranks) == StatsBase.tiedrank(parent(xvoff)) - end -end - @testset "_fold_around_median" begin @testset for sz in ((1000,), (1000, 4), (1000, 4, 8), (1000, 4, 8, 2)) x = rand(sz...) From 130dc8181290ca6920db9e044e3b32a5ab91ed43 Mon Sep 17 00:00:00 2001 From: Seth Axen Date: Mon, 26 Oct 2026 12:30:59 +0100 Subject: [PATCH 34/34] test: Remove unused StatsBase import --- test/utils.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/test/utils.jl b/test/utils.jl index 8d74fdac..db853627 100644 --- a/test/utils.jl +++ b/test/utils.jl @@ -3,7 +3,6 @@ using Test using OffsetArrays using Random using Statistics -using StatsBase: StatsBase @testset "unique_indices" begin @testset "indices=$(eachindex(inds))" for inds in [