Analyzing MCMC Output
Dootika Vats Nathan Robertson James M Flegal Galin L Jones
stat.CO stat.ME


Synchrotron x-ray scattering of UN and U_2N_3 epitaxial films
    G. H. Lander
    December 30, 2023
=============================================================







  §.§.§ Article Type:

Overview





  §.§.§ Abstract



  Markov chain Monte Carlo (MCMC) is a sampling-based method for
  estimating features of probability distributions. MCMC methods produce a serially correlated, yet representative, sample from
  the desired distribution.  As such it can be difficult to assess when
  the MCMC method is producing reliable results.  We present some
  fundamental methods for ensuring a reliable simulation
  experiment.  In particular, we present a workflow for output
  analysis in MCMC providing estimators, approximate sampling
  distributions, stopping rules, and visualization tools.
















§ INTRODUCTION
 


Markov chain Monte Carlo (MCMC) algorithms are an essential tool for
estimating features of probability distributions encountered in diverse
applications <cit.>.  The use of MCMC is commonly
identified with Bayesian settings, but it is also useful in other situations
<cit.>. 

Suppose, for our application, we have developed a probability
distribution F with support 𝒳⊆ℝ^d,
d ≥ 1.  Our goal is to use fixed, unknown features of F to make
inference about the population. For example, for
h : 𝒳→ℝ, we may be interested in the
expectation[If f is a probability function, then this notation avoids
  having separate formulas for the continuous case where
  μ_h = ∫_𝒳 h(x) f(x) dx and the discrete case
  where μ_h = ∑_x ∈𝒳 h(x) f(x) but it also
  applies to more general settings.]

    μ_h = _F[h(X)] = ∫_𝒳 h(x) F(dx)   .


The apparent simplicity of this hides an array of features, including
probabilities, means, moments, and marginal densities associated with
F.  Accordingly, we will typically want to use several expectations.
There are features of F that are not expectations, such as
quantiles, but we will defer discussion of this to
Section <ref>.  We collect all the features
of F we want in a p-dimensional vector, θ.

We use the notation ∼ to mean “distributed as,”
·∼ to mean “approximately
distributed as,” ≁ to mean “not distributed as,” and
≈ to mean “is approximately equal to.” Often F is
analytically intractable in the sense that directly calculating
θ is impossible and hence we may turn to estimating θ
using Monte Carlo sampling methods.  We consider only MCMC in which a
realization of a (aperiodic, Harris recurrent – for definitions see
<cit.>) Markov chain {X_1, X_2, …, X_n} is
produce so that for a sufficiently large Monte Carlo sample size n,
we have X_n ·∼ F. We will not
discuss how to construct or implement MCMC methods
<cit.>, but will instead assume that there is
already an efficient method for producing the MCMC sample, or output.
The output is then used to construct estimators of θ so that
for a sufficiently large Monte Carlo sample size, n, we have

    θ̂_n = θ̂(X_1, …, X_n) ≈θ .

In typical implementations of MCMC, the initial state, X_1,  is drawn from some other distribution than the target, that is, X_1  F, which
implies X_n  F for any finite Monte Carlo sample size n.
Moreover, there is an inherent serial correlation in the MCMC
sample. That is, the Markov chain draws are neither independent nor
identically distributed.  Thus, there are two common tasks in MCMC
output analysis (i) deciding when the simulation has produced a
representative sample from F, that is, when is n large enough so
that X_n ·∼ F, and (ii) when is n
sufficiently large to conclude θ̂_n≈θ.
Notice that the required number of draws may be different for each
task.  Task (i) is difficult and while there are some rigorous
approaches <cit.>, these are typically challenging to
implement in practically relevant settings.  This has led most
practitioners to approach this question by relying on graphical
summaries and ad hoc convergence diagnostics; see
Section <ref>.  Task (ii) is our main focus and is
fundamental to ensuring a reliable simulation experiment.  Classical
large-sample frequentist methods provide principled, practical
solutions for (ii).  However, since it is a Markov chain that is being
simulated, specialized techniques are required for their
implementation; this is covered in some detail in
Sections <ref>–<ref>.



In practice, addressing tasks (i) and (ii) can seem complicated to the
uninitiated.  Thus, we present a workflow for analyzing output from
MCMC addressing these challenges.  This is illustrated in the context of a
Bayesian example in Section <ref>. 



