diff --git a/Project.toml b/Project.toml index 98a82a85..820b5015 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.3.20" authors = ["David Widmann", "Seth Axen"] [workspace] @@ -8,24 +8,29 @@ projects = ["docs", "test"] [deps] AbstractFFTs = "621f4979-c628-5d54-868e-fcf4e3e8185c" -DataAPI = "9a962f9c-6df0-11e9-0e5d-c546b8b5ee8a" -DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8" -Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" +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" StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91" StatsFuns = "4c63d2b9-4356-54db-8cca-17b64c39e42c" + +[weakdeps] +Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" +MLJBase = "a7f614a8-145f-11e9-1d2a-a57a1082229d" +MLJModelInterface = "e80e1ace-859a-464e-9ed9-23947d8ae3ea" Tables = "bd369af6-aec1-5ad0-b16a-f7cc5008161c" +[extensions] +MLJBaseExt = ["Distributions", "MLJBase", "MLJModelInterface", "Tables"] + [compat] AbstractFFTs = "0.5, 1" -DataAPI = "1.12" -DataStructures = "0.18.13, 0.19" +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/MLJBaseExt.jl b/ext/MLJBaseExt.jl new file mode 100644 index 00000000..a22272f4 --- /dev/null +++ b/ext/MLJBaseExt.jl @@ -0,0 +1,156 @@ +module MLJBaseExt + +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: _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) + return 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 ecd049b7..75079ca8 100644 --- a/src/MCMCDiagnosticTools.jl +++ b/src/MCMCDiagnosticTools.jl @@ -1,14 +1,10 @@ module MCMCDiagnosticTools using AbstractFFTs: AbstractFFTs -using DataAPI: DataAPI -using DataStructures: DataStructures -using Distributions: Distributions -using MLJModelInterface: MLJModelInterface as MMI +using AliasTables: AliasTables using SpecialFunctions: SpecialFunctions using StatsBase: StatsBase using StatsFuns: StatsFuns, sqrt2 -using Tables: Tables using LinearAlgebra: LinearAlgebra using Random: Random @@ -35,4 +31,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/discretediag.jl b/src/discretediag.jl index 0c0071dc..034927ae 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) @@ -206,20 +206,21 @@ are independent random variables. 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() # 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] = rand(rng, sampler) 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) + rand(rng, sampler) end end end @@ -227,12 +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) - X[1] = StatsBase.sample(1:n) + rng = Random.default_rng() + X[1] = rand(rng, 1:length(samplers)) for i in 2:N - X[i] = StatsBase.wsample(1:n, vec(P[X[i - 1], :])) + X[i] = rand(rng, samplers[X[i - 1]]) end return X end @@ -303,18 +314,18 @@ 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) - sampler_phat = Distributions.sampler(Distributions.Categorical(phat)) + sampler_phat = AliasTables.AliasTable(phat) bstats = zeros(nsim) Y = Matrix{Int}(undef, t, d) for b in 1:nsim @@ -326,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 @@ -339,13 +351,14 @@ 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 + 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/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 diff --git a/src/rafterydiag.jl b/src/rafterydiag.jl index a6960838..fe888bda 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 diff --git a/src/rstar.jl b/src/rstar.jl index f4f2e6ff..c2e2c26e 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 """ @@ -149,6 +59,10 @@ statistic is returned (algorithm 1). If the classifier is probabilistic, i.e., i probabilities of classes, the scaled Poisson-binomial distribution of the ``R^*`` statistic 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 internally in the statistic. @@ -232,43 +146,12 @@ 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") - - # 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 +function _rstar_error_hint(io, exc, argtypes, kwargs) + 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 MLJBaseExt +function _rstar end diff --git a/src/utils.jl b/src/utils.jl index 0fb3f721..c93b9b53 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 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" 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/rstar.jl b/test/rstar.jl index 763dfe31..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 @@ -189,6 +188,35 @@ end @test_throws ArgumentError MCMCDiagnosticTools._rstar(1.0, rand(2), rand(2)) end + @testset "rstar without MLJBase" begin + # MLJBase is loaded here, so check in a new process that only loads + # MCMCDiagnosticTools + code = """ + using MCMCDiagnosticTools + err = try + rstar(nothing, randn(10, 2)) + nothing + catch e + e + end + println(Base.get_extension(MCMCDiagnosticTools, :MLJBaseExt) === nothing) + println(err isa MethodError) + 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("no method matching _rstar(", msg) + @test occursin("requires MLJBase to be loaded", msg) + + # with MLJBase loaded, other `MethodError`s get no hint + @test !occursin( + "requires MLJBase to be loaded", sprint(showerror, MethodError(rstar, (1,))) + ) + end + @testset "single chain: method ambiguity issue" begin samples = rand(1:5, N) rng = MersenneTwister(42)