Chapter 2 Code¶
We’ll use Python’s NumPy library in what follows, as well as Matplotlib for plotting.
import numpy as np import matplotlib.pyplot as plt
Here’s an implementation of the bisection algorithm. We will test it on this function:
f = lambda x: np.sin(4 * (x - 1/4)) + x + x**20 - 1
M = 1000 ϵ = 1e-8 α, β = 0, 1 i = 1 a, b = α, β while i <= M: c = (a + b) / 2 if abs(f(c)) < ϵ: print(c) break i += 1 if f(a) * f(c) < 0: b = c else: a = c if i > M: print("Failed to converge.")
User Defined Functions¶
Here’s the basic implementation of the function \(\tau\).
def tau(z, S, phi): """ Evaluates the function tau(z) given data S, phi, where S and phi are assumed to be arrays. """ a = 0 for i, x in enumerate(S): b = a + phi[i] if a < z <= b: return x a = b
Here’s a more efficient implementation.
def tau(z, S, phi): i = np.searchsorted(np.cumsum(phi), z) return S[i]
And here’s a closure that generates the function \(\tau\).
def tau_factory(S, phi): Φ = np.cumsum(phi) def tau(z): i = np.searchsorted(Φ, z) return S[i] return tau
We generate a function \(\tau\) that acts on \(z\) alone by calling the function factory:
phi = 0.2, 0.5, 0.3 S = 0, 1, 2 tau = tau_factory(S, phi)
tau(0.1) # Should be 0
All of these functions work as expected. To illustrate, here \(\tau\) is used to generate draws from a given distribution \(\phi\).
size = 100_000 draws = np.empty(size) for j in range(size): W = np.random.uniform() draws[j] = tau(W) # Compute fraction of draws with each possible value frequency = [np.mean(draws==j) for j in S]
Let’s check that the empirical frequency approximately coincides with the probabilities in \(\phi\).
fig, ax = plt.subplots() ax.bar(S, frequency, edgecolor='k', facecolor='b', alpha=0.25, label="histogram") ax.stem(S, phi, label='$\\phi$') ax.legend() plt.show()
Object Oriented Programming¶
Here’s a class that implements the function \(\tau\) as a method, as well as a method to generate draws from \(\phi\).
class Tau: def __init__(self, S, phi): self.S = S self.Φ = np.cumsum(phi) def tau(self, z): i = np.searchsorted(self.Φ, z) return self.S[i] def draw(self): W = np.random.uniform() return self.tau(W)
tau = Tau(S, phi)
for i in range(5): print(tau.draw())
0 1 1 1 1