§ MARKOV CHAIN THEORY AND STARTING VALUES
 

The dynamics of a time-homogeneous Markov chain are dictated by its
Markov transition kernel

    P^n(x, A) := (X_n+j∈ A | X_j =x) ,       n,j≥ 1

and P^1 ≡ P. If ν is the initial distribution (i.e.X_1 ∼ν), then ν P^n is the marginal distribution of
X_n. Markov chains for MCMC are constructed in such a way that the
target distribution F, is its stationary distribution. That is,
FP^n = F, so that if X_1 ∼ F, then each X_n ∼ F.  Of
course, producing X_1 ∼ F often is not possible in settings where
MCMC is relevant. However, if · denotes total variation
distance, under standard conditions <cit.>,

    ν P^n - F → 0  as n →∞.

This implies X_n d→ F as n →∞.  Moreover,
the total variation norm in (<ref>) is non-increasing in n
<cit.>.  Hence a representative, although correlated,
sample will be produced eventually, that is, for sufficiently large n.

Ideally, we would like to identify n^⋆ < ∞ such that if
n ≥ n^⋆, then ν P^n - F  < ϵ for some
ϵ >0 so that we could conclude
X_n ·∼ F.  Indeed, there are
methods for doing so <cit.>, but they are
often so conservative as to be of little practical use or are
difficult to apply <cit.>.  This has forced
practitioners to turn to other methods, the most common of which are
trace plots of the components of the simulation and so-called
convergence diagnostics <cit.>.  Perhaps the most
commonly used convergence diagnostic was developed by
<cit.>.  However, this diagnostic has been shown to
have severe limitations <cit.>;
see Section <ref> for more. In fact, many
convergence diagnostics, including the Gelman-Rubin diagnostic, were
developed to address task (ii) in Section <ref>, but
are often incorrectly understood to answer task (i).  In general, all
convergence diagnostics should be used with care since their very use
can introduce bias <cit.>.  Strictly speaking,
convergence does not occur for any finite n so all they can detect
is evidence of non-convergence.  Hence diagnostics can never guarantee
what we want because absence of evidence of non-convergence is not
evidence of convergence. For more discussion on convergence
diagnostics see <cit.>.

While the convergence in (<ref>) holds for any starting value,
some will be better than others since the rate of convergence can be
affected by the choice of starting value <cit.>. In
particular, starting in an area of low probability of F can lead to
slow convergence of the Markov chain <cit.>. If,
however, X_1∼ F, then FP^n - F = 0 for all n, and the
Markov chain produces exact draws from F (albeit still
correlated). Thus, starting values can have a substantial impact on
the quality of the samples. Choosing good starting values can save
considerable time in postprocessing and increase confidence in the
results. 



It is a truism that “Any point you don't mind having in a sample is a
good starting point.” <cit.>. While choosing a good
starting value may be difficult, it may not be impossible. In fact, it
may be possible to start from stationarity via perfect simulation
<cit.>, Bernoulli factories <cit.>, or
simple accept-reject samplers. When the target distribution is
low-dimensional or made low-dimensional by a linchpin variable trick
<cit.>, an accept-reject sampler that can produce one 
sample from the target may be available. Such a sampler is often too
computationally burdensome for a full Monte Carlo procedure, but may
provide a single draw at reasonable cost.

Starting values are also often obtained by finding a high probability
region via optimization.  When a closed-form expression is available, using the maximum
likelihood estimate can be computationally cheap.  Other optimization approaches are  certainly possible, but
finding a global optimum is in general difficult.  Even so, in many
situations, any value in a high probability region is a reasonable
starting value.  In Bayesian settings, practitioners may draw starting
values from the (proper) prior distributions of the parameters.



This can work particularly well when the prior distributions have been
chosen with care.  Finally, and particularly when implementing
component-wise MCMC methods <cit.> such as Gibbs
samplers, the parameter being updated first may not require a starting
value. So it is often a good idea to first update the parameter whose
starting value is least trustworthy.




§ ESTIMATION AND SAMPLING DISTRIBUTIONS
 
 Recall that given a realization
{X_1, X_2, …, X_n} of the Markov chain we construct an
estimator of θ,
θ̂_n = θ̂(X_1, …, X_n) so that
θ̂_n→θ almost surely as n →∞.
However, no matter how large the (finite) Monte Carlo sample size n,
there will be an unknown Monte Carlo error
θ̂_n - θ.  The approximate sampling distribution of
the Monte Carlo error is often available through a version of the
central limit theorem (CLT), which holds under moment conditions on
the functionals and Markov chain mixing conditions <cit.>,
both of which require theoretical study to verify in a given
application.  <cit.> give an accessible introduction to
this theory which has been applied in a number of practically
relevant settings <cit.>



 §.§ Means
 

