import jax.numpy as jnp

from jax import random

np.random.seed(142857)

@jit

def set_seed():

np.random.seed(142857)

set_seed()

rng = np.random.default_rng(142857)

```

Let's use Python to generate two beta distributions

:hide-output: false

@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.

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

```{code-cell} ipython3

:hide-output: false

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

:hide-output: false

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)

```

* Step one:

* use the numpy.random.choice function to flip an unfair coin that selects distribution $F$ with prob $\alpha$

* use the `rng.choice` method to flip an unfair coin that selects distribution $F$ with prob $\alpha$

and $G$ with prob $1 -\alpha$

* Step two:

```{code-cell} ipython3

@jit

def draw_lottery(p, N):

def draw_lottery(p, rng, N):

"Draw from the compound lottery directly."

draws = []

for i in range(0, N):

if np.random.rand()<=p:

draws.append(np.random.beta(F_a, F_b))

if rng.random()<=p:

draws.append(rng.beta(F_a, F_b))

else:

draws.append(np.random.beta(G_a, G_b))

draws.append(rng.beta(G_a, G_b))

return np.array(draws)

def draw_lottery_MC(p, N):

def draw_lottery_MC(p, rng, N):

"Draw from the compound lottery using the Monte Carlo trick."

xs = np.linspace(1e-8,1-(1e-8),10000)

CDF = p*sp.beta.cdf(xs, F_a, F_b) + (1-p)*sp.beta.cdf(xs, G_a, G_b)

Us = np.random.rand(N)

Us = rng.random(N)

draws = xs[np.searchsorted(CDF[:-1], Us)]

return draws

```

N = 100000

α = 0.0

sample1 = draw_lottery(α, N)

sample2 = draw_lottery_MC(α, N)

sample1 = draw_lottery(α, rng, N)

sample2 = draw_lottery_MC(α, rng, N)

# plot draws and density function

plt.hist(sample1, 50, density=True, alpha=0.5, label='direct draws')

Let's set a value of $\alpha$ and then watch how $\pi_t$ evolves.

```{code-cell} ipython3

def simulate_mixed(α, T=50, N=500):

def simulate_mixed(α, rng, T=50, N=500):

"""

Generate N sets of T observations of the likelihood ratio,

return as N x T matrix, when the true density is mixed h;α

"""

w_s = draw_lottery(α, N*T).reshape(N, T)

w_s = draw_lottery(α, rng, N*T).reshape(N, T)

l_arr = f(w_s) / g(w_s)

return l_arr

def plot_π_seq(α, π1=0.2, π2=0.8, T=200):

def plot_π_seq(α, rng, π1=0.2, π2=0.8, T=200):

"""

Compute and plot π_seq and the log likelihood ratio process

when the mixed distribution governs the data.

"""

l_arr_mixed = simulate_mixed(α, T=T, N=50)

l_arr_mixed = simulate_mixed(α, rng, T=T, N=50)

l_seq_mixed = np.cumprod(l_arr_mixed, axis=1)

T = l_arr_mixed.shape[1]

```

```{code-cell} ipython3

plot_π_seq(α = 0.6)

plot_π_seq(α = 0.6, rng=rng)

```

The above graph shows a sample path of the log likelihood ratio process as the blue dotted line, together with

Let's see what happens when we change $\alpha$

```{code-cell} ipython3

plot_π_seq(α = 0.2)

plot_π_seq(α = 0.2, rng=rng)

```

Evidently, $\alpha$ is having a big effect on the destination of $\pi_t$ as $t \rightarrow + \infty$

"Find limit of π sequence."

π_seq = np.zeros(T+1)

π_seq[0] = π_0

l_arr = simulate_mixed(α, T, N=1)[0]

l_arr = simulate_mixed(α, rng, T, N=1)[0]

for t in range(T):

π_seq[t+1] = update(π_seq[t], l_arr[t])

α = 0.8

# simulate data with true α

data = draw_lottery(α, 1000)

data = draw_lottery(α, rng, 1000)

sizes = [5, 20, 50, 200, 1000, 25000]

def model(w):

prior_params = [(1, 3), (1, 1), (3, 1)]

prior_means = [a/(a+b) for a, b in prior_params]

w_mix = draw_lottery(x_true, T_mix)

w_mix = draw_lottery(x_true, rng, T_mix)

```

```{code-cell} ipython3

Next, let's look at multiple simulations with a longer time horizon, all starting from a uniform prior.

```{code-cell} ipython3

set_seed()

rng = np.random.default_rng(142857)

n_paths = 20

T_long = 10_000

fig, ax = plt.subplots(figsize=(10, 5))

for j in range(n_paths):

w_path = draw_lottery(x_true, T_long)

w_path = draw_lottery(x_true, rng, T_long)

x_means = learn_x_bayesian(w_path, 1, 1) # Uniform prior

ax.plot(range(T_long + 1), x_means, alpha=0.5, linewidth=1)