public class PotentialScaleReduction extends Object
PotentialScaleReduction
class provides an online
computationa of Rhat, the potential scale reduction statistic for
measuring mixing and convergence of multiple Markov chain Monte
Carlo (MCMC) samplers.
At construction time, the number of estimators is specified.
There must be at least two estimators in order to compute Rhat.
Samples from the Markov chains are provided to this class via the
update(int,double)
method.
M
Markov chains with N
samples
each, with a sample being a floating point value y[m,n]
.
As usual, y[m,]
is the sequence of samples from the single
chain m
and y[,n]
are the n
th samples from
each chain. Unbiased mean and variance estimates are defined for
sequences in the usual way (see OnlineNormalEstimator
for
definitions). Using vector notation, mean'(y[,])
is the
average value of all samples whereas mean'(y[m,])
is the
average of samples from chain m
; similarly var'(y[,])
is the variance over all samples and var'(y[m,])
the variance of samples in chain m
.
The definition of the Rhat is:
whereRhat = sqrt(varHatPlus/W)
varHatPlus
is a weighted average of the withinchain
(W
) and betweenchain (B
) variances.
The betweenchain variance is defined byvarHatPlus = (N1)/N * W + 1/N * B.
The withinchain variance is the average of the unbiased withinchain variance estimates:B = N * var'(mean'(y[m,])) = N/(M1) * Σ_{m} (mean'(y[m,])  mean'(y[,]))^{2}.
This is the usual definition for chains in which there are the same number of samples. For the implementation here, we takeW = mean'(var'(y[m,])) = 1/M Σ_{m} var'(y[m,]).
N
to be the minimum of the numbers of samples in the chains. The
withinchain statistics mean'(y[m,])
and var'(y[m,])
are computed using all of the samples for chain m
. But the crosschain statistics are not normalized, so mean'(y[,])
is computed here as mean'(mean'(y[m,]))
.
mean'(y[m,])
and var'(y[m,])
, are available through the estimator returned by
estimator(int)
.
An estimator for the complete set of samples mean and variance,
mean'(y[,])
and var'(y[,])
, are available through
globalEstimator()
. Note that these are truly global
estimates, not the estimates used in asynchronous Rhat calculations
as defined in he previous section.
Constructor and Description 

PotentialScaleReduction(double[][] yss)
Construct a potential scale reduction for the specified matrix
Of estimates for each chain.

PotentialScaleReduction(int numChains)
Construct a potential scale reduction with the specified number
of Markov chains for input.

Modifier and Type  Method and Description 

OnlineNormalEstimator 
estimator(int chain)
Returns the estimator for the specified chain.

OnlineNormalEstimator 
globalEstimator()
Returns the estimator that pools the estimates across the
chains.

int 
numChains()
Returns the number of chains for this estimator.

double 
rHat()
Returns the Rhat statistic as defined in the class
documentation.

void 
update(int chain,
double y)
Provide a sample of the specified value from the specified chain.

public PotentialScaleReduction(int numChains)
numChains
 Number of Markov chains.IllegalArgumentException
 If the number of chains is less than 2.public PotentialScaleReduction(double[][] yss)
yss[m][n]
are for the n
th sample from chain m
. The chains
may have different numbers of samples.yss
 Matrix of estimates by chain and sample.IllegalStateException
 If the number of chains (length
of yss
) is less than 2.public int numChains()
public OnlineNormalEstimator estimator(int chain)
chain
 Index of chain.public OnlineNormalEstimator globalEstimator()
public void update(int chain, double y)
chain
 Chain from which sample was drawn.y
 Value of sample.IndexOutOfBoundsException
 If the chain is less than zero or
greater than or equal to the number of chains.public double rHat()