markov.py

# RA, 2020-01-01


from numpy import sqrt, exp
import numpy as np
from scipy.special import binom, erf
import matplotlib.pyplot as plt

# Values of b
bb = [0.1, 0.5, 1.5]

fig: plt.Figure
ax: plt.Axes
(fig, AX) = plt.subplots(1, len(bb), figsize=(5 * len(bb), 5))

for (b, ax) in zip(bb, AX):

	phi = (lambda t: exp(b * np.abs(t)))

	# Quadrature nodes count
	nn = np.arange(2, 21)

	# Exact values, i.e. E[Φ[Z]]
	I = exp((b ** 2) / 2) * (1 + erf(b / sqrt(2))) * np.ones(nn.shape)

	# Quadrature values, i.e. E[Φ[X]]
	Q = np.zeros(nn.shape)

	for n in nn:
		# Quadrature nodes
		x = sqrt(n) * np.linspace(-1, 1, n + 1)
		# Weights
		w = np.asarray([binom(n, m) / (2 ** n) for m in range(0, n + 1)])
		assert (abs(1 - sum(w)) <= 1e-8)

		# E[Φ[X]]
		Q[n == nn] = np.dot(w, phi(x))

	ax.plot(nn, I, 'b-', label="E[Φ[Z]]")
	ax.plot(nn, Q, 'rs', label="E[Φ[X]]")

	ax.set_xticks(nn)
	ax.set_xlabel("n")
	ax.set_title(F"b = {b}")
	ax.legend()
	ax.grid()

fig.savefig("markov.png", dpi=180, bbox_inches='tight', pad_inches=0)
plt.show()