In most settings we will want to estimate several expectations.  Let
h: 𝒳→ℝ^p be a function such that
_Fh(X)= μ_h is of interest. For example, to estimate the mean
vector of F, h will be the identity function. Estimation is
straightforward using a
sample mean since the Markov chain strong law ensures

    μ̅_n:= 1n∑_t = 1^n h(X_t) a.s.→μ_h  as n →∞ .

If a CLT holds, then there exists a p × p positive definite matrix,
Σ, such that as n→∞,

    √(n)(μ̅_n - μ_h) d→ N_p(0, Σ) .

Here, Σ encodes the covariance structure for h in the target distribution and the serial lag-covariance due to the Markov chain. More specifically,

    Σ = ∑_k=-∞^∞Cov_F(h(X_1), h(X_1+k)) .


The subscript F in (<ref>) means that the expectations
are calculated under the assumption that X_1 ∼ F.  This does
not mean that we need X_1 ∼ F for the CLT to hold.
Indeed, if the CLT holds for one initial distribution, then it holds
for every initial distribution <cit.>.

Under an independent sampling scheme,
Cov_F(h(X_1), h(X_1+k)) = 0 for all k  0, but the
Markov chains encountered in MCMC applications exhibit  serial
dependence. Thus, utilizing the sampling distribution for
μ̅_n to make large-sample inference requires specialized
methods for estimating Σ, which we discuss later.




 §.§ Quantiles
 

In addition to expectations, marginal quantiles are often of interest. Let
h:→ℝ, and for X ∼ F, set V = h(X). Let F_V be the
distribution of V, and for 0 < q < 1, the q-quantile of F_V is defined
as:

    ϕ_q = inf{x: F_V(x) ≥ q  } .

A natural estimator of ϕ_q is the ⌈ nq ⌉th order statistic. If {V_1,
V_2, …, V_n} is the transformed process, and {V_(1), V_(2), …,
V_(n)} are the order statistics, then an estimator of ϕ_q is

    ϕ̂_q:= V_(j+1) where  j < nq  ≤  j+1 ,

and ϕ̂_q a.s.→ϕ_q as n →∞.  An
approximate sampling distribution for ϕ̂_q is developed by
<cit.>. First, for any y define, 
    σ^2(y) =  ∑_k=-∞^∞Cov_F_V(I (V_1 ≤ y),   I (V_1+k≤ y)) ,

and let f_V be the density associated with F_V. Then, as n →∞,

    √(n) (ϕ̂_q - ϕ_q) d→N(0, σ^2(ϕ_q)/f_V(ϕ_q)^2) .

Here, we present a univariate sampling distribution for the quantiles, but
often multiple quantiles may be of interest. The joint sampling distributions
of multiple quantiles and of means and quantiles is available in
<cit.>.




 §.§ Functions of expectations
 


We are often interested in estimating functions of expectations, in
which case a plug-in type estimator is natural, and the approximate
sampling distribution may be obtained by an application of the delta
method. This includes estimating the covariance matrix of F,
Λ = _F(h(X_1)), with the sample covariance matrix,

    Λ_n := 1n∑_t = 1^n (h(X_t) -
    μ̅_n)(h(X_t) - μ̅_n)^T .

A delta method argument similar to the independent and identically
distributed setting yields an element-wise sampling distribution for
Λ_n.




§ ESTIMATING MONTE CARLO ERROR
 


The approximate sampling distributions of the previous section provide
the keys to assessing the reliability of the simulation effort, that
is, addressing task (ii) from
Section <ref>. Construction of confidence regions for
μ_h to address this problem has attracted substantial interest
<cit.>.  Suppose that Σ_n is an
estimator of Σ in the CLT (<ref>).  If
T^2_1-α, p , q denotes a 1-α quantile from a
Hotelling's T-squared distribution with dimensionality p and
degrees of freedom q, then it is straightforward to construct an
asymptotic 100(1-α)% confidence region

    C_α(n) = { n(μ̅_n - μ_h)^T Σ_n^-1
      (μ̅_n - μ_h) < T^2_1-α, p , q}  .

