@@ -132,8 +132,9 @@ Here is one solution:
132132``` {code-cell} ipython3
133133class Frequentist:
134134
135- def __init__(self, θ, n, I):
135+ def __init__(self, θ, n, I, rng=None ):
136136 self.θ, self.n, self.I = θ, n, I
137+ self.rng = rng or np.random.default_rng()
137138
138139 def binomial(self, k):
139140 '''Compute the theoretical probability.'''
@@ -143,7 +144,7 @@ class Frequentist:
143144 def draw(self):
144145 '''Draw n independent flips for I sequences.'''
145146 θ, n, I = self.θ, self.n, self.I
146- sample = np.random.rand( I, n)
147+ sample = self.rng.random(( I, n) )
147148 self.Y = (sample <= θ).astype(int)
148149
149150 def compute_fk(self, k):
@@ -165,9 +166,10 @@ class Frequentist:
165166```
166167
167168``` {code-cell} ipython3
169+ rng = np.random.default_rng(123)
168170θ, n, k, I = 0.7, 20, 10, 1_000_000
169171
170- freq = Frequentist(θ, n, I)
172+ freq = Frequentist(θ, n, I, rng=rng )
171173
172174freq.compare()
173175```
190192We'll vary $\theta$ from $0.01$ to $0.99$ and plot outcomes against $\theta$.
191193
192194``` {code-cell} ipython3
195+ rng = np.random.default_rng(234)
193196θ_low, θ_high, n_thetas = 0.01, 0.99, 50
194197thetas = np.linspace(θ_low, θ_high, n_thetas)
195198P = []
196199f_kI = []
197200for i in range(n_thetas):
198- freq = Frequentist(thetas[i], n, I)
201+ freq = Frequentist(thetas[i], n, I, rng=rng )
199202 freq.binomial(k)
200203 freq.draw()
201204 freq.compute_fk(k)
@@ -223,12 +226,13 @@ Now we fix $\theta=0.7, k=10, I=1,000,000$ and vary $n$ from $1$ to $100$.
223226Then we'll plot outcomes.
224227
225228``` {code-cell} ipython3
229+ rng = np.random.default_rng(345)
226230n_low, n_high, n_ns = 1, 100, 50
227231ns = np.linspace(n_low, n_high, n_ns, dtype='int')
228232P = []
229233f_kI = []
230234for i in range(n_ns):
231- freq = Frequentist(θ, ns[i], I)
235+ freq = Frequentist(θ, ns[i], I, rng=rng )
232236 freq.binomial(k)
233237 freq.draw()
234238 freq.compute_fk(k)
@@ -254,13 +258,14 @@ plt.show()
254258Now we fix $\theta=0.7, n=20, k=10$ and vary $\log(I)$ from $2$ to $6$.
255259
256260``` {code-cell} ipython3
261+ rng = np.random.default_rng(456)
257262I_log_low, I_log_high, n_Is = 2, 6, 200
258263log_Is = np.linspace(I_log_low, I_log_high, n_Is)
259264Is = np.power(10, log_Is).astype(int)
260265P = []
261266f_kI = []
262267for i in range(n_Is):
263- freq = Frequentist(θ, n, Is[i])
268+ freq = Frequentist(θ, n, Is[i], rng=rng )
264269 freq.binomial(k)
265270 freq.draw()
266271 freq.compute_fk(k)
406411``` {code-cell} ipython3
407412class Bayesian:
408413
409- def __init__(self, θ=0.4, n=1_000_000, α=0.5, β=0.5):
414+ def __init__(self, θ=0.4, n=1_000_000, α=0.5, β=0.5,
415+ rng=None):
410416 '''
411417 Parameters
412418 ----------
413419 θ : Probability of heads on each flip.
414420 n : Number of flips in the sequence.
415421 α, β : Parameters of the beta prior on θ.
422+ rng : NumPy random generator.
416423 '''
417424 self.θ, self.n, self.α, self.β = θ, n, α, β
425+ self.rng = rng or np.random.default_rng()
418426 self.prior = st.beta(α, β)
419427
420428 def draw(self):
421429 '''Simulate a sequence of n coin flips.'''
422- array = np.random.rand (self.n)
430+ array = self.rng.random (self.n)
423431 self.draws = (array < self.θ).astype(int)
424432
425433 def form_single_posterior(self, n_obs):
@@ -440,7 +448,8 @@ class Bayesian:
440448** d)** Please plot the posterior distribution for $\theta$ as a function of $\theta$ as $n$ grows from $1, 2, \ldots$.
441449
442450``` {code-cell} ipython3
443- bayes = Bayesian()
451+ rng = np.random.default_rng(567)
452+ bayes = Bayesian(rng=rng)
444453bayes.draw()
445454
446455n_obs_list = [1, 2, 3, 4, 5, 10, 20, 30, 50, 70,
@@ -523,7 +532,7 @@ ax.set_xlabel('Number of Observations', fontsize=11)
523532plt.show()
524533```
525534
526- Notice that in the graph above the posterior probability that $\theta \in [ .45, .55] $ typically exhibits a hump shape as $n$ increases.
535+ Notice that in the graph above the posterior probability that $\theta \in [ .45, .55] $ exhibits a hump shape as $n$ increases.
527536
528537Two opposing forces are at work.
529538
0 commit comments