tensiometer.chains_convergence

This file contains some functions to study convergence of the chains and to compare the two posteriors.

tensiometer.chains_convergence.GR_test(chains, param_names=None)

Function performing the Gelman Rubin (GR) test (described in Gelman and Rubin 92 and Brooks and Gelman 98) on a list of MCSamples or on a single MCSamples with different sub-chains. This test compares the variation of the mean across a pool of chains with the expected variation of the mean under the pdf that is being sampled. If we define the covariance of the mean as:

\[C_{ij} \equiv {\rm Cov}_c({\rm Mean}_s(\theta))_{ij}\]

and the mean covariance as:

\[M_{ij} \equiv {\rm Mean}_c[{\rm Cov}_s(\theta)_{ij}]\]

then we seek to maximize:

\[R-1 = {\rm max_{\theta}}\frac{C_{ij} \theta^i \theta^j} {M_{ij}\theta^i \theta^j}\]

where the subscript \(c\) means that the statistics is computed across chains while the subscrit \(s\) indicates that it is computed across samples. In this case the maximization is solved by finding the maximum eigenvalue of \(C M^{-1}\).

Parameters:

• chains – single or list of MCSamples
• param_names – names of the parameters involved in the test. By default uses all non-derived parameters.

Returns:

value of the GR test and corresponding parameter combination

tensiometer.chains_convergence.GR_test_from_samples(samples, weights)

Lower level function to perform the Gelman Rubin (GR) test. This works on a list of samples from different chains and corresponding weights. Refer to tensiometer.chains_convergence.GR_test() for more details of what this function is doing.

Parameters:

• samples – list of samples from different chains
• weights – weights of the samples for each chain

Returns:

value of the GR test and corresponding parameter combination

tensiometer.chains_convergence.GRn_test(chains, n, theta0=None, param_names=None, feedback=0, optimizer='ParticleSwarm', **kwargs)

Multi dimensional higher order moments convergence test. Compares the variation of a given moment among the population of chains with the expected variation of that quantity from the samples pdf.

We first build the \(k\) order tensor of parameter differences around a point \(\tilde{\theta}\):

\[Q^{(k)} \equiv Q_{i_1, \dots, i_k} \equiv (\theta_{i_1} -\tilde{\theta}_{i_1}) \cdots (\theta_{i_k} -\tilde{\theta}_{i_k})\]

then we build the tensor encoding its covariance across chains

\[V_M = {\rm Var}_c (E_s [Q^{(k)}])\]

which is a \(2k\) and then build the second tensor encoding the mean in chain moment:

\[M_V = {\rm Mean}_c (E_s[(Q^{(k)}-E_s[Q^{(k)}]) \otimes(Q^{(k)}-E_s[Q^{(k)}])])\]

where we have suppressed all indexes to not crowd the notation.

Then we maximize over parameters:

\[R_n -1 \equiv {\rm max}_\theta \frac{V_M \theta^{2k}}{M_V \theta^{2k}}\]

where \(\theta^{2k}\) is the tensor product of \(\theta\) for \(2k\) times.

Differently from the 2D case this problem has no solution in terms of eigenvalues of tensors so far and the solution is obtained by numerical minimization with the pymanopt library.

Parameters:

• chains – single or list of MCSamples
• n – order of the moment
• theta0 – center of the moments. By default equal to the mean
• param_names – names of the parameters involved in the test. By default uses all non-derived parameters.
• feedback – level of feedback. 0=no feedback, >0 increasingly chatty
• optimizer – choice of optimization algorithm for pymanopt. Default is ParticleSwarm, other possibility is TrustRegions.
• kwargs – keyword arguments for the optimizer.

Returns:

value of the GR moment test and corresponding parameter combination

tensiometer.chains_convergence.GRn_test_1D(chains, n, param_name, theta0=None)

One dimensional higher moments test. Compares the variation of a given moment among the population of chains with the expected variation of that quantity from the samples pdf.

This test is defined by:

\[R_n(\theta_0)-1 = \frac{{\rm Var}_c ({\rm Mean}_s(\theta-\theta_0)^n)}{{\rm Mean}_c ({\rm Var}_s(\theta-\theta_0)^n) }\]

where the subscript \(c\) means that the statistics is computed across chains while the subscrit \(s\) indicates that it is computed across samples.

Parameters:

• chains – single or list of MCSamples
• n – order of the moment
• param_name – names of the parameter involved in the test.
• theta0 – center of the moments. By default equal to the mean.

Returns:

value of the GR moment test and corresponding parameter combination (an array with one since this works in 1D)

tensiometer.chains_convergence.GRn_test_1D_samples(samples, weights, n, theta0=None)

Lower level function to compute the one dimensional higher moments test. This works on a list of samples from different chains and corresponding weights. Refer to tensiometer.chains_convergence.GRn_test_1D() for more details of what this function is doing.

Parameters:

• samples – list of samples from different chains
• weights – weights of the samples for each chain
• n – order of the moment
• theta0 – center of the moments. By default equal to the mean.

Returns:

value of the GR moment test and corresponding parameter combination (an array with one since this works in 1D)

tensiometer.chains_convergence.GRn_test_from_samples(samples, weights, n, theta0=None, feedback=0, optimizer='ParticleSwarm', **kwargs)

Lower level function to compute the multi dimensional higher moments test. This works on a list of samples from different chains and corresponding weights. Refer to tensiometer.chains_convergence.GRn_test() for more details of what this function is doing.

Parameters:

• samples – list of samples from different chains
• weights – weights of the samples for each chain
• n – order of the moment
• theta0 – center of the moments. By default equal to the mean
• feedback – level of feedback. 0=no feedback, >0 increasingly chatty
• optimizer – choice of optimization algorithm for pymanopt. Default is ParticleSwarm, other possibility is TrustRegions.
• kwargs – keyword arguments for the optimizer.

Returns:

value of the GR moment test and corresponding parameter combination