diff --git a/lectures/likelihood_bayes.md b/lectures/likelihood_bayes.md
index da81a560f..bac2fa97c 100644
--- a/lectures/likelihood_bayes.md
+++ b/lectures/likelihood_bayes.md
@@ -53,7 +53,7 @@ We'll begin by loading some Python modules.
 
 import matplotlib.pyplot as plt
 import numpy as np
-from numba import vectorize, jit, prange
+from numba import vectorize, jit
 from math import gamma
 import pandas as pd
 from scipy.integrate import quad
@@ -62,10 +62,7 @@ from scipy.integrate import quad
 import seaborn as sns
 colors = sns.color_palette()
 
-@jit
-def set_seed():
-    np.random.seed(142857)
-set_seed()
+rng = np.random.default_rng(142857)
 ```
 
 ## The Setting
@@ -163,7 +160,7 @@ g = jit(lambda x: p(x, G_a, G_b))
 
 ```{code-cell} ipython3
 @jit
-def simulate(a, b, T=50, N=500):
+def simulate(a, b, rng, T=50, N=500):
     '''
     Generate N sets of T observations of the likelihood ratio,
     return as N x T matrix.
@@ -175,7 +172,7 @@ def simulate(a, b, T=50, N=500):
     for i in range(N):
 
         for j in range(T):
-            w = np.random.beta(a, b)
+            w = rng.beta(a, b)
             l_arr[i, j] = f(w) / g(w)
 
     return l_arr
@@ -184,12 +181,12 @@ def simulate(a, b, T=50, N=500):
 We'll also use the following Python code to prepare some informative simulations
 
 ```{code-cell} ipython3
-l_arr_g = simulate(G_a, G_b, N=50000)
+l_arr_g = simulate(G_a, G_b, rng, N=50000)
 l_seq_g = np.cumprod(l_arr_g, axis=1)
 ```
 
 ```{code-cell} ipython3
-l_arr_f = simulate(F_a, F_b, N=50000)
+l_arr_f = simulate(F_a, F_b, rng, N=50000)
 l_seq_f = np.cumprod(l_arr_f, axis=1)
 ```
 
@@ -482,16 +479,16 @@ First, let's create a function to simulate data under the mixture timing protoco
 
 ```{code-cell} ipython3
 @jit
-def simulate_mixture_path(x_true, T):
+def simulate_mixture_path(x_true, T, rng):
     """
     Simulate T observations under mixture timing protocol.
     """
     w = np.empty(T)
     for t in range(T):
-        if np.random.rand() < x_true:
-            w[t] = np.random.beta(F_a, F_b)
+        if rng.random() < x_true:
+            w[t] = rng.beta(F_a, F_b)
         else:
-            w[t] = np.random.beta(G_a, G_b)
+            w[t] = rng.beta(G_a, G_b)
     return w
 ```
 
@@ -512,8 +509,8 @@ prior_params = [(1, 3), (1, 1), (3, 1)]
 prior_means = [a/(a+b) for a, b in prior_params]
 
 # Generate one path of observations from the mixture
-set_seed()
-w_mix = simulate_mixture_path(x_true, T_mix)
+rng = np.random.default_rng(142857)
+w_mix = simulate_mixture_path(x_true, T_mix, rng)
 ```
 
 ### Behavior of $\pi_t$ under wrong model
@@ -814,7 +811,7 @@ We'll plot a large sample of paths.
 
 ```{code-cell} ipython3
 @jit
-def martingale_simulate(π0, N=5000, T=200):
+def martingale_simulate(π0, rng, N=5000, T=200):
 
     π_path = np.empty((N,T+1))
     w_path = np.empty((N,T))
@@ -824,27 +821,27 @@ def martingale_simulate(π0, N=5000, T=200):
         π = π0
         for t in range(T):
             # draw w
-            if np.random.rand() <= π:
-                w = np.random.beta(F_a, F_b)
+            if rng.random() <= π:
+                w = rng.beta(F_a, F_b)
             else:
-                w = np.random.beta(G_a, G_b)
+                w = rng.beta(G_a, G_b)
             π = π*f(w)/g(w)/(π*f(w)/g(w) + 1 - π)
             π_path[n,t+1] = π
             w_path[n,t] = w
 
     return π_path, w_path
 
-def fraction_0_1(π0, N, T, decimals):
+def fraction_0_1(π0, rng, N, T, decimals):
 
-    π_path, w_path = martingale_simulate(π0, N=N, T=T)
+    π_path, w_path = martingale_simulate(π0, rng, N=N, T=T)
     values, counts = np.unique(np.round(π_path[:,-1], decimals=decimals), return_counts=True)
     return values, counts
 
-def create_table(π0s, N=10000, T=500, decimals=2):
+def create_table(π0s, rng, N=10000, T=500, decimals=2):
 
     outcomes = []
     for π0 in π0s:
-        values, counts = fraction_0_1(π0, N=N, T=T, decimals=decimals)
+        values, counts = fraction_0_1(π0, rng, N=N, T=T, decimals=decimals)
         freq = counts/N
         outcomes.append(dict(zip(values, freq)))
     table = pd.DataFrame(outcomes).sort_index(axis=1).fillna(0)
@@ -855,7 +852,7 @@ def create_table(π0s, N=10000, T=500, decimals=2):
 T = 200
 π0 = .5
 
-π_path, w_path = martingale_simulate(π0=π0, T=T, N=10000)
+π_path, w_path = martingale_simulate(π0=π0, rng=rng, T=T, N=10000)
 ```
 
 ```{code-cell} ipython3
@@ -910,7 +907,7 @@ $\pi_t$'s for various $t$'s.
 T = 200
 π0 = .3
 
-π_path3, w_path3 = martingale_simulate(π0=π0, T=T, N=10000)
+π_path3, w_path3 = martingale_simulate(π0=π0, rng=rng, T=T, N=10000)
 ```
 
 ```{code-cell} ipython3
@@ -963,7 +960,7 @@ The third column reports the fraction of $N = 10000$ simulations for which $\pi_
 
 ```{code-cell} ipython3
 # create table
-table = create_table(list(np.linspace(0,1,11)), N=10000, T=500)
+table = create_table(list(np.linspace(0,1,11)), rng, N=10000, T=500)
 table
 ```
 
@@ -989,15 +986,15 @@ Then we'll plot it.
 
 ```{code-cell} ipython3
 @jit
-def compute_cond_var(pi, mc_size=int(1e6)):
+def compute_cond_var(pi, rng, mc_size=int(1e6)):
     # create monte carlo draws
     mc_draws = np.zeros(mc_size)
 
-    for i in prange(mc_size):
-        if np.random.rand() <= pi:
-            mc_draws[i] = np.random.beta(F_a, F_b)
+    for i in range(mc_size):
+        if rng.random() <= pi:
+            mc_draws[i] = rng.beta(F_a, F_b)
         else:
-            mc_draws[i] = np.random.beta(G_a, G_b)
+            mc_draws[i] = rng.beta(G_a, G_b)
 
     dev = pi*f(mc_draws)/(pi*f(mc_draws) + (1-pi)*g(mc_draws)) - pi
     return np.mean(dev**2)
@@ -1006,7 +1003,7 @@ pi_array = np.linspace(0, 1, 40)
 cond_var_array = []
 
 for pi in pi_array:
-    cond_var_array.append(compute_cond_var(pi))
+    cond_var_array.append(compute_cond_var(pi, rng))
 
 fig, ax = plt.subplots()
 ax.plot(pi_array, cond_var_array)