bat / bat.jl Goto Github PK

A Bayesian Analysis Toolkit in Julia

License: Other

Julia 100.00%

julia statistics bayesian-inference bayesian-statistics mcmc machine-learning prior posterior

bat.jl's Introduction

Welcome to BAT, a Bayesian analysis toolkit in Julia.

BAT.jl offers a variety of posterior sampling, mode estimation and integration algorithms, supplemented by plotting recipes and I/O functionality.

BAT.jl originated as a rewrite/redesign of BAT, the Bayesian Analysis Toolkit in C++. BAT.jl now offer a different set of functionality and a wider variety of algorithms than its C++ predecessor.

Installation

To install BAT.jl, start Julia and run

julia> using Pkg
julia> pkg"add BAT"

!!! note

BAT.jl requires Julia >= v1.6, we recommend to use Julia >= v1.9 for optimal performance.

Documentation

Citing BAT.jl

When using BAT.jl for research, teaching or similar, please cite our work, see CITATION.bib.

bat.jl's People

Contributors

Stargazers

Watchers

bat.jl's Issues

Add CodeMeta file

See https://escape-ossr.gitlab.io/ossr-pages/page/contribute/publish_overview/ and https://codemeta.github.io/

bat_marginalize should come with an algorithm type and defaults

bat_marginalize should always return a NamedTuple (result = ...,), to leave room for additional, algorithm-specific output values. This is something like a standard for all bat_....() functions now. Sorry, forgot to mention that on the PR that added bat_marginalize.

CC @Cornelius-G, @VasylHafych

Tutorial is missing package installation

The first line here: https://bat.github.io/BAT.jl/dev/tutorial/#Tutorial-1 is failing because some packages are not installed.

`bat_report()` error: Tuple field type cannot be `Union{}`

Is this a bug? The samples seem sane...

using StatsBase
using Statistics
using Distributions
using BAT
using DensityInterface
using IntervalSets

counts = (pass = [3, 102, 3], fail = [3, 4, 1])

likelihood = let data = counts
    
    function logpdf_poisson(ν, μ)
        return logpdf(Poisson(μ + eps(μ)), ν)
    end
    
    logfuncdensity(function (parameters)
        p = parameters
        return (
            logpdf_poisson(counts.pass[1] + counts.pass[3], 2 * p.μb * p.ϵb) +
            logpdf_poisson(counts.fail[1] + counts.fail[3], 2 * p.μb * (1 - p.ϵb)) +
            logpdf_poisson(counts.pass[2], p.μb * p.ϵb + p.μs * p.ϵs) +
            logpdf_poisson(counts.fail[2], p.μb + p.μs - p.μb * p.ϵb - p.μs * p.ϵs)
        )      
    end)
end

prior = distprod(
    μb = 0..200,
    μs = 0..200,
    ϵb = 0..(1-eps()),
    ϵs = 0..(1-eps()),
)

posterior = PosteriorMeasure(likelihood, prior)
samples = bat_sample(posterior, MCMCSampling(mcalg = MetropolisHastings(), nsteps = 10^5, nchains = 4)).result
bat_report(samples)

Output:

Tuple field type cannot be Union{}

Stacktrace:
 [1] #s1#1
   @ ~/.julia/packages/TypedTables/ItVth/src/TypedTables.jl:22 [inlined]
 [2] var"#s1#1"(names::Any, T::Any, ::Any, a::Any)
   @ TypedTables ./none:0
 [3] (::Core.GeneratedFunctionStub)(::UInt64, ::LineNumberNode, ::Any, ::Vararg{Any})
   @ Core ./boot.jl:602
 [4] _table(nt::@NamedTuple{parameter::Vector{Symbol}, value::Vector{Union{}}})
   @ TypedTables ~/.julia/packages/TypedTables/ItVth/src/Table.jl:31
 [5] Table
   @ ~/.julia/packages/TypedTables/ItVth/src/Table.jl:24 [inlined]
 [6] fixed_parameter_table(smplv::StructArrays.StructVector{DensitySample{@NamedTuple{μb::Float64, μs::Float64, ϵb::Float64, ϵs::Float64}, Float64, Int64, BAT.MCMCSampleID, Nothing}, @NamedTuple{v::ValueShapes.ShapedAsNTArray{@NamedTuple{μb::Float64, μs::Float64, ϵb::Float64, ϵs::Float64}, 1, ArraysOfArrays.ArrayOfSimilarArrays{Float64, 1, 1, 2, ElasticArrays.ElasticMatrix{Float64, Vector{Float64}}}, ValueShapes.NamedTupleShape{(:μb, :μs, :ϵb, :ϵs), NTuple{4, ValueShapes.ValueAccessor{ValueShapes.ScalarShape{Real}}}, NamedTuple}}, logd::Vector{Float64}, weight::Vector{Int64}, info::StructArrays.StructVector{BAT.MCMCSampleID, @NamedTuple{chainid::Vector{Int32}, chaincycle::Vector{Int32}, stepno::Vector{Int64}, sampletype::Vector{Int64}}, Int64}, aux::Vector{Nothing}}, Int64})
   @ BAT ~/.julia/packages/BAT/pYfQ6/src/statistics/report.jl:63
 [7] bat_report!(md::Markdown.MD, smplv::StructArrays.StructVector{DensitySample{@NamedTuple{μb::Float64, μs::Float64, ϵb::Float64, ϵs::Float64}, Float64, Int64, BAT.MCMCSampleID, Nothing}, @NamedTuple{v::ValueShapes.ShapedAsNTArray{@NamedTuple{μb::Float64, μs::Float64, ϵb::Float64, ϵs::Float64}, 1, ArraysOfArrays.ArrayOfSimilarArrays{Float64, 1, 1, 2, ElasticArrays.ElasticMatrix{Float64, Vector{Float64}}}, ValueShapes.NamedTupleShape{(:μb, :μs, :ϵb, :ϵs), NTuple{4, ValueShapes.ValueAccessor{ValueShapes.ScalarShape{Real}}}, NamedTuple}}, logd::Vector{Float64}, weight::Vector{Int64}, info::StructArrays.StructVector{BAT.MCMCSampleID, @NamedTuple{chainid::Vector{Int32}, chaincycle::Vector{Int32}, stepno::Vector{Int64}, sampletype::Vector{Int64}}, Int64}, aux::Vector{Nothing}}, Int64})
   @ BAT ~/.julia/packages/BAT/pYfQ6/src/statistics/report.jl:95
 [8] bat_report(obj::StructArrays.StructVector{DensitySample{@NamedTuple{μb::Float64, μs::Float64, ϵb::Float64, ϵs::Float64}, Float64, Int64, BAT.MCMCSampleID, Nothing}, @NamedTuple{v::ValueShapes.ShapedAsNTArray{@NamedTuple{μb::Float64, μs::Float64, ϵb::Float64, ϵs::Float64}, 1, ArraysOfArrays.ArrayOfSimilarArrays{Float64, 1, 1, 2, ElasticArrays.ElasticMatrix{Float64, Vector{Float64}}}, ValueShapes.NamedTupleShape{(:μb, :μs, :ϵb, :ϵs), NTuple{4, ValueShapes.ValueAccessor{ValueShapes.ScalarShape{Real}}}, NamedTuple}}, logd::Vector{Float64}, weight::Vector{Int64}, info::StructArrays.StructVector{BAT.MCMCSampleID, @NamedTuple{chainid::Vector{Int32}, chaincycle::Vector{Int32}, stepno::Vector{Int64}, sampletype::Vector{Int64}}, Int64}, aux::Vector{Nothing}}, Int64})
   @ BAT ~/.julia/packages/BAT/pYfQ6/src/statistics/report.jl:17
 [9] top-level scope
   @ In[39]:1

