Matemáticas con código: pendientes, promedios y azar

import numpy as np
import matplotlib.pyplot as plt

def f(x):
    """La función de ejemplo: f(x) = x²."""
    return x ** 2

def pendiente_secante(f, a, h):
    """Pendiente de la recta que une (a, f(a)) con (a + h, f(a + h))."""
    return (f(a + h) - f(a)) / h

for h in [1, 0.1, 0.01, 0.001]:
    print(f"h = {h:<6} pendiente de la secante = "
          f"{pendiente_secante(f, 1, h):.4f}")

a = 1
xs = np.linspace(-0.5, 3.2, 200)
plt.plot(xs, f(xs), color="C0", linewidth=3, label="f(x) = x²")
for h, estilo in [(2, ":"), (1, "--"), (0.5, "-.")]:
    m = pendiente_secante(f, a, h)
    plt.plot(xs, f(a) + m * (xs - a), color="C2", linestyle=estilo,
             label=f"Secante, h = {h}: pendiente {m:g}")
    plt.scatter([a + h], [f(a + h)], color="C2", zorder=3)
# La tangente: la recta por (1, 1) con pendiente 2
plt.plot(xs, f(a) + 2 * (xs - a), color="C1", linewidth=2.5,
         label="Tangente: pendiente 2")
plt.scatter([a], [f(a)], color="C0", s=60, zorder=4)
plt.ylim(-1, 15)
plt.xlabel("x")
plt.ylabel("f(x)")
plt.title("Secantes que se acercan a la tangente")
plt.legend(fontsize=9, loc="upper left")
plt.show()

def derivada(f, x, h=1e-6):
    """Derivada numérica de f en x, con la diferencia centrada."""
    return (f(x + h) - f(x - h)) / (2 * h)

def cubo(x):
    return x ** 3

# La derivada exacta de x³ es 3x²: en x = 2 vale 12
for h in [0.1, 0.01, 0.001]:
    adelante = pendiente_secante(cubo, 2, h)
    centrada = derivada(cubo, 2, h)
    print(f"h = {h:<5}  hacia adelante: {adelante:.6f}   "
          f"centrada: {centrada:.6f}")
print(f"Centrada con h = 0.000001: {derivada(cubo, 2):.6f}")

def curva(x):
    """La curva de la actividad: x³/3 − x."""
    return x ** 3 / 3 - x

for x in [-2, -1, 0, 1, 2]:
    m = derivada(curva, x)
    if abs(m) < 1e-6:
        forma = "está plana"
    elif m > 0:
        forma = "sube"
    else:
        forma = "baja"
    print(f"x = {x:>2}: pendiente {round(m, 4):>5}  -> la curva {forma}")

def adentro(x):
    return 3 * x + 1

def afuera(u):
    return u ** 2

def compuesta(x):
    return afuera(adentro(x))

x0 = 1
u0 = adentro(x0)                      # u0 = 4
print("Engranaje de afuera, du² en u = 4:", round(derivada(afuera, u0), 4))
print("Engranaje de adentro, d(3x+1):    ", round(derivada(adentro, x0), 4))
print("Producto (regla de la cadena):    ",
      round(derivada(afuera, u0) * derivada(adentro, x0), 4))
print("Derivada numérica directa:        ", round(derivada(compuesta, x0), 4))

def g(x, y):
    """Un cuenco alargado: g(x, y) = x² + 3y²."""
    return x ** 2 + 3 * y ** 2

def gradiente_numerico(f, x, y, h=1e-6):
    """Las dos derivadas parciales de f en (x, y), juntas en un vector."""
    df_dx = (f(x + h, y) - f(x - h, y)) / (2 * h)
    df_dy = (f(x, y + h) - f(x, y - h)) / (2 * h)
    return np.array([df_dx, df_dy])

print("Gradiente numérico en (2, 1):",
      np.round(gradiente_numerico(g, 2, 1), 4))
print("Gradiente exacto (2x, 6y):   ", np.array([2 * 2, 6 * 1]))

punto = np.array([2.0, 1.0])
grad = gradiente_numerico(g, *punto)
paso = 1e-4                           # un paso muy corto
angulos = np.arange(360)              # una dirección por grado
rad = np.deg2rad(angulos)
direcciones = np.column_stack([np.cos(rad), np.sin(rad)])
subidas = [g(*(punto + paso * d)) - g(*punto) for d in direcciones]

