Analysis

We don’t only want to sample field configurations from the distribution given by the action, we want to measure observables on the resulting ensembles and get reliable estimates with understood uncertainties.

Thermalization

To get an unbiased estimate from an Ensemble generated from a Markov chain we have to guarantee that the set of configurations which we account for in our expectation values do not remember the initial state the chain was started from, as that state is generated by some process that is not representative of the actual distribution of interest (in a 'cold' start, it is from the all-zero configuration, for example).

Knowing how many configurations to cut is a judgement call, and you may be misled if the observables you consider thermalize very quickly; not every observable need thermalize at once. Moreover, one can never be completely confident that your samples are drawn from the basin of lowest action (in a global sense); perhaps the Markov chain simply has not reached a preferable basin of configurations.

To facilitate measuring only on configurations after a certain step in the Markov chain, the Ensemble provides the cut method, which returns another Ensemble without the leading configurations.

Handling Autocorrelation

There are also correlations from one configuration to the next. This introduces autocorrelation into the observables time series; a naive expectation value that does not account for the autocorrelation will produce underestimated uncertainties. Good algorithms for estimating the autocorrelation time are known [5].

We can measure the autocorrelation of a timeseries.

supervillain.analysis.autocorrelation(data, mean=None, _cutoff=1e-16)[source]

The autocorrelation function is

\[\begin{align} C(\tau) &= {\left\langle \Delta(t+\tau) \Delta(t) \right\rangle} / {\left\langle \Delta(t)^2 \right\rangle} & \Delta(t) &= \texttt{data}(t) - \texttt{mean} \end{align}\]

where the ⟨averages⟩ are over the time \(t\) and \(C\) is normalized to 1 at \(\tau=0\).

The integrated autocorrelation time \(\tau_{int}\) is

\[\tau_{int} = \int_{0}^{\tau_0} d\tau\; C(\tau)\]

where \(\tau_0\) is the first time where \(C\) is zero.

Note

As defined, \(t+\tau\) does not wrap around the end of the time series, because it makes no sense to say that the very end of a Markov chain influences the generation of the beginning. Nevertheless, this implementation does include correlations as though the Markov chain were periodic, to leverage Fourier acceleration of the convolution.

Parameters

data (timeseries) – The data to correlate
mean (float) – If None, compute the mean from the data. But, if you know something about the quantity you’re considering, you might want to impose a mean value rather than compute one from the data.
_cutoff (float) – If \(C(\tau=0)\) is less than the cutoff, there is a problem (for example, no fluctuations).

Returns

C (np.array) – The autocorrelation function \(C\), the same length as the data.
\(\tau_{int}\) (int) – The ceiling of the integrated autocorrelation time.

supervillain.analysis.autocorrelation_time(data, mean=None)[source]: Just like autocorrelation() but only returns \(\tau_{int}\).

Some simple ways of decreasing autocorrelation are to decimate your Markov Chain, only keeping every nᵗʰ configuration. The Ensemble provides the every method, which returns another Ensemble ensemble keeping configurations evenly spaced by n. A natural choice for n is the autocorrelation time.

Ensembles also have an autocorrelation_time(), which leverages the above autocorrelation_time() and understands which observables to include.

Blocking

class supervillain.analysis.Blocking(ensemble, width='auto')[source]

Rather than taking every() nth configuration we can instead average (or ‘block’) the observables from consecutive configurations together.

Any observable that the underlying ensemble supports can be evaluated in the same way; you can call blocking.ObservableOfInterest to get the blocked observable of interest.

Parameters

ensemble (supervillain.Ensemble) – The ensemble to block.
width (int or ‘auto’) – The number of samples that go into each block; if ‘auto’ set by the ensemble’s autocorrelation_time().

Ensemble: The ensemble underlying the blocking.

width: The width over which to average

drop: How many configurations are dropped from the start of the ensemble to make the blocking come out evenly.

blocks: How many blocks are in the blocking.

weight: The average weight of each block.

index: The average index of each block.

index_stride: The distance between blocks.

plot_history(axes, observable, label=None, histogram_label=None, bins=31, density=True, alpha=0.5, color=None, history_kwargs={})[source]: See also

Ensemble.plot_history.

The Bootstrap

Bootstrap resampling, bootstrapping, or “the bootstrap” is a resampling method used for uncertainty estimation. One draws, with replacement, from a sample drawn according to the distribution of interest. The idea is that each draw could have been what your samples were with the same likelihood as the ensemble you actually have, and that we can estimate uncertainties by looking at distributions of means of observables from these fictitious Markov chains.

class supervillain.analysis.Bootstrap(ensemble, draws=100)[source]

The bootstrap is a resampling technique for estimating uncertainties.

For samples with weights \(w\) the expectation value of an observable is

\[\left\langle O \right\rangle = \frac{\left\langle O w \right\rangle}{\left\langle w \right\rangle}\]

and an accurate bootstrap estimate of the left-hand side requires tracking the correlations between the numerator and denominator. Moreover, quoting correlated uncertainties requires resampling different observables in the same way.

Parameters

ensemble (Ensemble) – The ensemble to resample.
draws (int) – The number of times to resample.

Any Primary Observables that Ensemble supports can be called from the Bootstrap. Bootstrap uses getattr trickery under the hood to intercept calls and perform the weighted average transparently.

Each observable returns an array of the same dimension as the ensemble’s observable. However, rather than configurations first, draws are first.