The size of C_α(n) will then describe the precision of
estimation; we will discuss how to use this information in the
sequel. Thus, the main obstacle is estimating the variance in the
asymptotic distribution. There are a variety of estimators available
that broadly fit into three classes (1) spectral variance estimators
<cit.>,
(2) batch means estimators <cit.>, and
(3) initial sequence estimators <cit.>. Of these, the most popular are the batch means
estimators since they are easy to implement and are computationally
efficient.

The multivariate batch means estimator considers non-overlapping
batches of the Markov chain and constructs a sample covariance matrix
from the sample means of each batch. More formally, let n = ab where
a is the number of batches and b is the batch size. For
k = 0, …, a-1, define
Y̅_k:= b^-1∑_t=1^b h(X_kb + t). The batch means
estimator of Σ is,

    Σ̂_n,b:= ba-1∑_k = 0^a-1(Y̅_k - μ̅_n) (Y̅_k - μ̅_n )^T .

The asymptotic behavior of batch means estimators has been well
studied, however small sample performance of batch means estimators
can be suspect in the presence of high correlation. Recently,
carefully constructed linear combinations of batch means estimators
have been proposed for improving finite sample performance of
estimators of Σ <cit.>. Particularly, for an
even b, they obtain the following flat-top batch means estimator

    Σ̃_n = 2 Σ̂_n,b - Σ̂_n,b/2 ,

where Σ̂_n,b/2 is the batch means estimator from a batch
size of b/2. 

Ensuring the batch means estimator (or its variants) is strongly and
mean-square consistent requires that the batch size b and the number
of batches a must be chosen so that both increase to infinity as
n →∞. The choice of batch size b is critical to the
performance of the batch means estimator <cit.>. In particular, <cit.> show that the
mean-squared-optimal batch size is b ∝ n^1/3, where the
proportionality constant needs to be estimated separately, and its
size depends on the amount of serial correlation in the
chain. <cit.> present a parametric method of
estimating this proportionality constant, yielding an easily
implementable optimal batch size.








§ STOPPING THE SIMULATION
 


The justification for using MCMC experiments to estimate features of
F is asymptotic, but, in practice, the Monte Carlo sample size n
is finite.  Thus, the choice of n is crucial to ensuring the
reliability of the MCMC experiment.  There are several approaches that
can be used to terminate the simulation.  The simplest is a fixed-time
procedure where the practitioner specifies the Monte Carlo sample size
before the experiment begins. In this case, we can estimate the Monte
Carlo error and, if it is too large, then run the simulation longer.
This leads to a sequential fixed-volume or fixed relative-volume
approach <cit.> which terminates simulation at a
random time.  Discussion of these issues in the context of MCMC and
many more details about them can be found in <cit.>,
<cit.>, <cit.>, and
<cit.>.

We again focus on estimating a vector of expectations for the sake of
specificity.  A fruitful approach is to compare the volume of the
confidence region, Vol(C_α(n)), to the generalized
variance of the target F, that is, |Λ| where we recall
Λ = Var_F(h(X_1)) and | · | denotes the determinant.
In doing this, <cit.> show that simulating until
the effective sample size (ESS), defined as

    ESS:= n (|Λ||Σ|)^1/p  ,

is sufficiently large implies that the Monte Carlo error is
sufficiently small compared to the variability in the target
distribution.  ESS provides the number of independent and identically
distributed samples that will yield
the same generalized Monte Carlo error as the correlated sample.  As
discussed above, estimators of both Λ and Σ are
available, so that we can easily estimate it with

    ESS:= n (|Λ_n||Σ_n|)^1/p .

Here Σ_n is any of the estimators of Σ discussed in Section <ref>. Notice that ESS depends on the function h. It is thus important to first
establish the quantities of interest before estimating ESS. When p = 1, the
ESS defined above is the same as the univariate ESS of 
<cit.>.

A natural question to ask is what ESS is sufficient for reliable estimation? <cit.> provide a principled cutoff for the ESS based on the relative quality of estimation. Suppose we are interested in making 100(1 - α)% confidence regions of μ_h using μ̅_n such that the volume of the confidence region is an ϵth fraction of the variability of h under F. Then simulation can stop the first time 

    ESS≥2^2/pπ(p Γ(p/2))^2/pχ^2_1-α, pϵ^2:= M_α, ϵ, p .