print("Dirección en que g sube más:", angulos[np.argmax(subidas)],
      "grados")
angulo_grad = np.degrees(np.arctan2(grad[1], grad[0]))
print(f"Dirección del gradiente:     {angulo_grad:.1f} grados")
print(f"Mayor subida dividida por el paso: {max(subidas) / paso:.3f}")
print(f"Largo del gradiente:               {np.linalg.norm(grad):.3f}")

xs = np.linspace(-3.5, 3.5, 200)
ys = np.linspace(-2.2, 2.2, 200)
XX, YY = np.meshgrid(xs, ys)
plt.contour(XX, YY, g(XX, YY), levels=[1, 3, 7, 12, 19], colors="C2")

# Flechas del gradiente en una cuadrícula de puntos
xq, yq = np.meshgrid(np.arange(-2.5, 3, 1), [-1.2, 0, 1.2])
gx, gy = gradiente_numerico(g, xq, yq)
plt.quiver(xq, yq, gx, gy, color="C0", angles="xy",
           scale_units="xy", scale=10, width=0.007)
# El gradiente en (2, 1), sobre la curva de nivel g = 7
plt.quiver([2], [1], [grad[0]], [grad[1]], color="C1", angles="xy",
           scale_units="xy", scale=10, width=0.009)
plt.scatter([0], [0], marker="+", s=120, color="C0")
plt.gca().set_aspect("equal")
plt.xlabel("x")
plt.ylabel("y")
plt.title("Curvas de nivel de g y su gradiente")
plt.show()

# Minutos de viaje al colegio durante cinco días
sofia = np.array([20, 22, 19, 21, 18])
mateo = np.array([10, 35, 15, 30, 10])
print("Sofía:", sofia.mean(), "minutos en promedio")
print("Mateo:", mateo.mean(), "minutos en promedio")

for nombre, t in [("Sofía", sofia), ("Mateo", mateo)]:
    media = t.mean()
    varianza = np.mean((t - media) ** 2)
    desviacion = np.sqrt(varianza)
    print(f"{nombre}: media {media:.1f}   varianza {varianza:.1f}   "
          f"desviación {desviacion:.2f}")

print("Con NumPy:", np.var(mateo), round(np.std(mateo), 2))

z_sofia = (23 - sofia.mean()) / sofia.std()
z_mateo = (40 - mateo.mean()) / mateo.std()
print(f"23 minutos para Sofía: z = {z_sofia:.2f}")
print(f"40 minutos para Mateo: z = {z_mateo:.2f}")

z = (mateo - mateo.mean()) / mateo.std()
print("Mateo estandarizado:", np.round(z, 2))
print(f"Media: {abs(z.mean()):.2f}   desviación: {z.std():.2f}")

# Cuatro estudiantes: estatura (cm) y peso (kg)
tabla = np.array([[150, 45], [160, 52], [170, 61], [180, 70]],
                 dtype=float)
tabla_est = (tabla - tabla.mean(axis=0)) / tabla.std(axis=0)

print("Medias por columna:", tabla.mean(axis=0))
print("Desviaciones:", np.round(tabla.std(axis=0), 2))
print(np.round(tabla_est, 2))

rng = np.random.default_rng(7)
tiros = rng.integers(1, 7, size=100_000)
print("Primeros diez tiros:", tiros[:10])

for n in [10, 100, 1_000, 10_000, 100_000]:
    fraccion = np.mean(tiros[:n] == 6)
    print(f"{n:>7} tiros: fracción de seises = {fraccion:.4f}")
print(f"Probabilidad: 1/6 = {1 / 6:.4f}")

n_tiros = np.arange(1, len(tiros) + 1)
# Fracción de seises acumulada después de cada tiro
frecuencia = np.cumsum(tiros == 6) / n_tiros
plt.plot(n_tiros, frecuencia, color="C0", label="Fracción de seises")
plt.axhline(1 / 6, color="C2", linestyle="--", label="1/6")
plt.xscale("log")
plt.ylim(0, 0.5)
plt.xlabel("Cantidad de tiros (escala logarítmica)")
plt.ylabel("Fracción de seises")
plt.title("La frecuencia se acerca a 1/6")
plt.legend()
plt.show()

