markov.py

# RA, 2020-01-02


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

fig: plt.Figure
ax: plt.Axes
(fig, ax) = plt.subplots()

# Quadrature nodes count
nn = np.arange(1, 41)

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

	# Values of b
	bb = np.unique(np.hstack([np.linspace(0.1, 10, 101), np.sqrt(np.arange(1, 11))]))

	# Exact values, i.e. E[Φ(Z)]
	EZ = (2 / bb)

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

	# Phi values, i.e. Phi(b)
	Phi = np.zeros(bb.shape)

	for b in bb:
		# phi = (lambda t: exp(b * np.abs(t)))
		phi = (lambda t: exp(-b * (np.abs(t) - b)) * (np.abs(t) >= b) / (1 / sqrt(2 * pi) * exp(-1 / 2 * (t ** 2))))

		# E[Φ(X)]
		EX[b == bb] = np.dot(w, phi(x))

		# Φ(b)
		Phi[b == bb] = phi(b)

	# Nasty factor
	C = EX / EZ

	# Exact probability P[X >= b]
	P = np.asarray([sum(w[abs(x) >= b]) for b in bb])

	# Markov estimate
	Q = EX / Phi

	ax.plot(bb, C, '.k')
	ax.set_title("C for different n and b")

	ax.set_xlabel("b")
	ax.grid()

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