API
Summary statistics
PosteriorStats.SummaryStats — Typestruct SummaryStats{D, V<:(AbstractVector)}A container for a column table of values computed by summarize.
This object implements the Tables and TableTraits column table interfaces. It has a custom show method.
SummaryStats behaves like an OrderedDict of columns, where the columns can be accessed using either Symbols or a 1-based integer index.
name::String: The name of the collection of summary statistics, used as the table title in display.data::Any: The summary statistics for each parameter. It must implement the Tables interface.parameter_names::AbstractVector: Names of the parameters
SummaryStats([name::String,] data[, parameter_names])
-SummaryStats(data[, parameter_names]; name::String="SummaryStats")Construct a SummaryStats from tabular data with optional stats name and param_names.
data must not contain a column :parameter, as this is reserved for the parameter names, which are always in the first column.
PosteriorStats.default_diagnostics — Functiondefault_diagnostics(focus=Statistics.mean; kwargs...)Default diagnostics to be computed with summarize.
The value of focus determines the diagnostics to be returned:
Statistics.mean:mcse_mean,mcse_std,ess_tail,ess_bulk,rhatStatistics.median:mcse_median,ess_tail,ess_bulk,rhat
PosteriorStats.default_stats — Functiondefault_stats(focus=Statistics.mean; prob_interval=0.94, kwargs...)Default statistics to be computed with summarize.
The value of focus determines the statistics to be returned:
Statistics.mean:mean,std,hdi_94%Statistics.median:median,mad,eti_94%
If prob_interval is set to a different value than the default, then different HDI and ETI statistics are computed accordingly. hdi refers to the highest-density interval, while eti refers to the equal-tailed interval.
PosteriorStats.default_summary_stats — Functiondefault_summary_stats(focus=Statistics.mean; kwargs...)Combinatiton of default_stats and default_diagnostics to be used with summarize.
PosteriorStats.summarize — Functionsummarize(data, stats_funs...; name="SummaryStats", [var_names]) -> SummaryStatsCompute the summary statistics in stats_funs on each param in data.
stats_funs is a collection of functions that reduces a matrix with shape (draws, chains) to a scalar or a collection of scalars. Alternatively, an item in stats_funs may be a Pair of the form name => fun specifying the name to be used for the statistic or of the form (name1, ...) => fun when the function returns a collection. When the function returns a collection, the names in this latter format must be provided.
If no stats functions are provided, then those specified in default_summary_stats are computed.
var_names specifies the names of the parameters in data. If not provided, the names are inferred from data.
To support computing summary statistics from a custom object, overload this method specifying the type of data.
See also SummaryStats, default_summary_stats, default_stats, default_diagnostics.
Examples
Compute mean, std and the Monte Carlo standard error (MCSE) of the mean estimate:
julia> using Statistics, StatsBase
+SummaryStats(data[, parameter_names]; name::String="SummaryStats")Construct a SummaryStats from tabular data with optional stats name and param_names.
data must not contain a column :parameter, as this is reserved for the parameter names, which are always in the first column.
PosteriorStats.default_diagnostics — Functiondefault_diagnostics(focus=Statistics.mean; kwargs...)Default diagnostics to be computed with summarize.
The value of focus determines the diagnostics to be returned:
Statistics.mean:mcse_mean,mcse_std,ess_tail,ess_bulk,rhatStatistics.median:mcse_median,ess_tail,ess_bulk,rhat
PosteriorStats.default_stats — Functiondefault_stats(focus=Statistics.mean; prob_interval=0.94, kwargs...)Default statistics to be computed with summarize.
The value of focus determines the statistics to be returned:
Statistics.mean:mean,std,hdi_94%Statistics.median:median,mad,eti_94%
If prob_interval is set to a different value than the default, then different HDI and ETI statistics are computed accordingly. hdi refers to the highest-density interval, while eti refers to the equal-tailed interval.
PosteriorStats.default_summary_stats — Functiondefault_summary_stats(focus=Statistics.mean; kwargs...)Combinatiton of default_stats and default_diagnostics to be used with summarize.
PosteriorStats.summarize — Functionsummarize(data, stats_funs...; name="SummaryStats", [var_names]) -> SummaryStatsCompute the summary statistics in stats_funs on each param in data.
stats_funs is a collection of functions that reduces a matrix with shape (draws, chains) to a scalar or a collection of scalars. Alternatively, an item in stats_funs may be a Pair of the form name => fun specifying the name to be used for the statistic or of the form (name1, ...) => fun when the function returns a collection. When the function returns a collection, the names in this latter format must be provided.
If no stats functions are provided, then those specified in default_summary_stats are computed.
var_names specifies the names of the parameters in data. If not provided, the names are inferred from data.
To support computing summary statistics from a custom object, overload this method specifying the type of data.
See also SummaryStats, default_summary_stats, default_stats, default_diagnostics.
Examples
Compute mean, std and the Monte Carlo standard error (MCSE) of the mean estimate:
julia> using Statistics, StatsBase
julia> x = randn(1000, 4, 3) .+ reshape(0:10:20, 1, 1, :);
@@ -37,7 +37,7 @@
q5 q25 q50 q75 q95
1 -1.61 -0.668 0.00447 0.653 1.64
2 8.41 9.34 10.0 10.7 11.6
- 3 18.4 19.3 20.0 20.6 21.6Credible intervals
PosteriorStats.hdi — Functionhdi(samples::AbstractVecOrMat{<:Real}; [prob, sorted, method]) -> IntervalSets.ClosedInterval
+ 3 18.4 19.3 20.0 20.6 21.6Credible intervals
PosteriorStats.hdi — Functionhdi(samples::AbstractVecOrMat{<:Real}; [prob, sorted, method]) -> IntervalSets.ClosedInterval
hdi(samples::AbstractArray{<:Real}; [prob, sorted, method]) -> Array{<:IntervalSets.ClosedInterval}Estimate the highest density interval (HDI) of samples for the probability prob.
The HDI is the minimum width Bayesian credible interval (BCI). That is, it is the smallest possible interval containing (100*prob)% of the probability mass.[Hyndman1996] This implementation uses the algorithm of [ChenShao1999].
Arguments
samples: an array of shape(draws[, chains[, params...]]). If multiple parameters are present, a marginal HDI is computed for each.
Keywords
prob: the probability mass to be contained in the HDI. Default is0.94.sorted=false: iftrue, the input samples are assumed to be sorted.method::Symbol: the method used to estimate the HDI. Available options are::unimodal: Assumes a unimodal distribution (default). Bounds are entries insamples.:multimodal: Fits a density estimator tosamplesand finds the HDI of the estimated density. For continuous data, the density estimator is a kernel density estimate (KDE) computed usingkde_reflected. For discrete data, a histogram with bin width 1 is used.:multimodal_sample: Like:multimodal, but uses the density estimator to find the entries insampleswith the highest density and computes the HDI from those points. This can correct for inaccuracies in the density estimator.
is_discrete::Union{Bool,Nothing}=nothing: Specify if the data is discrete (integer-valued). Ifnothing, it's automatically determined.kwargs: For continuous data and multimodalmethods, remaining keywords are forwarded tokde_reflected.
Returns
intervals: Ifsamplesis a vector or matrix, then a singleIntervalSets.ClosedIntervalis returned for:unimodalmethod, or a vector ofIntervalSets.ClosedIntervalfor multimodal methods. For higher dimensional inputs, an array with the shape(params...,)is returned, containing marginal HDIs for each parameter.
Any default value of prob is arbitrary. The default value of prob=0.94 instead of a more common default like prob=0.95 is chosen to remind the user of this arbitrariness.
Examples
Here we calculate the 83% HDI for a normal random variable:
julia> x = randn(2_000);
julia> hdi(x; prob=0.83)
@@ -52,7 +52,7 @@
julia> hdi(x; method=:multimodal)
2-element Vector{IntervalSets.ClosedInterval{Float64}}:
-1.8882401079608677 .. 2.0017686164555037
- 2.9839268475847436 .. 6.9235952578447275PosteriorStats.hdi! — Functionhdi!(samples::AbstractArray{<:Real}; [prob, sorted])A version of hdi that partially sorts samples in-place while computing the HDI.
PosteriorStats.eti — Functioneti(samples::AbstractVecOrMat{<:Real}; [prob, kwargs...]) -> IntervalSets.ClosedInterval
+ 2.9839268475847436 .. 6.9235952578447275PosteriorStats.hdi! — Functionhdi!(samples::AbstractArray{<:Real}; [prob, sorted])A version of hdi that partially sorts samples in-place while computing the HDI.
PosteriorStats.eti — Functioneti(samples::AbstractVecOrMat{<:Real}; [prob, kwargs...]) -> IntervalSets.ClosedInterval
eti(samples::AbstractArray{<:Real}; [prob, kwargs...]) -> Array{<:IntervalSets.ClosedInterval}Estimate the equal-tailed interval (ETI) of samples for the probability prob.
The ETI of a given probability is the credible interval wih the property that the probability of being below the interval is equal to the probability of being above it. That is, it is defined by the (1-prob)/2 and 1 - (1-prob)/2 quantiles of the samples.
Arguments
samples: an array of shape(draws[, chains[, params...]]). If multiple parameters are present
Keywords
prob: the probability mass to be contained in the ETI. Default is0.94.kwargs: remaining keywords are passed toStatistics.quantile.
Returns
intervals: Ifsamplesis a vector or matrix, then a singleIntervalSets.ClosedIntervalis returned. Otherwise, an array with the shape(params...,), is returned, containing a marginal ETI for each parameter.
Any default value of prob is arbitrary. The default value of prob=0.94 instead of a more common default like prob=0.95 is chosen to reminder the user of this arbitrariness.
Examples
Here we calculate the 83% ETI for a normal random variable:
julia> x = randn(2_000);
julia> eti(x; prob=0.83)
@@ -62,7 +62,7 @@
3-element Vector{IntervalSets.ClosedInterval{Float64}}:
-1.951006825019686 .. 1.9011666217153793
3.048993174980314 .. 6.9011666217153795
- 8.048993174980314 .. 11.90116662171538PosteriorStats.eti! — Functioneti!(samples::AbstractArray{<:Real}; [prob, kwargs...])A version of eti that partially sorts samples in-place while computing the ETI.
LOO and WAIC
PosteriorStats.AbstractELPDResult — Typeabstract type AbstractELPDResultAn abstract type representing the result of an ELPD computation.
Every subtype stores estimates of both the expected log predictive density (elpd) and the effective number of parameters p, as well as standard errors and pointwise estimates of each, from which other relevant estimates can be computed.
Subtypes implement the following functions:
PosteriorStats.PSISLOOResult — TypeResults of Pareto-smoothed importance sampling leave-one-out cross-validation (PSIS-LOO).
See also: loo, AbstractELPDResult
estimates: Estimates of the expected log pointwise predictive density (ELPD) and effective number of parameters (p)pointwise: Pointwise estimatespsis_result: Pareto-smoothed importance sampling (PSIS) results
PosteriorStats.WAICResult — TypeResults of computing the widely applicable information criterion (WAIC).
See also: waic, AbstractELPDResult
estimates: Estimates of the expected log pointwise predictive density (ELPD) and effective number of parameters (p)pointwise: Pointwise estimates
PosteriorStats.elpd_estimates — Functionelpd_estimates(result::AbstractELPDResult; pointwise=false) -> (; elpd, elpd_mcse, lpd)Return the (E)LPD estimates from the result.
PosteriorStats.information_criterion — Functioninformation_criterion(elpd, scale::Symbol)Compute the information criterion for the given scale from the elpd estimate.
scale must be one of (:deviance, :log, :negative_log).
information_criterion(result::AbstractELPDResult, scale::Symbol; pointwise=false)Compute information criterion for the given scale from the existing ELPD result.
scale must be one of (:deviance, :log, :negative_log).
If pointwise=true, then pointwise estimates are returned.
PosteriorStats.loo — Functionloo(log_likelihood; reff=nothing, kwargs...) -> PSISLOOResult{<:NamedTuple,<:NamedTuple}Compute the Pareto-smoothed importance sampling leave-one-out cross-validation (PSIS-LOO). [Vehtari2017][LOOFAQ]
log_likelihood must be an array of log-likelihood values with shape (chains, draws[, params...]).
Keywords
reff::Union{Real,AbstractArray{<:Real}}: The relative effective sample size(s) of the likelihood values. If an array, it must have the same data dimensions as the corresponding log-likelihood variable. If not provided, then this is estimated usingMCMCDiagnosticTools.ess.kwargs: Remaining keywords are forwarded to [PSIS.psis].
See also: PSISLOOResult, waic
Examples
Manually compute $R_\mathrm{eff}$ and calculate PSIS-LOO of a model:
julia> using ArviZExampleData, MCMCDiagnosticTools
+ 8.048993174980314 .. 11.90116662171538PosteriorStats.eti! — Functioneti!(samples::AbstractArray{<:Real}; [prob, kwargs...])A version of eti that partially sorts samples in-place while computing the ETI.
LOO and WAIC
PosteriorStats.AbstractELPDResult — Typeabstract type AbstractELPDResultAn abstract type representing the result of an ELPD computation.
Every subtype stores estimates of both the expected log predictive density (elpd) and the effective number of parameters p, as well as standard errors and pointwise estimates of each, from which other relevant estimates can be computed.
Subtypes implement the following functions:
PosteriorStats.PSISLOOResult — TypeResults of Pareto-smoothed importance sampling leave-one-out cross-validation (PSIS-LOO).
See also: loo, AbstractELPDResult
estimates: Estimates of the expected log pointwise predictive density (ELPD) and effective number of parameters (p)pointwise: Pointwise estimatespsis_result: Pareto-smoothed importance sampling (PSIS) results
PosteriorStats.WAICResult — TypeResults of computing the widely applicable information criterion (WAIC).
See also: waic, AbstractELPDResult
estimates: Estimates of the expected log pointwise predictive density (ELPD) and effective number of parameters (p)pointwise: Pointwise estimates
PosteriorStats.elpd_estimates — Functionelpd_estimates(result::AbstractELPDResult; pointwise=false) -> (; elpd, elpd_mcse, lpd)Return the (E)LPD estimates from the result.
PosteriorStats.information_criterion — Functioninformation_criterion(elpd, scale::Symbol)Compute the information criterion for the given scale from the elpd estimate.
scale must be one of (:deviance, :log, :negative_log).
information_criterion(result::AbstractELPDResult, scale::Symbol; pointwise=false)Compute information criterion for the given scale from the existing ELPD result.
scale must be one of (:deviance, :log, :negative_log).
If pointwise=true, then pointwise estimates are returned.
PosteriorStats.loo — Functionloo(log_likelihood; reff=nothing, kwargs...) -> PSISLOOResult{<:NamedTuple,<:NamedTuple}Compute the Pareto-smoothed importance sampling leave-one-out cross-validation (PSIS-LOO). [Vehtari2017][LOOFAQ]
log_likelihood must be an array of log-likelihood values with shape (chains, draws[, params...]).
Keywords
reff::Union{Real,AbstractArray{<:Real}}: The relative effective sample size(s) of the likelihood values. If an array, it must have the same data dimensions as the corresponding log-likelihood variable. If not provided, then this is estimated usingMCMCDiagnosticTools.ess.kwargs: Remaining keywords are forwarded to [PSIS.psis].
See also: PSISLOOResult, waic
Examples
Manually compute $R_\mathrm{eff}$ and calculate PSIS-LOO of a model:
julia> using ArviZExampleData, MCMCDiagnosticTools
julia> idata = load_example_data("centered_eight");
@@ -79,7 +79,7 @@
Pareto shape (k) diagnostic values:
Count Min. ESS
(-Inf, 0.5] good 7 (87.5%) 151
- (0.5, 0.7] okay 1 (12.5%) 446PosteriorStats.waic — Functionwaic(log_likelihood::AbstractArray) -> WAICResult{<:NamedTuple,<:NamedTuple}Compute the widely applicable information criterion (WAIC).[Watanabe2010][Vehtari2017][LOOFAQ]
log_likelihood must be an array of log-likelihood values with shape (chains, draws[, params...]).
See also: WAICResult, loo
Examples
Calculate WAIC of a model:
julia> using ArviZExampleData
+ (0.5, 0.7] okay 1 (12.5%) 446PosteriorStats.waic — Functionwaic(log_likelihood::AbstractArray) -> WAICResult{<:NamedTuple,<:NamedTuple}Compute the widely applicable information criterion (WAIC).[Watanabe2010][Vehtari2017][LOOFAQ]
log_likelihood must be an array of log-likelihood values with shape (chains, draws[, params...]).
See also: WAICResult, loo
Examples
Calculate WAIC of a model:
julia> using ArviZExampleData
julia> idata = load_example_data("centered_eight");
@@ -88,7 +88,7 @@
julia> waic(log_like)
WAICResult with estimates
elpd elpd_mcse p p_mcse
- -31 1.4 0.9 0.33Model comparison
PosteriorStats.ModelComparisonResult — TypeModelComparisonResultResult of model comparison using ELPD.
This struct implements the Tables and TableTraits interfaces.
Each field returns a collection of the corresponding entry for each model:
name: Names of the models, if provided.rank: Ranks of the models (ordered by decreasing ELPD)elpd_diff: ELPD of a model subtracted from the largest ELPD of any modelelpd_diff_mcse: Monte Carlo standard error of the ELPD differenceweight: Model weights computed withweights_methodelpd_result:AbstactELPDResults for each model, which can be used to access useful stats like ELPD estimates, pointwise estimates, and Pareto shape values for PSIS-LOOweights_method: Method used to compute model weights withmodel_weights
PosteriorStats.compare — Functioncompare(models; kwargs...) -> ModelComparisonResultCompare models based on their expected log pointwise predictive density (ELPD).
The ELPD is estimated either by Pareto smoothed importance sampling leave-one-out cross-validation (LOO) or using the widely applicable information criterion (WAIC). We recommend loo. Read more theory here - in a paper by some of the leading authorities on model comparison dx.doi.org/10.1111/1467-9868.00353
Arguments
models: aTuple,NamedTuple, orAbstractVectorwhose values are eitherAbstractELPDResultentries or any argument toelpd_method.
Keywords
weights_method::AbstractModelWeightsMethod=Stacking(): the method to be used to weight the models. Seemodel_weightsfor detailselpd_method=loo: a method that computes anAbstractELPDResultfrom an argument inmodels.sort::Bool=true: Whether to sort models by decreasing ELPD.
Returns
ModelComparisonResult: A container for the model comparison results. The fields contain a similar collection tomodels.
Examples
Compare the centered and non centered models of the eight school problem using the defaults: loo and Stacking weights. A custom myloo method formates the inputs as expected by loo.
julia> using ArviZExampleData
+ -31 1.4 0.9 0.33Model comparison
PosteriorStats.ModelComparisonResult — TypeModelComparisonResultResult of model comparison using ELPD.
This struct implements the Tables and TableTraits interfaces.
Each field returns a collection of the corresponding entry for each model:
name: Names of the models, if provided.rank: Ranks of the models (ordered by decreasing ELPD)elpd_diff: ELPD of a model subtracted from the largest ELPD of any modelelpd_diff_mcse: Monte Carlo standard error of the ELPD differenceweight: Model weights computed withweights_methodelpd_result:AbstactELPDResults for each model, which can be used to access useful stats like ELPD estimates, pointwise estimates, and Pareto shape values for PSIS-LOOweights_method: Method used to compute model weights withmodel_weights
PosteriorStats.compare — Functioncompare(models; kwargs...) -> ModelComparisonResultCompare models based on their expected log pointwise predictive density (ELPD).
The ELPD is estimated either by Pareto smoothed importance sampling leave-one-out cross-validation (LOO) or using the widely applicable information criterion (WAIC). We recommend loo. Read more theory here - in a paper by some of the leading authorities on model comparison dx.doi.org/10.1111/1467-9868.00353
Arguments
models: aTuple,NamedTuple, orAbstractVectorwhose values are eitherAbstractELPDResultentries or any argument toelpd_method.
Keywords
weights_method::AbstractModelWeightsMethod=Stacking(): the method to be used to weight the models. Seemodel_weightsfor detailselpd_method=loo: a method that computes anAbstractELPDResultfrom an argument inmodels.sort::Bool=true: Whether to sort models by decreasing ELPD.
Returns
ModelComparisonResult: A container for the model comparison results. The fields contain a similar collection tomodels.
Examples
Compare the centered and non centered models of the eight school problem using the defaults: loo and Stacking weights. A custom myloo method formates the inputs as expected by loo.
julia> using ArviZExampleData
julia> models = (
centered=load_example_data("centered_eight"),
@@ -118,7 +118,7 @@
rank elpd elpd_mcse elpd_diff elpd_diff_mcse weight p ⋯
non_centered 1 -31 1.4 0 0.0 0.52 0.9 ⋯
centered 2 -31 1.4 0.06 0.067 0.48 0.9 ⋯
- 1 column omittedPosteriorStats.model_weights — Functionmodel_weights(elpd_results; method=Stacking())
+ 1 column omittedPosteriorStats.model_weights — Functionmodel_weights(elpd_results; method=Stacking())
model_weights(method::AbstractModelWeightsMethod, elpd_results)Compute weights for each model in elpd_results using method.
elpd_results is a Tuple, NamedTuple, or AbstractVector with AbstractELPDResult entries. The weights are returned in the same type of collection.
Stacking is the recommended approach, as it performs well even when the true data generating process is not included among the candidate models. See [YaoVehtari2018] for details.
See also: AbstractModelWeightsMethod, compare
Examples
Compute Stacking weights for two models:
julia> using ArviZExampleData
julia> models = (
@@ -139,10 +139,10 @@
:non_centered => 1.0Now we compute BootstrappedPseudoBMA weights for the same models:
julia> model_weights(elpd_results; method=BootstrappedPseudoBMA()) |> pairs
pairs(::NamedTuple) with 2 entries:
:centered => 0.483727
- :non_centered => 0.516273The following model weighting methods are available
PosteriorStats.AbstractModelWeightsMethod — Typeabstract type AbstractModelWeightsMethodAn abstract type representing methods for computing model weights.
Subtypes implement model_weights(method, elpd_results).
PosteriorStats.BootstrappedPseudoBMA — Typestruct BootstrappedPseudoBMA{R<:Random.AbstractRNG, T<:Real} <: AbstractModelWeightsMethodModel weighting method using pseudo Bayesian Model Averaging using Akaike-type weighting with the Bayesian bootstrap (pseudo-BMA+)[YaoVehtari2018].
The Bayesian bootstrap stabilizes the model weights.
BootstrappedPseudoBMA(; rng=Random.default_rng(), samples=1_000, alpha=1)
-BootstrappedPseudoBMA(rng, samples, alpha)Construct the method.
rng::Random.AbstractRNG: The random number generator to use for the Bayesian bootstrapsamples::Int64: The number of samples to draw for bootstrappingalpha::Real: The shape parameter in the Dirichlet distribution used for the Bayesian bootstrap. The default (1) corresponds to a uniform distribution on the simplex.
See also: Stacking
PosteriorStats.PseudoBMA — Typestruct PseudoBMA <: AbstractModelWeightsMethodModel weighting method using pseudo Bayesian Model Averaging (pseudo-BMA) and Akaike-type weighting.
PseudoBMA(; regularize=false)
-PseudoBMA(regularize)Construct the method with optional regularization of the weights using the standard error of the ELPD estimate.
This approach is not recommended, as it produces unstable weight estimates. It is recommended to instead use BootstrappedPseudoBMA to stabilize the weights or Stacking. For details, see [YaoVehtari2018].
See also: Stacking
PosteriorStats.Stacking — Typestruct Stacking{O<:Optim.AbstractOptimizer} <: AbstractModelWeightsMethodModel weighting using stacking of predictive distributions[YaoVehtari2018].
Stacking(; optimizer=Optim.LBFGS(), options=Optim.Options()
-Stacking(optimizer[, options])Construct the method, optionally customizing the optimization.
optimizer::Optim.AbstractOptimizer: The optimizer to use for the optimization of the weights. The optimizer must support projected gradient optimization via amanifoldfield.options::Optim.Options: The Optim options to use for the optimization of the weights.
See also: BootstrappedPseudoBMA
Predictive checks
PosteriorStats.loo_pit — Functionloo_pit(y, y_pred, log_weights; kwargs...) -> Union{Real,AbstractArray}Compute leave-one-out probability integral transform (LOO-PIT) checks.
Arguments
y: array of observations with shape(params...,)y_pred: array of posterior predictive samples with shape(draws, chains, params...).log_weights: array of normalized log LOO importance weights with shape(draws, chains, params...).
Keywords
is_discrete: If not provided, then it is set totrueiff elements ofyandy_predare all integer-valued. Iftrue, then data are smoothed usingsmooth_datato make them non-discrete before estimating LOO-PIT values.kwargs: Remaining keywords are forwarded tosmooth_dataif data is discrete.
Returns
pitvals: LOO-PIT values with same size asy. Ifyis a scalar, thenpitvalsis a scalar.
LOO-PIT is a marginal posterior predictive check. If $y_{-i}$ is the array $y$ of observations with the $i$th observation left out, and $y_i^*$ is a posterior prediction of the $i$th observation, then the LOO-PIT value for the $i$th observation is defined as
\[P(y_i^* \le y_i \mid y_{-i}) = \int_{-\infty}^{y_i} p(y_i^* \mid y_{-i}) \mathrm{d} y_i^*\]
The LOO posterior predictions and the corresponding observations should have similar distributions, so if conditional predictive distributions are well-calibrated, then all LOO-PIT values should be approximately uniformly distributed on $[0, 1]$.[Gabry2019]
Examples
Calculate LOO-PIT values using as test quantity the observed values themselves.
julia> using ArviZExampleData
+ :non_centered => 0.516273The following model weighting methods are available
PosteriorStats.AbstractModelWeightsMethod — Typeabstract type AbstractModelWeightsMethodAn abstract type representing methods for computing model weights.
Subtypes implement model_weights(method, elpd_results).
PosteriorStats.BootstrappedPseudoBMA — Typestruct BootstrappedPseudoBMA{R<:Random.AbstractRNG, T<:Real} <: AbstractModelWeightsMethodModel weighting method using pseudo Bayesian Model Averaging using Akaike-type weighting with the Bayesian bootstrap (pseudo-BMA+)[YaoVehtari2018].
The Bayesian bootstrap stabilizes the model weights.
BootstrappedPseudoBMA(; rng=Random.default_rng(), samples=1_000, alpha=1)
+BootstrappedPseudoBMA(rng, samples, alpha)Construct the method.
rng::Random.AbstractRNG: The random number generator to use for the Bayesian bootstrapsamples::Int64: The number of samples to draw for bootstrappingalpha::Real: The shape parameter in the Dirichlet distribution used for the Bayesian bootstrap. The default (1) corresponds to a uniform distribution on the simplex.
See also: Stacking
PosteriorStats.PseudoBMA — Typestruct PseudoBMA <: AbstractModelWeightsMethodModel weighting method using pseudo Bayesian Model Averaging (pseudo-BMA) and Akaike-type weighting.
PseudoBMA(; regularize=false)
+PseudoBMA(regularize)Construct the method with optional regularization of the weights using the standard error of the ELPD estimate.
This approach is not recommended, as it produces unstable weight estimates. It is recommended to instead use BootstrappedPseudoBMA to stabilize the weights or Stacking. For details, see [YaoVehtari2018].
See also: Stacking
PosteriorStats.Stacking — Typestruct Stacking{O<:Optim.AbstractOptimizer} <: AbstractModelWeightsMethodModel weighting using stacking of predictive distributions[YaoVehtari2018].
Stacking(; optimizer=Optim.LBFGS(), options=Optim.Options()
+Stacking(optimizer[, options])Construct the method, optionally customizing the optimization.
optimizer::Optim.AbstractOptimizer: The optimizer to use for the optimization of the weights. The optimizer must support projected gradient optimization via amanifoldfield.options::Optim.Options: The Optim options to use for the optimization of the weights.
See also: BootstrappedPseudoBMA
Predictive checks
PosteriorStats.loo_pit — Functionloo_pit(y, y_pred, log_weights; kwargs...) -> Union{Real,AbstractArray}Compute leave-one-out probability integral transform (LOO-PIT) checks.
Arguments
y: array of observations with shape(params...,)y_pred: array of posterior predictive samples with shape(draws, chains, params...).log_weights: array of normalized log LOO importance weights with shape(draws, chains, params...).
Keywords
is_discrete: If not provided, then it is set totrueiff elements ofyandy_predare all integer-valued. Iftrue, then data are smoothed usingsmooth_datato make them non-discrete before estimating LOO-PIT values.kwargs: Remaining keywords are forwarded tosmooth_dataif data is discrete.
Returns
pitvals: LOO-PIT values with same size asy. Ifyis a scalar, thenpitvalsis a scalar.
LOO-PIT is a marginal posterior predictive check. If $y_{-i}$ is the array $y$ of observations with the $i$th observation left out, and $y_i^*$ is a posterior prediction of the $i$th observation, then the LOO-PIT value for the $i$th observation is defined as
\[P(y_i^* \le y_i \mid y_{-i}) = \int_{-\infty}^{y_i} p(y_i^* \mid y_{-i}) \mathrm{d} y_i^*\]
The LOO posterior predictions and the corresponding observations should have similar distributions, so if conditional predictive distributions are well-calibrated, then all LOO-PIT values should be approximately uniformly distributed on $[0, 1]$.[Gabry2019]
Examples
Calculate LOO-PIT values using as test quantity the observed values themselves.
julia> using ArviZExampleData
julia> idata = load_example_data("centered_eight");
@@ -182,7 +182,7 @@
"Hotchkiss" 0.435094
"Lawrenceville" 0.220627
"St. Paul's" 0.775086
- "Mt. Hermon" 0.296706PosteriorStats.r2_score — Functionr2_score(y_true::AbstractVector, y_pred::AbstractArray) -> (; r2, r2_std)$R²$ for linear Bayesian regression models.[GelmanGoodrich2019]
Arguments
y_true: Observed data of lengthnoutputsy_pred: Predicted data with size(ndraws[, nchains], noutputs)
Examples
julia> using ArviZExampleData
+ "Mt. Hermon" 0.296706PosteriorStats.r2_score — Functionr2_score(y_true::AbstractVector, y_pred::AbstractArray) -> (; r2, r2_std)$R²$ for linear Bayesian regression models.[GelmanGoodrich2019]
Arguments
y_true: Observed data of lengthnoutputsy_pred: Predicted data with size(ndraws[, nchains], noutputs)
Examples
julia> using ArviZExampleData
julia> idata = load_example_data("regression1d");
@@ -193,4 +193,4 @@
julia> r2_score(y_true, y_pred) |> pairs
pairs(::NamedTuple) with 2 entries:
:r2 => 0.683197
- :r2_std => 0.0368838Utilities
PosteriorStats.kde_reflected — Functionkde_reflected(data::AbstractVector{<:Real}; bounds=extrema(data), kwargs...)Compute the boundary-corrected kernel density estimate (KDE) of data using reflection.
For $x \in (l, u)$, the reflected KDE has the density
\[\hat{f}_R(x) = \hat{f}(x) + \hat{f}(2l - x) + \hat{f}(2u - x),\]
where $\hat{f}$ is the usual KDE of data. This is equivalent to augmenting the original data with 2 additional copies of the data reflected around each bound, computing the usual KDE, trimming the KDE to the bounds, and renormalizing.
Any non-finite bounds are ignored. Remaining kwargs are passed to KernelDensity.kde. The default bandwidth is estimated using the Improved Sheather-Jones (ISJ) method [Botev2010].
PosteriorStats.smooth_data — Functionsmooth_data(y; dims=:, interp_method=CubicSpline, offset_frac=0.01)Smooth y along dims using interp_method.
interp_method is a 2-argument callabale that takes the arguments y and x and returns a DataInterpolations.jl interpolation method, defaulting to a cubic spline interpolator.
offset_frac is the fraction of the length of y to use as an offset when interpolating.
- Hyndman1996Rob J. Hyndman (1996) Computing and Graphing Highest Density Regions, Amer. Stat., 50(2): 120-6. DOI: 10.1080/00031305.1996.10474359 jstor.
- ChenShao1999Ming-Hui Chen & Qi-Man Shao (1999) Monte Carlo Estimation of Bayesian Credible and HPD Intervals, J Comput. Graph. Stat., 8:1, 69-92. DOI: 10.1080/10618600.1999.10474802 jstor.
- Vehtari2017Vehtari, A., Gelman, A. & Gabry, J. Practical Bayesian model evaluation using leave-one-out cross-validation and WAIC. Stat Comput 27, 1413–1432 (2017). doi: 10.1007/s11222-016-9696-4 arXiv: 1507.04544
- LOOFAQAki Vehtari. Cross-validation FAQ. https://mc-stan.org/loo/articles/online-only/faq.html
- Watanabe2010Watanabe, S. Asymptotic Equivalence of Bayes Cross Validation and Widely Applicable Information Criterion in Singular Learning Theory. 11(116):3571−3594, 2010. https://jmlr.csail.mit.edu/papers/v11/watanabe10a.html
- Vehtari2017Vehtari, A., Gelman, A. & Gabry, J. Practical Bayesian model evaluation using leave-one-out cross-validation and WAIC. Stat Comput 27, 1413–1432 (2017). doi: 10.1007/s11222-016-9696-4 arXiv: 1507.04544
- LOOFAQAki Vehtari. Cross-validation FAQ. https://mc-stan.org/loo/articles/online-only/faq.html
- YaoVehtari2018Yuling Yao, Aki Vehtari, Daniel Simpson, and Andrew Gelman. Using Stacking to Average Bayesian Predictive Distributions. 2018. Bayesian Analysis. 13, 3, 917–1007. doi: 10.1214/17-BA1091 arXiv: 1704.02030
- YaoVehtari2018Yuling Yao, Aki Vehtari, Daniel Simpson, and Andrew Gelman. Using Stacking to Average Bayesian Predictive Distributions. 2018. Bayesian Analysis. 13, 3, 917–1007. doi: 10.1214/17-BA1091 arXiv: 1704.02030
- YaoVehtari2018Yuling Yao, Aki Vehtari, Daniel Simpson, and Andrew Gelman. Using Stacking to Average Bayesian Predictive Distributions. 2018. Bayesian Analysis. 13, 3, 917–1007. doi: 10.1214/17-BA1091 arXiv: 1704.02030
- YaoVehtari2018Yuling Yao, Aki Vehtari, Daniel Simpson, and Andrew Gelman. Using Stacking to Average Bayesian Predictive Distributions. 2018. Bayesian Analysis. 13, 3, 917–1007. doi: 10.1214/17-BA1091 arXiv: 1704.02030
- Gabry2019Gabry, J., Simpson, D., Vehtari, A., Betancourt, M. & Gelman, A. Visualization in Bayesian Workflow. J. R. Stat. Soc. Ser. A Stat. Soc. 182, 389–402 (2019). doi: 10.1111/rssa.12378 arXiv: 1709.01449
- GelmanGoodrich2019Andrew Gelman, Ben Goodrich, Jonah Gabry & Aki Vehtari (2019) R-squared for Bayesian Regression Models, The American Statistician, 73:3, 307-9, DOI: 10.1080/00031305.2018.1549100.
- Botev2010Kernel density estimation via diffusion. Z. I. Botev, J. F. Grotowski, and D. P. Kroese (2010) Annals of Statistics, Volume 38, Number 5, pages 2916-2957. doi: 10.1214/10-AOS799