23.b)

import random as rnd
import matplotlib.pyplot as plt
import math


def p(x, σ, μ):
    return 1/(math.sqrt(2*math.pi)*σ) * math.exp(-1/(2*σ**2) * (x-μ)**2)


def aceptar(a, b, σ, μ):
    xt = 0
    while 1:
        xt = rnd.random()*(b - a) + a
        r = rnd.random()/(math.sqrt(2*math.pi)*σ)
        if r < p(xt, σ, μ):
            break
    return xt


n = 100_000
b = 8
a = -8
σ = math.sqrt(1.5)
μ = 2
x = []

for i in range(n):
    x.append(aceptar(a, b, σ, μ))

plt.subplot(2, 1, 1)
plt.hist(x, 100)
plt.ylabel("yo")

y = []
for i in range(n):
    y.append(rnd.normalvariate(μ, σ))

plt.subplot(2, 1, 2)
plt.hist(y, 100)
plt.ylabel("random::normalvariate")
plt.show()