Each draw is a weighted average over the resampled weight, as shown above, and is therefore an estimator for the expectation value. These are guaranteed (by the central limit theorem) to be normally distributed as long as you have not sinned. To get an uncertainty estimate one need only take the mean() for a central value and std() for the uncertainty on the mean.

plot_band(axis, observable, color=None)[source]

Plots the single-number-valued observable as a horizontal band.

Parameters

axis (matplotlib.pyplot.axis) – The axis on which to plot.
observable (string) – Name of the observable or derived quantity.
color (matplotlib color) – See the matplotlib color API. Defaults to the previously-used color.

plot_correlator(axis, correlator, offset=0.0, irrep='A1', multiplier=1.0, linestyle='none', marker='o', markerfacecolor='none', **kwargs)[source]

Plots the space-dependent correlator against \(\Delta x\) on the axis. Plotting options and kwargs are forwarded.

Parameters

axis (matplotlib.pyplot.axis) – The axis on which to plot.
correlator (string) – Name of the observable or derived quantity.
offset (float) – Horizontal displacement, good for visually separating two correlators.
irrep (string) – Project the correlator to an irrep() understood by the lattice.
multiplier (float) – Rescales the observable by an overall constant.

estimate(observable)[source]

Parameters: observable (string) – Name of the observable or derived quantity
Returns: A tuple with the central value and uncertainty estimate for the observable. Need not be scalars, if the observable has other indices, the pieces of the tuples have those indices.
Return type: tuple

Uncertainty

class supervillain.analysis.uncertain.Uncertain(mean, uncertainty)[source]

Measurements from experiments and theoretical computations may have uncertainty. A common notation for symmetric errors, explained on Wikipedia and by NIST, is to write the uncertainty as digits in parentheses that indicate an uncertainty on the corresponding least significant digits.

So, for example, since the mass of the electron is measured to be

\[(0.51099895000±0.00000000015) \textrm{MeV}/c^2\]

we can instead write \(0.51099895000(15)\) MeV\(/c^2\).

The Uncertain class takes a mean and uncertainty and helps produce this formatted shorthand so that

electron_mass = Uncertain(0.51099895000, 0.00000000015) # MeV / c^2
print(electron_mass)
# +5.1099895000(15) × 10^-1

automatically leveraging scientific notation.

There are a few exceptional cases. First, if the uncertainty is 0, the resulting string is just the mean with no uncertainty.

print(Uncertain(3.14159,0))
# +3.14159

Second, if the uncertainty is greater than the absolute value of the mean, we resort to an explicit ± notation.

print(Uncertain(1,10))
# (+1 ± 10)

You can just specify how many digits of uncertainty to show. The default shown above is u2; u1 shows one digit in the uncertainty,

print(f'{electron_mass:u1})
# 5.109989500(2) × 10^-1

We can also format the mean with .precision as any normal float in an fstring. If the uncertainty is too small on that scale, we still show a (0) uncertainty indicator.

print(f'{electron_mass:.3}')
# 5.110(0) × 10^-1

But, you cannot specify both .precision and u[digits].

We can use scientific E notation by passing e in the format string.

print(f'{electron_mass:eu3}')
# 5.10998950000(150)e-1

To mandate a sign even for positive numbers, add a +.

print(f'{electron_mass:+eu3}')
# +5.10998950000(150)e-1

classmethod from_string(string)[source]

Parses uncertain quantities from a string.

electron_mass = Uncertain.from_string('9.1093837015(28)e-31') # kg
print(electron_mass.mean, electron_mass.uncertainty)
# 9.1093837015e-31 2.8000000000000004e-40

Comparing Results

We can get at-a-glance comparisons between ensembles. In this example we generate two ensembles from the same action and algorithm and compare their results (which ought to match, with sufficient statistics!).

(Source code, png, hires.png, pdf)

supervillain.analysis.comparison_plot.setup(observables=('ActionDensity', 'InternalEnergyDensity', 'InternalEnergyDensitySquared', 'WindingSquared'))[source]

The return values are the same as those of matplotlib.pyplot.subplots.

Parameters

observables (iterable) – The observables you wish to compare.

Returns

fig (matplotlib.pyplot.figure) – A new figure for drawing comparisons.
ax (array of axes) – Axes in the figure. One row per observable. Three columns, one for the Monte Carlo history, one for a histogram, and one for bootstraps. Even if setting up for only 1 observable, the array is two-dimensional.

supervillain.analysis.comparison_plot.bootstraps(ax, boots, labels=None, observables=('ActionDensity', 'InternalEnergyDensity', 'InternalEnergyDensitySquared', 'WindingSquared'))[source]

One row per observable, for each bootstrap object, calls plot_history() on the underlying ensemble, plot_band() on this history, and puts a bootstrap histogram in the third column.

Parameters

ax (array of axes)
boots (iterable of Bootstraps)
labels (iterable of strings)
observables (iterable of strings)

supervillain.analysis.comparison_plot.histories(ax, ensembles, labels=None, observables=('ActionDensity', 'InternalEnergyDensity', 'InternalEnergyDensitySquared', 'WindingSquared'))[source]

Calls plot_history() for each observable and row in ax. Labels the trace with the corresponding label and the autocorrelation_time() of the observable. That makes it good for ‘raw’ ensembles; ensembles that are properly decorrelated will have a very short τ.

Parameters

ax (array of axes) – As returned from setup().
ensembles (iterable of Ensembles) – Each will have its Monte Carlo history plotted and histogrammed for each observable.
labels (iterable of strings) – Names for the legend, one per ensemble.
observables (iterable of strings)