<cit.> showed that as ϵ→ 0, stopping simulation according to (<ref>) will yield confidence regions with coverage probabilities converging to 1 - α. In practice, we do not check criterion (<ref>) until after some minimum n^* steps to avoid premature termination due to poor early estimates of Σ and Λ. A reasonable choice of n^* is M_α, ϵ, p, which can be calculated a priori since all quantities are known.

The lower bound in (<ref>) is in the spirit of sample
size calculations for a desired half-width of a confidence interval
for a simple test of means, where ϵ serves the same purpose
as a (relative) half-width. The relative tolerance level, ϵ,
can be chosen according to the level of precision in the estimates
desired, typically ϵ≤ .05.  Once ESS
reaches the cutoff, users can stop simulating since they can be
100(1-α)% confident that estimates are within an ϵth level of tolerance relative to the target distribution.


The Gelman-Rubin-Brooks diagnostic of <cit.> and
<cit.> is a popular tool for assessing convergence
in task (ii). The original diagnostic estimates Σ, by the
sample covariance of mean vectors from m independent chains. Since
m is often small, the original Gelman-Rubin-Brooks statistic has
high variability <cit.>. However, a more robust statistic is achieved if the
estimators of Σ described in Section <ref>
are used in the Gelman-Rubin-Brooks diagnostic
<cit.>.
In addition, if R̂^p denotes the improved Gelman-Rubin-Brooks statistic, then <cit.> showed that 

    R̂^p ≈√(1 + 1ESS) .

Thus, there is a direct relationship between ESS and the
Gelman-Rubin-Brooks statistic.  The above is specifically for a single
chain, while a multiple chain statistic is provided in
<cit.>. 




The Gelman-Rubin-Brooks diagnostic suggests that the simulation be
terminated when R̂^p is below a pre-defined
cutoff. <cit.> used the bound in
(<ref>) to obtain a cutoff for R̂^p : simulation
stops the first time

    R̂^p ≤√(1 + 1M_α, ϵ, p) .

However, we recommend using ESS as it is more naturally interpretable
and has lower variability in readily available software.  Notice that
both R̂^p and ESS are aimed at assessing
estimation and not ensuring a representative sample; that is, aimed at
addressing task (ii) and not task (i) from
Section <ref>.





§ EXTENSIONS



The practice of thinning (or subsampling) the Markov chain is common.
Thinning a Markov chain refers to using every mth observation in the
simulation.  This is often done to reduce autocorrelation, but the
resulting process is still a Markov chain with kernel corresponding to
P^m. <cit.>, <cit.>, and
<cit.> showed that since this comes at the cost of the
number of usable samples, the variance of a thinned Monte Carlo
estimator is always larger than the variance of the original Monte
Carlo estimator.

There are, however, situations where thinning is worthwhile. <cit.>
and <cit.>  argue that when post-processing on the raw MCMC data is
expensive, it may be computationally efficient to thin the samples. That is,
if the function h is costly to evaluate relative to the time it takes to get
more samples, thinning is beneficial. Thinning can also be useful in
high-dimensional problems where storing the full MCMC output calls for large
memory requirements.  When the original MCMC sample is thinned, all output
analysis procedures then apply to the thinned MCMC sample.

Simulating multiple parallel chains is also common practice, and can
often be useful in parallel computing environments. However, multiple
short runs can be misleading and can retain large estimation bias
<cit.>, thus we encourage users to run
each chain for as long as they would with a single chain.





Another issue that has received little attention is estimation of higher order
moments. These fit naturally into the discussion in
Section <ref>, but bring new practical challenges
since, for the same Monte Carlo sample size, estimation quality reduces
drastically as moments increase.

Our discussion has focused on the time-homogeneous Markov chains
encountered in standard MCMC applications. The theoretical and
practical tools required for other simulation techniques can be quite
different. For adaptive MCMC, <cit.> provides estimators for
the asymptotic variance of the Monte Carlo
estimator. <cit.> provide uncertainty quantification
for trans-dimensional MCMC methods. However, the literature for these
processes is not as rich as traditional MCMC, and would benefit from
further work.

















§ EXAMPLE WORKFLOW
 