N = 1_000_000                         # personas simuladas
p_enfermedad = 0.01                   # 1 de cada 100
sensibilidad = 0.90                   # P(positivo | enferma)
falso_positivo = 0.05                 # P(positivo | sana)

rng = np.random.default_rng(1)
enferma = rng.random(N) < p_enfermedad
azar = rng.random(N)
# A cada persona le sale positivo con la probabilidad que le toca
positivo = np.where(enferma, azar < sensibilidad, azar < falso_positivo)

print("Personas con resultado positivo:", positivo.sum())
print(f"P(enferma | positivo) simulada: {enferma[positivo].mean():.4f}")
p_positivo = (sensibilidad * p_enfermedad
              + falso_positivo * (1 - p_enfermedad))
print(f"P(enferma | positivo) con Bayes: "
      f"{sensibilidad * p_enfermedad / p_positivo:.4f}")

rng = np.random.default_rng(5)
estaturas = rng.normal(165, 7, size=100_000)
print(f"Media: {estaturas.mean():.2f}   "
      f"desviación: {estaturas.std():.2f}")

z = (estaturas - estaturas.mean()) / estaturas.std()
for k in [1, 2, 3]:
    print(f"Fracción con |z| menor que {k}: {np.mean(np.abs(z) < k):.3f}")

def densidad_normal(x, mu, sigma):
    """La fórmula de la campana."""
    return (np.exp(-(x - mu) ** 2 / (2 * sigma ** 2))
            / (sigma * np.sqrt(2 * np.pi)))

rng = np.random.default_rng(3)
fig, ejes = plt.subplots(2, 2, figsize=(4.8, 4.8))
for eje, k in zip(ejes.flat, [1, 2, 5, 20]):
    # 100 000 tiradas de k dados, sumando cada tirada
    sumas = rng.integers(1, 7, size=(100_000, k)).sum(axis=1)
    valores, cuentas = np.unique(sumas, return_counts=True)
    eje.bar(valores, cuentas / len(sumas), width=0.8, color="C0")
    eje.set_title("1 dado" if k == 1 else f"Suma de {k} dados")
    eje.set_xticks([k, 6 * k])          # la menor y la mayor suma
    eje.set_yticks([])

# Un dado tiene media 3,5 y varianza 35/12; con k dados, k veces eso
xs = np.linspace(20, 120, 300)
eje.plot(xs, densidad_normal(xs, 3.5 * k, np.sqrt(35 / 12 * k)),
         color="C1", linewidth=2)
fig.suptitle("Sumar dados forma una campana")
plt.show()

print("e =", np.e)
print("exp(1) =", np.exp(1), "  log(e) =", np.log(np.e))
print("log(exp(5)) =", np.log(np.exp(5)))
print("exp(10) =", round(np.exp(10), 1), "  log(10) =",
      round(np.log(10), 4))

xs = np.linspace(-3, 3, 200)
plt.plot(xs, np.exp(xs), color="C0", label="exp(x)")
xs_pos = np.linspace(0.02, 9, 300)     # log solo acepta positivos
plt.plot(xs_pos, np.log(xs_pos), color="C0", linestyle="--",
         label="log(x)")
plt.plot([-3, 9], [-3, 9], color="C2", linestyle=":",
         label="y = x")
plt.xlim(-3, 9)
plt.ylim(-3, 9)
plt.gca().set_aspect("equal")
plt.xlabel("x")
plt.title("exp y log son inversas")
plt.legend(loc="center left", bbox_to_anchor=(1, 0.5))
plt.show()

print("log(6 · 8)      =", round(np.log(6 * 8), 6))
print("log(6) + log(8) =", round(np.log(6) + np.log(8), 6))

# La probabilidad de 1000 sucesos independientes, cada uno de 1 %
probabilidades = np.full(1000, 0.01)
print("Producto directo:", np.prod(probabilidades))
print("Suma de logaritmos:", round(np.sum(np.log(probabilidades)), 2))

for x in [0, 1, 2]:
    print(f"x = {x}: exp(x) = {np.exp(x):.4f}   "
          f"pendiente = {derivada(np.exp, x):.4f}")