"mk" autocorrelation calculation

[grahamgower]

Hi @minaskar,

I looked at the autocorrelation stuff you've just added, and it looks really useful. Thanks!

zeus/zeus/autocorr.py

Lines 54 to 65 in 658ad94

	if method == 'mk':
	# Minas Karamanis method
	f = _autocorr_func_1d(y.reshape((-1), order='C'))
	elif method == 'dfm':
	# Daniel Forman-Mackey method
	f = np.zeros(y.shape[1])
	for yy in y:
	f += _autocorr_func_1d(yy)
	f /= len(y)
	else:
	# Goodman-Weary method
	f = _autocorr_func_1d(np.mean(y, axis=0))

I had a question/comment about the "mk" calculation. It looks like you're wanting to calculate the autocorrelation over each of the chains, which makes more sense to me than averaging over the chains (dfm) or averaging before calculating autocorrelation (Goodman-Weare). But with the mk method, by simply reshaping the 2d y array into a 1d array, the autocorrelation will be artificially lower than it should be, on account of also look at the correlation across the boundary points between chains. I.e. if the chains each have length n, then the autocorrelation of y2 = y.reshape((-1)) will be artifically lower at each position in y2 that is an integral multiple of n.

Wouldn't it be more appropriate to instead take ffts of each chain separately and average over the chains in frequency domain?

[minaskar]

Hi @grahamgower,

Thanks for your comment.

I think I agree with what you say about the "MK" recipe although I believe the effect is very minor. I haven't seen a case in which it actually affects the outcome. Please let me know if you can think of anything where this could be a problem. A nice property of this estimator, apart from its lower variance (see the figure), is that it can help diagnose violations of ergodicity. E.g. If the walkers are unable to jump from mode to mode in a multimodal case the the MK estimate of the autocorrelation time would be equal to the length of the chain.

I've tried what you suggested and it seems to give similar results to the Goodman-Weare estimator as you can see in the figure bellow:

[grahamgower]

Thanks for your insights. I guess all these estimators are asymptotically unbiased when the chain length tends to inf.

I've tried what you suggested and it seems to give similar results to the Goodman-Weare estimator as you can see in the figure bellow:

Ah! This is obvious in hindsight because the fourier transform preserves linearity. Sorry, I should admit that I haven't looked at a fourier transform in quite some years (and never used one outside a classroom).

[grahamgower]

Here's another simple alternative that seems reasonable: https://arxiv.org/abs/2009.01799. They just take the ACvF of mean-centered chains, where the mean is calculated globally over all chains.

[minaskar]

Thanks! That looks interesting, I'll have a look.

I'm not 100% sure that this applies to chains like the ones that zeus generates (from interacting walkers) but I'll have to read the paper first.