Pullback of DistributionTransform

Hi everyone,

I got a problem regarding the Zygote.pullback and BAT.DistributionTransform, which somehow does no longer work. I use Julia 1.7 with a freshly installed BAT version (should probably be the main branch). Here is a minimal example:

using BAT
using Zygote
prior = NamedTupleDist(
ξ = BAT.StandardMvNormal(length(D))
)
pos = rand(prior)
bwd_trafo = BAT.DistributionTransform(Normal, prior)
fwd_trafo = inv(bwd_trafo)
pos_t = bwd_trafo(pos)
pullback(fwd_trafo,pos_t)

Has anyone an idea what might be going on here? Here are the first lines of the error message, let me know in case you need the full thing:

ERROR: type DataType has no field mutable
Stacktrace:
[1] getproperty
@ ./Base.jl:37 [inlined]
[2] adjoint
@ ~/.julia/packages/Zygote/zowrf/src/lib/lib.jl:281 [inlined]
[3] _pullback
@ ~/.julia/packages/ZygoteRules/AIbCs/src/adjoint.jl:65 [inlined]
[4] _pullback
@ ~/.julia/packages/ValueShapes/p5OTT/src/array_shape.jl:21 [inlined]
[5] _pullback(ctx::Zygote.Context, f::Type{ArrayShape{Float64, 1}}, args::Tuple{Int64})
@ Zygote ~/.julia/packages/Zygote/zowrf/src/compiler/interface2.jl:0
[6] _pullback
@ ~/.julia/packages/ValueShapes/p5OTT/src/array_shape.jl:27 [inlined]
[7] _pullback(ctx::Zygote.Context, f::Type{ArrayShape{Float64}}, args::Tuple{Int64})
@ Zygote ~/.julia/packages/Zygote/zowrf/src/compiler/interface2.jl:0

Doesn't seem to work with Julia 1.6.2?

When running the paper example, I get

UndefVarError: n not defined

