gistfile1.py

import sys

from matplotlib.pyplot import *
from numpy import *
from numpy.random import *



def autocorr(X, N, d=1):
	"""
	Estimates autocorrelation from a sample of a possibly multivariate
	stationary Markov chain.
	"""

	X = X - mean(X, 1).reshape(-1, 1)
	v = mean(sum(square(X), 0))

	# autocovariance
	A = [v]

	for t in range(1, N + 1, d):
		A.append(mean(sum(X[:, :-t] * X[:, t:], 0)))

	# normalize by variance
	return hstack(A) / v



def main(argv):
	ion()

	# hyperparameters
	dim = 4
	num_samples = 10000
	sigmas = arange(0.2, 3.0, 0.05)

	for sigma in sigmas:
		### JOINT ACCEPT/REJECT

		X = randn(dim, num_samples)

		for i in range(1, num_samples):
			# propose step
			X[:, i] = X[:, i - 1] + sigma * X[:, i]

			# accept/reject step
			if rand() > exp(0.5 * (sum(square(X[:, i - 1]), 0) - sum(square(X[:, i]), 0))):
				X[:, i] = X[:, i - 1]

		### INDEPENDENT ACCEPT/REJECT

		Y = randn(dim, num_samples)

		for i in range(1, num_samples):
			# propose steps
			Y[:, i] = Y[:, i - 1] + sigma * Y[:, i]

			# accept/reject steps
			for j in range(dim):
				if rand() > exp(0.5 * (square(Y[j, i - 1]) - square(Y[j, i]))):
					Y[j, i] = Y[j, i - 1]

		plot(autocorr(X, 100), 'r')
		plot(autocorr(Y, 100), 'b')
		draw()

	xlabel('Time')
	ylabel('Correlation')
	title('Autocorrelation')
	legend(['Jointly', 'Independently'], loc='upper right')

	raw_input()

	return 0



if __name__ == '__main__':
	sys.exit(main(sys.argv))