Using an example, we now illustrate a workflow for MCMC output
analysis integrated with useful visualizations. We use a Bayesian
reliability model to assess the reliability of liquid crystal display
(LCD) projectors. We recognize the target distribution F, establish
the function of interest h, choose an appropriate MCMC sampler,
choose starting values, visually assess the quality of the MCMC
sampler using a preliminary Monte Carlo sample, implement the stopping
rules for h, report point estimates, and provide appropriate
graphical tools.

Consider the LCD projector data of <cit.>. To test the manufacturer's claim of expected lamp life in an LCD projector being 1500 hours, identical lamps were placed in 31 projectors for various models and their time to failure was recorded. The data is presented in Table <ref>.


For i = 1, …, 31, let t_i denote the observed failure time for each lamp. <cit.> assumed that the t_i's are a realization from,

    T_i ∼Weibull(λ, β) ,
 
where λ > 0 is the scale parameter and β >0 is the shape parameter. Interest is in estimating the mean time to failure (MTTF) and the reliability function at t = 1500. Under the Weibull likelihood, the MTTF is

    MTTF = λ^-1/β Γ(1 + 1β) ,

and the reliability function is

    R(t) = exp{-λ t^β} .

<cit.> further assume priors
λ∼Gamma(2.5, 2350)  and β∼Gamma(1, 1), where each is parameterized by a shape and rate parameter.  The density of the posterior is

    f(λ, β| T) ∝λ^32.5β^31( ∏_i=1^31 t_i )^β - 1exp{ - λ∑_i=1^31 t_i^β}exp{-  β}exp{- 2350 λ}  ,
 
where the normalizing constant is unknown and we use component-wise MCMC methods <cit.> to sample from this distribution. The component λ will be updated by a Gibbs step and β will be updated by a Metropolis-Hastings step. The full conditional distribution of λ is

    λ|β, T ∼Gamma( 33.5, 2350 + ∑_i=1^31 t_i^β) .

The full conditional distribution for β is not available in closed-form, so we implement a Metropolis-Hastings step, yielding a Metropolis-within-Gibbs sampler <cit.>. The proposal distribution is a N( ·, .1^2), which yields an approximately optimal acceptance probability as suggested by <cit.>. We update λ first, followed by β, thus a starting value is only needed for β. We start from the MLE of β which is approximately 1.12.

We are interested in estimating MTTF and R(1500) to contest the
manufacturers' claims. Thus, the function whose posterior mean is of
interest h is

    h(λ, β) = ([ λ^-1/β Γ(1 + 1β);      exp{-λ t^β} ]) .


Here p = 2 and setting the relative tolerance to be
ϵ = .05 and α=.05 yields a minimum ESS of 7529 using
(<ref>). We first run the MCMC sampler for n^* = 7529  steps as
a check to see whether the sampler is exploring the state space
adequately and mixing well. Any issues with the sampler may be
addressed in such preliminary steps before a long run for a final
analysis is reported. Figure <ref> shows the trace plots
of λ and β, which indicates that there are no obvious
mixing problems.






This initial run yields ESS = 1196 which is
noticeably smaller than its cutoff. Seeing this, we run the sampler
for 92471 more steps yielding an overall sample size of 1e5. The final
ESS = 11834 > 7529 which indicates we are 95%
confident that the simulation has terminated with a relative tolerance
smaller than ϵ = .05.

Figure <ref> shows the autocorrelation in MTTF  and R(1500) along with a cross-correlation plot and the Monte Carlo confidence region for μ̅_n. MTTF has little autocorrelation and R(1500) has significant autocorrelation over lag 50. A multivariate analysis of the MCMC output is critically important in this example as seen by the cross-correlation plot. The cross-correlation plot at zero indicates the correlation between MTTF and R(1500) in the posterior distribution. This correlation along with the lag correlations would be completely ignored by a univariate analysis. The resulting confidence region constructed using the batch means estimator and the sampling distribution in (<ref>) captures this cross-correlation structure.



The final estimates of MTTF and R(1500) are 596.8 and 0.073
respectively. The MTTF of 596.8 is significantly far from 1500 with a
95% credible interval of (434, 834). The reliability function at a
failure time of 1500 is also low with a 95% credible interval of
(0.020, 0.163). The marginal density plots with the respective mean
and quantiles are given in Figure <ref>, along with
simultaneous uncertainty plotted for all 6 estimates using the methods
of <cit.>. The error bands around the estimates are
nearly indistinguishable from the estimates themselves, indicating a
small Monte Carlo error.



apalike