Stacktrace:
[1] _ndims(#unused#::Type{Tuple{}})
@ TypedTables ~/.julia/packages/TypedTables/zfbS2/src/TypedTables.jl:26
[2] _ndims(#unused#::NamedTuple{(), Tuple{}})
@ TypedTables ~/.julia/packages/TypedTables/zfbS2/src/TypedTables.jl:25
[3] _table(nt::NamedTuple{(), Tuple{}})
@ TypedTables ~/.julia/packages/TypedTables/zfbS2/src/Table.jl:31
[4] Table(ts::Tables.CopiedColumns{CSV.File}; kwargs::Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
@ TypedTables ~/.julia/packages/TypedTables/zfbS2/src/Table.jl:24
[5] Table
@ ~/.julia/packages/TypedTables/zfbS2/src/Table.jl:24 [inlined]
[6] |>(x::Tables.CopiedColumns{CSV.File}, f::Type{Table})
@ Base ./operators.jl:858
[7] read(source::String, sink::Type; copycols::Bool, kwargs::Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
@ CSV ~/.julia/packages/CSV/22psr/src/CSV.jl:67
[8] read(source::String, sink::Type)
@ CSV ~/.julia/packages/CSV/22psr/src/CSV.jl:64
[9] top-level scope
@ In[27]:97
[10] eval
@ ./boot.jl:360 [inlined]
[11] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)
@ Base ./loading.jl:1116

Use auto-correlation for sample splitting in AHMI

For splitting samples in two (and also for sub-sets in the covariance estimates), we should split samples in chucks of correlation lengths

"My first BAT script" example documentation.

Need at least basic example of working with core functionality. Probably better to use Jupyter notebook with julia kernel plugin since Github could render it.

Replacement of empty subsets

In https://github.com/bat/BAT.jl/blob/master/src/integration/ahmi/harmonic_mean_integration.jl#L453 integral estimates that are nan or inf (presumably when the hyper-rectangles are empty) are replaced by the estimate from the full sample. This should not be the case

how to get(build) docs

ERROR: Load

Error: makedocs encountered an error. Terminating build
BatDocsBuildProblem.txt

Order of Samples not preserved in IO via HDF5.jl

The example from the test set "hdf5" fails, if the order of the parameters in the prior is not alphabetically sorted.

MWE

I just changed a and b in example_posterior:

using Test
using BAT
using StableRNGs
using Distributions
using ValueShapes

import HDF5


function BAT.example_posterior()
    rng = StableRNGs.StableRNG(0x4cf83495c736cac2)
    prior = NamedTupleDist(
        b = [4.2, 3.3], # switched with `a`
        a = Exponential(),  
        c = Normal(1, 3),
        d = [Weibull(), Weibull()],
        e = Beta(),
        f = MvNormal([0.3, -2.9], [1.7 0.5; 0.5 2.3])
    )
    n = totalndof(varshape(prior))
    A = randn(rng, n, n)
    likelihood = varshape(prior)(MvNormal(A * A'))
    PosteriorDensity(likelihood, prior)
end

mktempdir() do tmp_datadir
    results_filename = joinpath(tmp_datadir, "results.hdf5")
    samples = bat_sample(BAT.example_posterior(), MCMCSampling(nsteps = 1000, strict = false)).result
    bat_write(results_filename, samples)
    samples2 = bat_read(results_filename).result
    @test samples == samples2
end

Add BinningAlgorithm algo type

We currently handle choice of default binning differently in bat_marginalmode and bat_marginalize. We should add BinningAlgorithm as a new super-type for auto-binning algorithms (and also offer user-specified fixed binning).

TagBot trigger issue

This issue is used to trigger TagBot; feel free to unsubscribe.

If you haven't already, you should update your TagBot.yml to include issue comment triggers.
Please see this post on Discourse for instructions and more details.

If you'd like for me to do this for you, comment TagBot fix on this issue.
I'll open a PR within a few hours, please be patient!

Type instability in bat_sample

Hi @VasylHafych, Hi @oschulz. I apologize for the long message but I want to tell you what I did and where did I get stuck

As you asked me on Friday I'm creating an issue to explain what I found regarding type instability in bat_sample. The problem is simple:

This works: @inferred(bat_sample(MvNormal([0.4, 0.6], [2.0 1.2; 1.2 3.0]), IIDSampling(nsamples=10^6)))
This doesn't work: iid_distribution = NamedTupleDist(a = mixture_model,), @inferred(bat_sample(iid_distribution, IIDSampling(nsamples=10^6)))
This doesn't work: @inferred(bat_sample(posterior, ps)), ps = PartitionedSampling(sampler = mcmc, npartitions=4, exploration_sampler=mcmc_exp, integrator = ahmi, nmax_resampling=5). For some posterior.

The output of, e.g., @code_warntype bat_sample(iid_distribution, IIDSampling(nsamples=10^6)) is:

`Variables
#self#::Core.Const(BAT.bat_sample)
target::NamedTupleDist{(:a,), Tuple{MixtureModel{Multivariate, Continuous, MvNormal, Categorical{Float64, Vector{Float64}}}}, Tuple{ValueAccessor{ArrayShape{Real, 1}}}}
algorithm::IIDSampling

Body::NamedTuple{(:result, :optargs, :kwargs), _A} where _A<:Tuple
1 ─ nothing
│ %2 = Core.NamedTuple()::Core.Const(NamedTuple())
│ %3 = Base.pairs(%2)::Core.Const(Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}())
│ %4 = BAT.:(var"#bat_sample#166")(%3, #self#, target, algorithm)::NamedTuple{(:result, :optargs, :kwargs), _A} where _A<:Tuple
└── return %4`

where Body::NamedTuple{(:result, :optargs, :kwargs), _A} where _A<:Tuple is highlighted in red both before 1 - nothing and in %4.

After calling @which bat_sample(iid_distribution, IIDSampling(nsamples=10^6))I get that the executed code is PATH/TO/DIR/.julia/dev/BAT/src/algotypes/sampling_algorithm.jl:57 which is this code:

@inline function bat_sample(target::AnySampleable, algorithm::AbstractSamplingAlgorithm; kwargs...)
rng = bat_default_withinfo(bat_sample, Val(:rng), target)
bat_sample(rng, target, algorithm; kwargs...)
end

so it's just choosing the random number generator and I guess is not the source of the instability. Is at this point where I'm stuck. I tried to @code_warntype the implementations of IIDSampling and PartitionedSampling but I didn't get anywhere probably because I'm doing something wrong.

Add support for UltraNest

I am just opening this to have a place to track and discuss, no rush.

UltraNest is a general-purpose Bayesian inference library based on nested sampling, documented at https://johannesbuchner.github.io/UltraNest/

I think @oschulz's idea was to create a Julia package in a better form than the example Julia script at https://github.com/JohannesBuchner/UltraNest/tree/master/languages/julia . I don't know how Julia packages are distributed.

example of DataFrames data

I want more examples. I want an example of processing DataFrames data.

BAT dependencies, package load time(s) and modularization

This issue is intended to keep track of modularizing (splitting up) BAT to reduce package load time(s) and increase flexibility.

The current high load time of BAT is mainly due to it's dependencies. Dependency optiions and load-time cost, preliminary analysis:

Unavoidable expensive core deps:

Distributions (includes StatsBase)
StaticArrays (indirect through many packages)

Unavoidable non-negligible deps:

ArraysOfArrays
ValueShapes
RecipesBase
Possibly some ArrayInterface packages (ArrayInterfaceCore is now very lightweight)

Hard-to-avoid non-negligible direct/indirect deps

BangBang
Transducers (non-negligible cost on top of BangBang)
TerminalLoggers
PrettyTables (may be avoidable)
DataStructures

Cheap deps:

StructArrays: Very cheap on top of StaticArrays
AbstractDifferentiation: Almost free on top of ChainRulesCore

Autodiff choices:

ForwardDiff: Quite cheap on top of StaticArrays and Distributions
FiniteDifferences: Basically free on top of StaticArrays and Distributions
FiniteDiff: Instead of FiniteDifferences, would pull in ArrayInterface
Zygote: Expensive
(Future option) Diffractor: Quite cheap on top of StaticArrays

Math deps:

LinearMaps: About 100 ms, no deps

Statistics deps:

LogarithmicNumbers: Very cheap on top of unavoidable deps
MeasureBase: Not that expensive on top of unavoidable deps (uses LogarithmicNumbers)

Optimizer deps:

NLSolversBase: Cheap itself, but depends on ArrayInterface, ChainRulesCore, SpecialFunctions and ForwardDiff
Optim: Cheap itself, but pulls in ForwardDiff and ArrayInterface via FiniteDiff,
BAT would need only Nelder-Mead and L-BFGS for core functionality though
LBFGSB: Cheap, but indroduces binary deps
Manopt: Would be chop on top of StaticArrays if it didn't pull in ColorSchemes

Sampler deps:

AbstractMCMC: Cheap on top of StaticArrays, Distributions, Transducers, TerminalLoggers
AdvancedHMC: Quite cheap on top of AbstractMCMC
PSIS: Cheap on top of unavoidable deps

Integator deps:

CUBA: Cheap, but indroduces binary deps
MonteCarloIntegration: Almost free on top of unavoidable deps
HCubature: Almost free on top of unavoidable deps
QuadGK: Cheap, but only univariate

Deps to avoid:

Folds: Significant load time cost on top of Transducers
SciMLBase: Really expensive
GalacticOptim: Not cheap even on top of SciMLBase, also load-time interplay with SciMLBase
Optimization: Over 1000 ms load time on top of it's heavy deps (Optim, GPUArrays, RecursiveArrayTools, SciMLBase)

Expensive deps to get rid of:

DoubleFloats
Polynomials
... ?

Non-negligible deps to get rid of:

KernelDensity
PrettyTables (maybe, see above)
AdaptiveRejectionSampling (cost due to ForwardDiff though)
... ?

Trim function of AHMI not working as expected

The function defined here https://github.com/bat/BAT.jl/blob/master/src/integration/ahmi/util.jl#L2 is supposed to cut out the central 68% of the integral estimates, i.e. reject outliers. If less than 3 estimates are provided it does not make a cut.

Test:

BAT._trim( [-1., 0, 1, 2, 3]) = ([1, 5], [2, 3, 4]) looks fine,
However BAT._trim([-10., 0, 1, 2, 3]) = ([1, 2, 5], [3, 4]) while expecting ([1, 5], [2, 3, 4])
and BAT._trim([0., 0., 1, 2, 3]) = ([1, 5], [2, 3, 4]) while expecting ([5], [1, 2, 3, 4])

I propose the following as a replacement:

function _trim(a)
    
    if length(a) < 3
        return collect(range(1,stop=length(a)))
    end
        
    lower = quantile(a, 0.15865)
    upper = quantile(a, 0.84135)
    indices = findall(x-> x >= lower && x <= upper, a)
    return (setdiff(range(1,stop=length(a)), indices), indices)
end```

HMC sampler returns repeated samples

Currently, the HMC sampler seems to repeats samples the proposal are rejected in the MH step:

samples = bat_sample(PosteriorDensity(v -> LogDVal(1), NamedTupleDist(x = Normal())), 10^5, AHMC()).result
all(isone, samples.weight) == true
length(samples.v) > length(unique(samples.v.x))

If we get repeated samples from AdvancedHMC, we want to increase the weight of the sample in steps of one instead of repeating it in the output. AHMI, for example, needs unique samples.

Bridge sampling with non-integer weights

When using the bridge sampler on samples generated with non-integer weights, it crashes (stack trace below).

It is in particular the nsamples=Int(sum(second_batch.weight)) here:

BAT.jl/src/integration/bridge_sampling_integration.jl

Line 134 in e4fcc53

 proposal_samples = bat_sample(proposal_density,IIDSampling(nsamples=Int(sum(second_batch.weight)))).result 

I suggest for instance to use floor(Int,sum(second_batch.weight)) instead

ArithmeticError: <PyCall.jlwrap (in a Julia function called from Python)
JULIA: InexactError: Int64(207.62753894835834)
Stacktrace:
 [1] Int64
   @ ./float.jl:723 [inlined]
 [2] bridge_sampling_integral(target_density::BAT.PosteriorDensity{Float64, Float64, BAT.ReshapedDensity{BAT.TransformedDensity{BAT.DensityWithShape{BAT.LFDensity{PyObject}, ScalarShape{Real}}, BAT.DistributionTransform{BAT.StandardUvNormal{Float64}, Distributions.Uniform{Float64}, ScalarShape{Real}, ScalarShape{Real}}, BAT.TDNoCorr, ScalarShape{Real}}, ArrayShape{Real, 1}}, BAT.DistributionDensity{Distributions.Product{Distributions.Continuous, BAT.StandardUvNormal{Float64}, FillArrays.Fill{BAT.StandardUvNormal{Float64}, 1, Tuple{Base.OneTo{Int64}}}}, BAT.HyperRectBounds{Float64}}, ArrayShape{Real, 1}, BAT.HyperRectBounds{Float32}}, target_samples::StructArrays.StructVector{BAT.DensitySample{Vector{Float64}, Float64, Float64, Nothing, Nothing}, NamedTuple{(:v, :logd, :weight, :info, :aux), Tuple{ArrayOfSimilarArrays{Float64, 1, 1, 2, SubArray{Float64, 2, Matrix{Float64}, Tuple{UnitRange{Int64}, Base.Slice{Base.OneTo{Int64}}}, false}}, Vector{Float64}, Vector{Float64}, Vector{Nothing}, Vector{Nothing}}}, Int64}, strict::Bool, ess_alg::BAT.EffSampleSizeFromAC{BAT.GeyerAutocorLen})
   @ BAT ~/.julia/packages/BAT/Ws127/src/integration/bridge_sampling_integration.jl:134
 [3] bat_integrate_impl(target::BAT.SampledDensity{BAT.PosteriorDensity{Float64, Float64, BAT.DensityWithShape{BAT.LFDensity{PyObject}, ScalarShape{Real}}, BAT.DistributionDensity{Distributions.Uniform{Float64}, BAT.HyperRectBounds{Float64}}, ScalarShape{Real}, BAT.HyperRectBounds{Float64}}, StructArrays.StructVector{BAT.DensitySample{Float64, Float64, Float64, Nothing, Nothing}, NamedTuple{(:v, :logd, :weight, :info, :aux), Tuple{SubArray{Float64, 1, Matrix{Float64}, Tuple{Int64, Base.Slice{Base.OneTo{Int64}}}, true}, Vector{Float64}, Vector{Float64}, Vector{Nothing}, Vector{Nothing}}}, Int64}, BAT.UnknownSampleGenerator}, algorithm::BAT.BridgeSampling{BAT.PriorToGaussian, BAT.EffSampleSizeFromAC{BAT.GeyerAutocorLen}})
   @ BAT ~/.julia/packages/BAT/Ws127/src/integration/bridge_sampling_integration.jl:33
 [4] bat_integrate(target::BAT.SampledDensity{BAT.PosteriorDensity{Float64, Float64, BAT.DensityWithShape{BAT.LFDensity{PyObject}, ScalarShape{Real}}, BAT.DistributionDensity{Distributions.Uniform{Float64}, BAT.HyperRectBounds{Float64}}, ScalarShape{Real}, BAT.HyperRectBounds{Float64}}, StructArrays.StructVector{BAT.DensitySample{Float64, Float64, Float64, Nothing, Nothing}, NamedTuple{(:v, :logd, :weight, :info, :aux), Tuple{SubArray{Float64, 1, Matrix{Float64}, Tuple{Int64, Base.Slice{Base.OneTo{Int64}}}, true}, Vector{Float64}, Vector{Float64}, Vector{Nothing}, Vector{Nothing}}}, Int64}, BAT.UnknownSampleGenerator}, algorithm::BAT.BridgeSampling{BAT.PriorToGaussian, BAT.EffSampleSizeFromAC{BAT.GeyerAutocorLen}})
   @ BAT ~/.julia/packages/BAT/Ws127/src/algotypes/integration_algorithm.jl:45
 [5] invokelatest(::Any, ::Any, ::Vararg{Any, N} where N; kwargs::Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
   @ Base ./essentials.jl:708
 [6] invokelatest(::Any, ::Any, ::Vararg{Any, N} where N)
   @ Base ./essentials.jl:706
 [7] _pyjlwrap_call(f::Function, args_::Ptr{PyCall.PyObject_struct}, kw_::Ptr{PyCall.PyObject_struct})
   @ PyCall ~/.julia/packages/PyCall/7a7w0/src/callback.jl:28
 [8] pyjlwrap_call(self_::Ptr{PyCall.PyObject_struct}, args_::Ptr{PyCall.PyObject_struct}, kw_::Ptr{PyCall.PyObject_struct})
   @ PyCall ~/.julia/packages/PyCall/7a7w0/src/callback.jl:44>

Handling periodic parameters

I would like to understand how you handle the transformation of periodic parameters like an angle to an unbounded domain internally. Couldn't find any mention in the docs and also couldn't find an example or something in the https://github.com/bat/BAT.jl/blob/main/src/transforms/distribution_transform.jl file but I might have missed that.

Thanks

HMC warmup not implemented correctly?

Hi, I was updating some old notebooks and I noticed that in BAT v2.0.3 the stepsize reached In a simple example of fitting the mean of a Gaussian distribution seems off in comparison with Stan. I had noticed this and fixed it back before the v2 release but it seems to have popped up again (see fig below for a comparison of HMC in BAT and Stan).

I'll try to find time to check this in more detail but just wanted to bring it to your attention. Ideally I could make a CI test for this behaviour - but can at least share my code if someone else has more time.

Can not sample with current master (5369acd)

I just used the current master branch 5369acd and when trying to sample, I get the following error:
(also for fd5f305, with 060a49b however it works fine)

MethodError: no method matching isless(::ValueShapes.NamedTupleShape{(:a,),Tuple{ValueShapes.ValueAccessor{ValueShapes.ScalarShape{Real}}}}, ::ValueShapes.NamedTupleShape{(:a,),Tuple{ValueShapes.ValueAccessor{ValueShapes.ScalarShape{Real}}}})
Closest candidates are:
  isless(!Matched::Missing, ::Any) at missing.jl:87
  isless(!Matched::Intervals.Endpoint{T,Intervals.Direction{:Right}(),B} where B where T, ::Any) at /home/user/.julia/packages/Intervals/cpG6o/src/endpoint.jl:161
  isless(!Matched::Intervals.Endpoint{T,Intervals.Direction{:Left}(),B} where B where T, ::Any) at /home/user/.julia/packages/Intervals/cpG6o/src/endpoint.jl:171
  ...
in top-level scope at base/util.jl:155
in bat_sample at julia/BAT.jl/src/algotypes/sampling_algorithm.jl:110
in #bat_sample#70 at julia/BAT.jl/src/algotypes/sampling_algorithm.jl:111 
in bat_sample at julia/BAT.jl/src/algotypes/sampling_algorithm.jl:97 
in #bat_sample#68 at julia/BAT.jl/src/algotypes/sampling_algorithm.jl:97 
in bat_sample_impl at julia/BAT.jl/src/samplers/mcmc/mcmc_sample.jl:139
in #bat_sample_impl#159 at julia/BAT.jl/src/samplers/mcmc/mcmc_sample.jl:145
in #mcmc_sample at base/none 
in #mcmc_sample#155 at julia/BAT.jl/src/samplers/mcmc/mcmc_sample.jl:46 
in Tuple at julia/BAT.jl/src/samplers/mcmc/mcmc_sample.jl:9 
in MCMCSpec at julia/BAT.jl/src/samplers/mcmc/algorithms/mh_sampler.jl:142 
in BAT.MHIterator at julia/BAT.jl/src/densities/posterior_density.jl:46
in  at base/none
in  at base/none
in #logvalof#30 at julia/BAT.jl/src/densities/abstract_density.jl:144
in get_shaped_variate at julia/BAT.jl/src/densities/abstract_density.jl:185
in <= at base/operators.jl:317
in < at base/operators.jl:268

`HamiltonianMC()` is failed on the tutorial model

Hello, everyone,

I'm trying out BAT.jl out of curiosity.
I found that HamiltonianMC() does not work on the tutorial model.
More specifically, in the tutorial, the following syntax is used to sample from the posterior of the model.

samples = bat_sample(posterior, MCMCSampling(mcalg = MetropolisHastings(), nsteps = 10^5, nchains = 4)).result

Changing mcalg = HamiltonianMC(), return an error.

ERROR: Density evaluation with logdensityof failed at [NaN, NaN, NaN, NaN, NaN], must not evaluate to NaN, density is DistMeasure(objectid = 0x81ff2597fbd0fd5d, varshape = ValueShapes.ArrayShape{Real, 1}((5,)))

It would be greatly appreciated if someone could help with this problem.

Thank you,

Eigenfunction Error

Hi,
Seems like function eig() has to be replaced with eigvals().

BAT.jl/src/samplers/mcmc/algorithms/mh_tuner.jl

Line 108 in 4e41e73

ev = eig(Σ_new)[1]

The first one is crashing tunning cycle and giving an error.

Vasyl

Update License

Had a look at the license file, I feel flattered to be mentioned in there but I guess instead the actual copyright holders should appear. In any case, my old LMU email has been retired

Plot histogram

Time to create issues! I know this isn't a glorious first issue like Ubuntu #1 but here we go. On the dev branch 465f200 with julia v0.6, I cannot plot a histogram of the samples from this minimal example

using Distributions, PDMats, StatsBase, IntervalSets
using Base.Test
using BAT, BAT.Logging
import Plots

set_log_level!(BAT, LOG_TRACE)


Σ = PDMat([1.0 1.5; 1.5 4.0])
density = MvDistDensity(MvNormal(Σ))

algorithm = MetropolisHastings{Int}()
bounds = HyperRectBounds([-5, -8], [5, 8], reflective_bounds)

chainspec = MCMCSpec(algorithm, density, bounds)

samples, sampleids, stats = @time rand(
    chainspec,
    100000,
    4,
)
Plots.plot(samples)

julia> plot(samples, (1,2), seriestype=:histogram2d, nbins=150)
ERROR: MethodError: no method matching fit(::Type{StatsBase.Histogram{Int64,N,E} where E where N}, ::Tuple{Array{Float64,1},Array{Float64,1}}, ::Array{Int64,1}, ::Tuple{StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}},StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}}; closed=:left)
Closest candidates are:
  fit(::Type{StatsBase.Histogram{T,N,E} where E where N}, ::Tuple{Vararg{AbstractArray{T,1} where T,N}}, ::StatsBase.AbstractWeights{W,T,V} where V<:AbstractArray{T,1} where T<:Real, ::Tuple{Vararg{AbstractArray{T,1} where T,N}}; closed) where {T, N, W} at /home/beaujean/.julia/v0.6/StatsBase/src/hist.jl:289
  fit(::Type{StatsBase.Histogram}, ::Any...; kwargs...) at /home/beaujean/.julia/v0.6/StatsBase/src/hist.jl:339
  fit(::StatsBase.StatisticalModel, ::Any...) at /home/beaujean/.julia/v0.6/StatsBase/src/statmodels.jl:92 got unsupported keyword argument "closed"
  ...
Stacktrace:
 [1] #fit#95(::Array{Any,1}, ::Function, ::Type{StatsBase.Histogram}, ::Tuple{Array{Float64,1},Array{Float64,1}}, ::Vararg{Any,N} where N) at /home/beaujean/.julia/v0.6/StatsBase/src/hist.jl:339
 [2] (::StatsBase.#kw##fit)(::Array{Any,1}, ::StatsBase.#fit, ::Type{StatsBase.Histogram}, ::Tuple{Array{Float64,1},Array{Float64,1}}, ::Array{Int64,1}, ::Tuple{StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}},StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}}) at ./<missing>:0
 [3] #_make_hist#232(::Bool, ::Array{Int64,1}, ::Function, ::Tuple{Array{Float64,1},Array{Float64,1}}, ::Int64) at /home/beaujean/.julia/v0.6/Plots/src/recipes.jl:551
 [4] (::Plots.#kw##_make_hist)(::Array{Any,1}, ::Plots.#_make_hist, ::Tuple{Array{Float64,1},Array{Float64,1}}, ::Int64) at ./<missing>:0
 [5] macro expansion at /home/beaujean/.julia/v0.6/Plots/src/recipes.jl:651 [inlined]
 [6] apply_recipe(::Dict{Symbol,Any}, ::Type{Val{:histogram2d}}, ::Array{Float64,1}, ::Array{Float64,1}, ::Void) at /home/beaujean/.julia/v0.6/RecipesBase/src/RecipesBase.jl:287
 [7] _process_seriesrecipe(::Plots.Plot{Plots.GRBackend}, ::Dict{Symbol,Any}) at /home/beaujean/.julia/v0.6/Plots/src/pipeline.jl:406
 [8] _plot!(::Plots.Plot{Plots.GRBackend}, ::Dict{Symbol,Any}, ::Tuple{BAT.DensitySampleVector{Float64,Float64,Int64},Tuple{Int64,Int64}}) at /home/beaujean/.julia/v0.6/Plots/src/plot.jl:231
 [9] #plot#145(::Array{Any,1}, ::Function, ::BAT.DensitySampleVector{Float64,Float64,Int64}, ::Vararg{Any,N} where N) at /home/beaujean/.julia/v0.6/Plots/src/plot.jl:56
 [10] (::RecipesBase.#kw##plot)(::Array{Any,1}, ::RecipesBase.#plot, ::BAT.DensitySampleVector{Float64,Float64,Int64}, ::Tuple{Int64,Int64}, ::Vararg{Tuple{Int64,Int64},N} where N) at ./<missing>:0

Custom CIs and/or printing all CI values via bat_report

Hello!

I am a new user of BAT, and currently use bat_report(samples) for the credible intervals to be printed on screen. However, that only prints out the 68% CI, not the others.

Also, is it possible to change the default CIs ?

If you could add some functions to customise CI quantiles and print out all CIs on screen (all CIs shown in the marginalised plot), it could be of great help.

Type of Sample Weights

Currently, the sample weights outputed by BAT are of type Array{Int64}. While standard Julia Stats functions require weights type such as (AnalyticWeights, FrequencyWeights, ProbabilityWeights, ...).
I suggest using FrequencyWeights for Metropolis-Hastings and ProbabilityWeights for accept/reject sampling.

BAT Precomiler Error: ZeroVector not defined

Hi all,

I would like to use BAT to do some analysis on some of my spectroscopy data from AWAKE.

I have been going through the tutorial and when I try to precompile BAT I get an error about ZeroVector not being defined.

After doing some searching I found that ZeroVector is supposed to live in the Distributions module, which is a part of JuliaStats, or maybe even it was supposed to be in Base but it looks like it is stripped out for some reason. I am using version 1.4.0.

Any help is appreciated.

I included my Stacktrace:

using BAT
[ Info: Precompiling BAT [c0cd4b16-88b7-57fa-983b-ab80aecada7e]
ERROR: LoadError: LoadError: LoadError: UndefVarError: ZeroVector not defined
Stacktrace:
[1] getproperty(::Module, ::Symbol) at .\Base.jl:26
[2] top-level scope at C:\Users\moody.juliapro\JuliaPro_v1.4.0-1\packages\BAT\D0e4J\src\distributions\distribution_functions.jl:6
[3] include(::Module, ::String) at .\Base.jl:377
[4] include(::String) at C:\Users\moody.juliapro\JuliaPro_v1.4.0-1\packages\BAT\D0e4J\src\BAT.jl:5
[5] top-level scope at C:\Users\moody.juliapro\JuliaPro_v1.4.0-1\packages\BAT\D0e4J\src\distributions\distributions.jl:3
[6] include(::Module, ::String) at .\Base.jl:377
[7] include(::String) at C:\Users\moody.juliapro\JuliaPro_v1.4.0-1\packages\BAT\D0e4J\src\BAT.jl:5
[8] top-level scope at C:\Users\moody.juliapro\JuliaPro_v1.4.0-1\packages\BAT\D0e4J\src\BAT.jl:53
[9] include(::Module, ::String) at .\Base.jl:377
[10] top-level scope at none:2
[11] eval at .\boot.jl:331 [inlined]
[12] eval(::Expr) at .\client.jl:449
[13] top-level scope at .\none:3
in expression starting at C:\Users\moody.juliapro\JuliaPro_v1.4.0-1\packages\BAT\D0e4J\src\distributions\distribution_functions.jl:6
in expression starting at C:\Users\moody.juliapro\JuliaPro_v1.4.0-1\packages\BAT\D0e4J\src\distributions\distributions.jl:3
in expression starting at C:\Users\moody.juliapro\JuliaPro_v1.4.0-1\packages\BAT\D0e4J\src\BAT.jl:53
ERROR: Failed to precompile BAT [c0cd4b16-88b7-57fa-983b-ab80aecada7e] to C:\Users\moody.juliapro\JuliaPro_v1.4.0-1\compiled\v1.4\BAT\ILbBM_qbTru.ji.
Stacktrace:
[1] error(::String) at .\error.jl:33
[2] compilecache(::Base.PkgId, ::String) at .\loading.jl:1272
[3] _require(::Base.PkgId) at .\loading.jl:1029
[4] require(::Base.PkgId) at .\loading.jl:927
[5] require(::Module, ::Symbol) at .\loading.jl:922

julia> VERSION
v"1.4.0"

Linear Regression with BAT

First and foremost, thank you very much for your work and the package.

I do have experience with STAN.jl and Turing.jl, and wanted to try out BAT.jl.
However, I fail to understand the syntax to define 'likelihoods' and 'priors'.

In the previous mentioned Packages, one can define the distribution of a parameter with the convenient tilde-syntax.

The following code, shows how one could implement the posterior for a simple linear regression in Turing.
Would you help me to translate it to the BAT-syntax?

Lets assume one uses the Packages:

using Turing
using Distributions

and defines the posterior to be:

@model function posterior(x, y, yerror)

    # define prior probabilities for the two unknown parameters
    intercept ~ Normal(0, 1)
    slope ~ Normal(0, 1)

    # define the function which should be able to describe the data
    mu = slope * x .+ intercept

    # the likelihood function is assumed to me a multivariate Normal distribution
    y ~ MvNormal(mu, yerror)
end

too strong version req. of dependencies

I tried to add BAT.jl to one of my projects, which was not possible due to a version conflict in the dependencies. After checking the Project.toml of BAT.jl it seems to me that the required versions are a little bit toooo strong. For every dependency only some specific version numbers are allowed, which seems a little bit over the top?! 🤔

bat_sample error

Hi,
I followed the tutorial to define a custom log-likelihood to give that to bat_sample(). My likelihood seems to work fine, however, bat_sample() doesn't work for me. After some time I get the error below:

AssertionError: length(chains) == nchains Stacktrace: [1] mcmc_init!(rng::Random123.Philox4x{UInt64, 10}, algorithm::MetropolisHastings{BAT.MvTDistProposal, RepetitionWeighting{Int64}, AdaptiveMHTuning}, density::PosteriorDensity{BAT.TransformedDensity{BAT.GenericDensity{var"#36#37"{Vector{Float64}, typeof(pdf)}}, BAT.DistributionTransform{ValueShapes.StructVariate{NamedTuple{(:offset, :resolution, :k)}}, BAT.InfiniteSpace, BAT.MixedSpace, BAT.StandardMvNormal{Float64}, NamedTupleDist{(:offset, :resolution, :k), Tuple{Product{Continuous, Uniform{Float64}, Vector{Uniform{Float64}}}, Uniform{Float64}, Uniform{Float64}}, Tuple{ValueAccessor{ArrayShape{Real, 1}}, ValueAccessor{ScalarShape{Real}}, ValueAccessor{ScalarShape{Real}}}}}, BAT.TDNoCorr}, BAT.DistributionDensity{BAT.StandardMvNormal{Float64}, BAT.HyperRectBounds{Float32}}, BAT.HyperRectBounds{Float32}, ArrayShape{Real, 1}}, nchains::Int64, init_alg::MCMCChainPoolInit, tuning_alg::AdaptiveMHTuning, nonzero_weights::Bool, callback::Function) @ BAT ~/.julia/packages/BAT/XvOy6/src/samplers/mcmc/chain_pool_init.jl:160

Any idea what could be the issue?

Reproducer:

using Random, LinearAlgebra, Statistics, Distributions, StatsBase, Plots
using BAT, IntervalSets, ValueShapes, TypedTables 
import SpecialFunctions
using Profile, ProfileView, QuadGK

##  Error function.
function my_erf(x,coeff,base)
    return coeff/2*(1 .+ SpecialFunctions.erf.(x)) .+ base
[mcmc_BAT_Julia.pdf](https://github.com/bat/BAT.jl/files/11219143/mcmc_BAT_Julia.pdf)

end
## Step function.
function my_step(x,coeff,base)
    return coeff/2*(1 .+ sign.(x)) .+ base
end
## Model (sum of two errors or steps).
function count(p::NamedTuple{(:offset, :resolution, :k)}, x)
    step1_coeff = 6+p.k
    step2_coeff = 2-p.k
    if p.resolution> 0.0000001
        step1_x = (x .- p.offset[1])/(sqrt(2)*p.resolution)
        step1 = my_erf(step1_x,step1_coeff,4)
        
        step2_x = (x .- p.offset[2])/(sqrt(2)*p.resolution)
        step2 = my_erf(step2_x,step2_coeff,0.)
    else
        step1_x = (x .- p.offset[1])
        step1 = my_step(step1_x,step1_coeff,4)
            
        step2_x = (x .- p.offset[2])
        step2 = my_step(step2_x,step2_coeff,0)
    end
    return step1+step2
end 
## Integral of model
function get_integral(p::NamedTuple{(:offset, :resolution, :k)},low, high)
    total,error = quadgk(x -> count(p,x),low,high)
    return total
end
## PDF 
function pdf(p::NamedTuple{(:offset, :resolution, :k)},x,low, high)
    total = get_integral(p,low,high)
    value = count(p,x)
    return value/total
end
## Sampler
function sampler(p::NamedTuple{(:offset, :resolution, :k)},n,low,high)
    max = count(p,high)
    arr = Array{Float64}(undef, n)
    i = 1
    while i<=n
        y = rand()*max
        x = rand()*(high-low)+low
        if y <= count(p, x)
            arr[i] = x
            i+=1
        end
    end
    return arr
end  

true_par_values = (offset = [99, 150], resolution = 5, k = 0)
arr = sampler(true_par_values,10000,0,500)
## Unbinned log-likelihood 
likelihood = let data = arr, f =pdf
    params -> begin
        function event_log_likelihood(event)
            log(f(params,event,0,500))
        end
        return LogDVal(mapreduce(event_log_likelihood, +, data))
    end
end
## Flat priors
prior = NamedTupleDist(
    offset = [Uniform(50, 150), Uniform(80, 220)],
    resolution = Uniform(0,20),
    k = Uniform(-6,2)
)
## Posterior 
posterior = PosteriorDensity(likelihood, prior)
## running bat_sample
samples = bat_sample(posterior, MCMCSampling(mcalg = MetropolisHastings(), nsteps = 10^3)).result

Use histogram as both prior and input data?

as of:

we can use histogram as prior, in HEP (as @oschulz knows), the input data is also histogram (with same binning of course), and each uncertainty variation as one additional histogram.

Can we run likelihood for parameter extraction in such case?

Error of iid sampler with a truncated mixture distribution

Hi @oschulz,

I found one bug while testing BAT iid sampler. I am not sure if this is a problem of BAT, but if so it might be useful to fix it to be able to test AHMI properly. Here is a code example:

using Distributions, BAT, ValueShapes

mixture = MixtureModel(Normal[
   Normal(-1.0, 0.5),
   Normal(1.0, 0.5),],)

dist = product_distribution([truncated(mixture, -2, 2), truncated(Normal(0.0,1.0), -2, 2), ])

iid_distributions = NamedTupleDist(a = dist,)

samples_iid = bat_sample(iid_distributions, 10^5).result;

the error message is

ERROR: MethodError: no method matching iterate(::MixtureModel{Univariate,Continuous,Normal,Float64})
Closest candidates are:
  iterate(::Core.SimpleVector) at essentials.jl:603
  iterate(::Core.SimpleVector, ::Any) at essentials.jl:603
  iterate(::ExponentialBackOff) at error.jl:253
  ...
Stacktrace:
 [1] copyto!(::Array{Float64,1}, ::MixtureModel{Univariate,Continuous,Normal,Float64}) at ./abstractarray.jl:721
 [2] _collect(::UnitRange{Int64}, ::MixtureModel{Univariate,Continuous,Normal,Float64}, ::Base.HasEltype, ::Base.HasLength) at ./array.jl:609
 [3] collect(::MixtureModel{Univariate,Continuous,Normal,Float64}) at ./array.jl:603

By the way, MCMC sampling works correctly with this distribution. Only iid sampling of truncated distribution fails .

Multil Template Fitter planned?

Hey, I wanted to ask if you planned on implementing the BCMultiTemplateFitter from the C++ version. Or if there is anything similar in the Julia version. This would be very useful when fitting multiple spectral templates from simulation to experimental reference data.

MLE best-fit parameter disagreeing with `pyhf`

Related to:

#318

using PyCall, BAT, Distributions

py"""
import pyhf
spec = {
	'channels': [
		{
			'name': 'signal',
			'samples': [
				{
					'name': 'signal',
					'data': [2,3,4,5],
					'modifiers': [
						{'name': 'mu', 'type': 'normfactor', 'data': None}
					],
				},
				{
					'name': 'bkg1',
					'data': [30,19,9,4],
					'modifiers': [
						{ 
							"name": "theta", 
							"type": "histosys", 
							"data": {
								"hi_data": [31,21,12,7], 
								"lo_data": [29,17,6,1]
							}
						}
					],
				},
			],
		}
	]
}
pdf_pyhf = pyhf.Model(spec)
"""

pdf_pyhf = py"pdf_pyhf"
v_data = [34,22,13,11] # observed data
data_pyhf = [v_data; pdf_pyhf.config.auxdata]

prior = BAT.NamedTupleDist(
    μ = Uniform(0, 4),
    θ = Normal()
)

function likelihood_pyhf(v)
	(;μ, θ) = v
	LogDVal(only(pdf_pyhf.logpdf([μ, θ], data_pyhf)))
end

@assert likelihood((;μ=1.3, θ=-0.07)).logval ≈ likelihood_pyhf((;μ=1.3, θ=-0.07)).logval

Result

pyhf = pyimport("pyhf")
(μ̂_pyhf, θ̂_pyhf), nll_pyhf = pyhf.infer.mle.fit(data_pyhf, pdf_pyhf, return_fitted_val=true)
# μ = 1.30654, θ = -0.0603666

posterior_BATpyhf = PosteriorDensity(likelihood_pyhf, prior)
best_fit_BATpyhf = bat_findmode(posterior_BATpyhf).result
# μ = 1.2872543725923415, θ = -0.030559943549792377

Sampling with a single chain

Sampling with nchains=1 results in an AssertionError:

JULIA: AssertionError: length(chains) == nchains
Stacktrace:
  [1] mcmc_init!(rng::Random123.Philox4x{UInt64, 10}, algorithm::BAT.MetropolisHastings{BAT.MvTDistProposal, BAT.RepetitionWeighting{Int64}, BAT.AdaptiveMHTuning}, density::BAT.PosteriorDensity{Float64, Float64, BAT.ReshapedDensity{BAT.TransformedDensity{BAT.DensityWithShape{BAT.LFDensity{PyObject}, ScalarShape{Real}}, BAT.DistributionTransform{BAT.StandardUvNormal{Float64}, Distributions.Uniform{Float64}, ScalarShape{Real}, ScalarShape{Real}}, BAT.TDNoCorr, ScalarShape{Real}}, ArrayShape{Real, 1}}, BAT.DistributionDensity{Distributions.Product{Distributions.Continuous, BAT.StandardUvNormal{Float64}, FillArrays.Fill{BAT.StandardUvNormal{Float64}, 1, Tuple{Base.OneTo{Int64}}}}, BAT.HyperRectBounds{Float64}}, ArrayShape{Real, 1}, BAT.HyperRectBounds{Float32}}, nchains::Int64, init_alg::BAT.MCMCChainPoolInit, tuning_alg::BAT.AdaptiveMHTuning, nonzero_weights::Bool, callback::Function)
    @ BAT ~/.julia/packages/BAT/Ws127/src/samplers/mcmc/chain_pool_init.jl:161
  [2] bat_sample_impl(rng::Random123.Philox4x{UInt64, 10}, target::BAT.PosteriorDensity{Float64, Float64, BAT.DensityWithShape{BAT.LFDensity{PyObject}, ScalarShape{Real}}, BAT.DistributionDensity{Distributions.Uniform{Float64}, BAT.HyperRectBounds{Float64}}, ScalarShape{Real}, BAT.HyperRectBounds{Float64}}, algorithm::BAT.MCMCSampling{BAT.MetropolisHastings{BAT.MvTDistProposal, BAT.RepetitionWeighting{Int64}, BAT.AdaptiveMHTuning}, BAT.PriorToGaussian, BAT.MCMCChainPoolInit, BAT.MCMCMultiCycleBurnin, BAT.BrooksGelmanConvergence, typeof(BAT.nop_func)})
    @ BAT ~/.julia/packages/BAT/Ws127/src/samplers/mcmc/mcmc_sample.jl:52
  [3] #bat_sample#159
    @ ~/.julia/packages/BAT/Ws127/src/algotypes/sampling_algorithm.jl:45 [inlined]
  [4] bat_sample
    @ ~/.julia/packages/BAT/Ws127/src/algotypes/sampling_algorithm.jl:45 [inlined]
  [5] #bat_sample#161
    @ ~/.julia/packages/BAT/Ws127/src/algotypes/sampling_algorithm.jl:59 [inlined]
  [6] bat_sample(target::BAT.PosteriorDensity{Float64, Float64, BAT.DensityWithShape{BAT.LFDensity{PyObject}, ScalarShape{Real}}, BAT.DistributionDensity{Distributions.Uniform{Float64}, BAT.HyperRectBounds{Float64}}, ScalarShape{Real}, BAT.HyperRectBounds{Float64}}, algorithm::BAT.MCMCSampling{BAT.MetropolisHastings{BAT.MvTDistProposal, BAT.RepetitionWeighting{Int64}, BAT.AdaptiveMHTuning}, BAT.PriorToGaussian, BAT.MCMCChainPoolInit, BAT.MCMCMultiCycleBurnin, BAT.BrooksGelmanConvergence, typeof(BAT.nop_func)})
    @ BAT ~/.julia/packages/BAT/Ws127/src/algotypes/sampling_algorithm.jl:58
  [7] invokelatest(::Any, ::Any, ::Vararg{Any, N} where N; kwargs::Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ Base ./essentials.jl:708
  [8] invokelatest(::Any, ::Any, ::Vararg{Any, N} where N)
    @ Base ./essentials.jl:706
  [9] _pyjlwrap_call(f::Function, args_::Ptr{PyCall.PyObject_struct}, kw_::Ptr{PyCall.PyObject_struct})
    @ PyCall ~/.julia/packages/PyCall/7a7w0/src/callback.jl:28
 [10] pyjlwrap_call(self_::Ptr{PyCall.PyObject_struct}, args_::Ptr{PyCall.PyObject_struct}, kw_::Ptr{PyCall.PyObject_struct})
    @ PyCall ~/.julia/packages/PyCall/7a7w0/src/callback.jl:44>

How does BAT differ from Turing?

Hi! Found this on Github and I'd like to say, BAT.jl looks really interesting, and I'm always happy to see new tools for Bayesian analysis in Julia! I have a quick question -- how does BAT.jl differ from other PPLs in Julia, like Turing, Soss, or Gen?

Plotting Warnings

Hi @oschulz @Cornelius-G. It looks like plotting recipes are now throwing a huge number of warning messages. For example, if I run a tutorial from the documentation I am getting:

plot(samples, 3)

┌ Warning: Attribute alias `alpha` detected in the user recipe defined for the signature (::StructArrays.StructArray{DensitySample{NamedTuple{(:a, :mu, :sigma),Tuple{Array{Float64,1},Array{Float64,1},Float64}},Float64,Int64,BAT.MCMCSampleID,Nothing},1,NamedTuple{(:v, :logd, :weight, :info, :aux),Tuple{ShapedAsNTArray{NamedTuple{(:a, :mu, :sigma),Tuple{Array{Float64,1},Array{Float64,1},Float64}},1,ArraysOfArrays.ArrayOfSimilarArrays{Float64,1,1,2,ElasticArrays.ElasticArray{Float64,2,1,Array{Float64,1}}},NamedTupleShape{(:a, :mu, :sigma),Tuple{ValueAccessor{ArrayShape{Real,1}},ValueAccessor{ArrayShape{Real,1}},ValueAccessor{ScalarShape{Real}}}}},Array{Float64,1},Array{Int64,1},StructArrays.StructArray{BAT.MCMCSampleID,1,NamedTuple{(:chainid, :chaincycle, :stepno, :sampletype),Tuple{Array{Int32,1},Array{Int32,1},Array{Int64,1},Array{Int64,1}}},Int64},Array{Nothing,1}}},CartesianIndex{1}}, ::Int64). To ensure expected behavior it is recommended to use the default attribute `seriesalpha`.
└ @ Plots /Users/vhafych/.julia/packages/Plots/bx41w/src/pipeline.jl:15
┌ Warning: Attribute alias `color` detected in the user recipe defined for the signature (::Histogram{Int64,2,Tuple{StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}},StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}}}, ::Tuple{Int64,Int64}). To ensure expected behavior it is recommended to use the default attribute `seriescolor`.
└ @ Plots /Users/vhafych/.julia/packages/Plots/bx41w/src/pipeline.jl:15
....

Can you check if you have the same problem? I am using Plots v1.1.3, warnings are present for both gr() and pyplot() backends.

Mapping between names and indices fails for constant prior parameters

I just noticed that when setting parameters of the prior constant, the mapping between parameter names and indices for the samples fails (which is needed in the plotting recipes). This is because in the shape of the samples the constant parameter is there, while there are no entries for it in the vector of unshaped.(samples).
I guess I could fix this by adding some checks to see whether a shape contains a ConstValueShape in the function asindex(vs::NamedTupleShape, name::Symbol)in valueshapes_utils.jl, but I was wondering whether it is a good idea to drop the fixed parameter from the unshaped samples instead of explicitly filling the vectors with the fixed value. I guess it might be needed when using the samples for further computations.

Support pyhf

Hello,
I believe the previous version of BAT had an interface for working with RooFit likelihoods (that can be stored in a RooFit workspace in a ROOT file). I am curious if there are plans to provide such a module, perhaps in a separate repository. I know there is a generic API for likelihoods.

If I am correct that this was supported in the previous version, it might be worth mentioning the current status of this type of interface in BAT.jl.

Thanks!

MCMC diagnostic plots

Hi,

I am reaching out to inquire about the possibility of incorporating methods for plotting diagnostic figures related to MCMC tuning in this project (e.g. trace and acf plots).

While exploring the project, I came across an example from four years ago that seems to demonstrate a MCMCDiagnostics method (Example Link). However, I could not find this method in the documentation, and my attempts to utilize it based on the example were unsuccessful.

Could it be possible to have documentation or an updated example on the MCMCDiagnostics method? Alternatively, if this feature is not present, would it be possible to consider adding methods for generating trace and acf plots during MCMC tuning?

Thanks in advance.

AHMI fails with multiple workers

Hi @oschulz! As I promised, here is an example of code that shows that AHMI fails when multiple workers are used.

using Distributed

Distributed.addprocs(3)

@everywhere begin

    using ValueShapes
    using ArraysOfArrays
    using Distributions
    using BAT

    samples = bat_sample(NamedTupleDist(a = Normal(1,1)), 10^5).result
    hmi_data_iid = BAT.HMIData(unshaped.(samples))
    BAT.hm_integrate!(hmi_data_iid)
    
end

I am getting the following error

ERROR: LoadError: On worker 3:
ArgumentError: List of selected workers must be sorted and contain no duplicates
workpart at /Users/vhafych/.julia/packages/ParallelProcessingTools/7dVGh/src/workpartition.jl:96
hm_update_integrationvolumes_dataset! at /Users/vhafych/MPP-Server/gitrepos/BAT.jl/src/integration/ahmi/hm_integration_rectangle.jl:81
hm_create_integrationvolumes! at /Users/vhafych/MPP-Server/gitrepos/BAT.jl/src/integration/ahmi/hm_integration_rectangle.jl:30
#hm_integrate!#131 at /Users/vhafych/MPP-Server/gitrepos/BAT.jl/src/integration/ahmi/harmonic_mean_integration.jl:55
hm_integrate! at /Users/vhafych/MPP-Server/gitrepos/BAT.jl/src/integration/ahmi/harmonic_mean_integration.jl:37 [inlined] (repeats 2 times)
top-level scope at /Users/vhafych/MPP-Server/AHMI_workers_issue/test.jl:14
eval at ./boot.jl:331
#101 at /Users/sabae/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.4/Distributed/src/process_messages.jl:290
run_work_thunk at /Users/sabae/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.4/Distributed/src/process_messages.jl:79
run_work_thunk at /Users/sabae/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.4/Distributed/src/process_messages.jl:88
#94 at ./task.jl:358

...and 2 more exception(s).

Stacktrace:
 [1] sync_end(::Array{Any,1}) at ./task.jl:316
 [2] macro expansion at ./task.jl:335 [inlined]
 [3] remotecall_eval(::Module, ::Array{Int64,1}, ::Expr) at /Users/sabae/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.4/Distributed/src/macros.jl:217
....

Can you please have a look at it?