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
Bisection¶
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
The implementation:
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.")
0.408293504267931
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
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