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

Las exactas: 2 · 3 = 6; −1 · 2⁻² = −0,25; 1/0,5 = 2.

def inversa(x):
    return 1 / x

casos = [("x²", f, 3, 6), ("1/x", inversa, 2, -0.25),
         ("log(x)", np.log, 0.5, 2)]
print(f"{'':>13}{'exacta':>11}{'adelante':>11}{'centrada':>11}"
      f"{'fina':>11}")
for nombre, fun, x, exacta in casos:
    adelante = pendiente_secante(fun, x, 0.01)
    centrada = derivada(fun, x, 0.01)
    fina = derivada(fun, x)            # h = 0.000001
    print(f"{nombre + ' en ' + str(x):>13}{exacta:>11.6f}"
          f"{adelante:>11.6f}{centrada:>11.6f}{fina:>11.6f}")

                  exacta   adelante   centrada       fina
      x² en 3   6.000000   6.010000   6.000000   6.000000
     1/x en 2  -0.250000  -0.248756  -0.250006  -0.250000
log(x) en 0.5   2.000000   1.980263   2.000267   2.000000

Con h = 0,01, la secante hacia adelante se equivoca ya en la segunda o tercera cifra decimal; la centrada, recién en la cuarta o más allá. Con x² la centrada da el valor exacto: con una parábola, la diferencia centrada nunca se equivoca. Con h = 0,000001, la centrada coincide con la exacta en las seis cifras que mostramos. La de log(x) en 0,5 es la que más error tiene, porque cerca de 0 el logaritmo se curva mucho.

for x in [-1.5, 0.5, 1.5]:
    print(f"x = {x:>4}: pendiente {derivada(curva, x):+.4f}")

def altura(t):
    return 80 - 5 * t ** 2

print(f"Derivada de la altura en t = 2: {derivada(altura, 2):.4f}")

x = -1.5: pendiente +1.2500
x =  0.5: pendiente -0.7500
x =  1.5: pendiente +1.2500
Derivada de la altura en t = 2: -20.0000

a) En −1,5 y en 1,5 la pendiente es +1,25: la curva sube, y por cada unidad que avanzas a la derecha sube unas 1,25 unidades. En 0,5 la pendiente es −0,75: la curva baja, y menos empinada. Coincide con la derivada exacta x² − 1.

b) La derivada es −20. El signo negativo dice que la altura disminuye: la piedra cae. El valor dice qué tan rápido: en ese instante pierde 20 metros por cada segundo, es decir, cae a 20 m/s. La derivada de la posición respecto del tiempo es la velocidad.

a) La derivada de afuera es exp(u) y la de adentro es −2x. Multiplicando: la derivada es −2x · exp(−x²).

b) En x = 1: −2 · exp(−1) ≈ −0,7358.

def campana(x):
    return np.exp(-x ** 2)

print(f"Con la regla de la cadena: {-2 * 1 * np.exp(-1):.4f}")
print(f"Derivada numérica:         {derivada(campana, 1):.4f}")

Con la regla de la cadena: -0.7358
Derivada numérica:         -0.7358

c) −2x · exp(−x²) vale cero solo cuando x = 0, porque la exponencial nunca vale cero. Ahí la campana está plana: es su punto más alto, un máximo.

a) Para ∂h/∂x, y es un número fijo: la derivada de x² es 2x, la de x·y es y (como la de 3x es 3) y 2y² no cambia, así que su derivada es 0. Queda ∂h/∂x = 2x + y. Del mismo modo, ∂h/∂y = x + 4y. En (1, 1), el gradiente es (3, 5).

def h_xy(x, y):
    return x ** 2 + x * y + 2 * y ** 2

p = np.array([1.0, 1.0])
grad_h = gradiente_numerico(h_xy, *p)
print("Gradiente numérico en (1, 1):", np.round(grad_h, 4))

# Vectores de largo 1: a favor, en contra y perpendicular
unitario = grad_h / np.linalg.norm(grad_h)
perpendicular = np.array([-unitario[1], unitario[0]])
for nombre, d in [("a favor", unitario), ("en contra", -unitario),
                  ("perpendicular", perpendicular)]:
    cambio = h_xy(*(p + 0.01 * d)) - h_xy(*p)
    print(f"Paso {nombre:>13}: h cambia {cambio:+.5f}")
print(f"Largo del gradiente × 0.01:  {np.linalg.norm(grad_h) * 0.01:.5f}")

Gradiente numérico en (1, 1): [3. 5.]
Paso       a favor: h cambia +0.05853
Paso     en contra: h cambia -0.05809
Paso perpendicular: h cambia +0.00008
Largo del gradiente × 0.01:  0.05831

c) A favor del gradiente, h sube casi exactamente «largo del gradiente × largo del paso»: es la mayor subida posible con ese paso. En contra, baja casi lo mismo: la mayor bajada posible. En la dirección perpendicular el cambio es casi cero: caminas a lo largo de la curva de nivel. (No es cero exacto porque las curvas de nivel están curvadas; con un paso más chico, el cambio perpendicular se achica mucho más rápido que los otros dos.)

a) La media es 40/8 = 5. Las distancias a la media son −3, −1, −1, −1, 0, 0, 2 y 4; sus cuadrados suman 9 + 1 + 1 + 1 + 0 + 0 + 4 + 16 = 32. La varianza es 32/8 = 4 y la desviación estándar, √4 = 2.

b) Cada z es la distancia a la media dividida por 2: −1,5; −0,5; −0,5; −0,5; 0; 0; 1; 2. Quien durmió 2 horas está a 1,5 desviaciones por debajo del promedio.

import pandas as pd

datos = np.array([2, 4, 4, 4, 5, 5, 7, 9])
print("Media:", datos.mean(), "  varianza:", datos.var(),
      "  desviación:", datos.std())
print("z:", (datos - datos.mean()) / datos.std())
print(f"pandas, por defecto: {pd.Series(datos).std():.4f}")
print(f"pandas con ddof=0:   {pd.Series(datos).std(ddof=0):.4f}")

Media: 5.0   varianza: 4.0   desviación: 2.0
z: [-1.5 -0.5 -0.5 -0.5  0.   0.   1.   2. ]
pandas, por defecto: 2.1381
pandas con ddof=0:   2.0000

c) NumPy da lo mismo que la cuenta a mano. pandas, por defecto, divide por n − 1 = 7 en lugar de por 8: √(32/7) ≈ 2,138. Con ddof=0 divide por n y coincide con NumPy. Ninguna de las dos está «mal», pero hay que usar la misma en todo un trabajo.

a) De los 100 enfermos, sacan dos positivos 100 · 0,9 · 0,9 = 81. De los 9900 sanos, 9900 · 0,05 · 0,05 = 24,75 (en promedio). La probabilidad es 81/(81 + 24,75) ≈ 0,766: con dos positivos, la sospecha sube del 15 % a casi el 77 %. Es la regla de Bayes aplicada dos veces: lo que salió del primer positivo (15,4 %) es el punto de partida del segundo.

N = 1_000_000
rng = np.random.default_rng(2)
enferma = rng.random(N) < 0.01
# Cada prueba se decide con su propio número al azar
prob_pos = np.where(enferma, 0.90, 0.05)
prueba1 = rng.random(N) < prob_pos
prueba2 = rng.random(N) < prob_pos
dos_positivos = prueba1 & prueba2

print("Personas con dos positivos:", dos_positivos.sum())
print(f"Simulación: {enferma[dos_positivos].mean():.3f}")
formula = (0.01 * 0.9 ** 2) / (0.01 * 0.9 ** 2 + 0.99 * 0.05 ** 2)
print(f"Fórmula:    {formula:.3f}")

Personas con dos positivos: 10422
Simulación: 0.768
Fórmula:    0.766

La simulación (0,768) y la fórmula (0,766) coinciden en dos cifras: solo unas 10 000 personas sacaron dos positivos, así que la simulación tiene un margen de unas milésimas. Fíjate que prob_pos da a cada persona la probabilidad de positivo que le toca según esté enferma o sana; comparar rng.random(N) < prob_pos hace una Bernoulli distinta para cada una.

a) 179 está a (179 − 165)/7 = 2 desviaciones por encima de la media. Entre −2 y +2 desviaciones queda el 95 %; el 5 % restante se reparte por igual entre las dos colas, así que más de 179 cm mide alrededor del 2,5 % (el valor exacto es 2,3 %). Entre 158 y 172 cm es estar a menos de una desviación de la media: alrededor del 68 %.

rng = np.random.default_rng(8)
alturas = rng.normal(165, 7, size=100_000)
print(f"Más de 179 cm:       {np.mean(alturas > 179):.3f}")
print(f"Entre 158 y 172 cm:  {np.mean((alturas > 158) & (alturas < 172)):.3f}")

Más de 179 cm:       0.023
Entre 158 y 172 cm:  0.682

Las fracciones simuladas quedan muy cerca de 0,023 y de 0,683, los valores exactos de la normal.

a) Si a = e^u y b = e^v, entonces a/b = e^{u − v}, así que log(a/b) = u − v = log(a) − log(b). Igual que el producto se vuelve suma, la división se vuelve resta.

print("log(50 / 4)        =", round(np.log(50 / 4), 6))
print("log(50) - log(4)   =", round(np.log(50) - np.log(4), 6))

directa = (1 / 6) ** 500
L = 500 * np.log(1 / 6)               # log de un producto = suma de logs
print("Directa:", directa)
print(f"Con logaritmos: log(p) = {L:.2f}")
print(f"En potencia de 10: p = 10 elevado a {L / np.log(10):.2f}")

log(50 / 4)        = 2.525729
log(50) - log(4)   = 2.525729
Directa: 0.0
Con logaritmos: log(p) = -895.88
En potencia de 10: p = 10 elevado a -389.08

b) La cuenta directa da 0: el número es tan chico que la computadora no lo puede guardar. Con logaritmos no hay problema: log(p) ≈ −895,88. Es lo que hace un modelo cuando combina muchas probabilidades.

c) p ≈ 10⁻³⁸⁹: escrito con decimales, es un 0, una coma y 389 ceros antes de la primera cifra distinta de cero. La computadora solo guarda números de hasta unos 10⁻³⁰⁸.

caras = np.arange(1, 7)
justo = np.full(6, 1 / 6)
cargado = np.array([0.15, 0.15, 0.15, 0.15, 0.15, 0.25])

def tirar(n, probs, rng):
    """Tira n veces un dado con esas probabilidades."""
    return rng.choice(caras, size=n, p=probs)

def fracciones(tiros):
    """Fracción de veces que salió cada cara, de 1 a 6."""
    return np.bincount(tiros, minlength=7)[1:] / len(tiros)

rng = np.random.default_rng(10)

# Paso 2: pocas tiradas engañan
for n in [60, 6000]:
    for nombre, probs in [("justo", justo), ("cargado", cargado)]:
        print(f"{n:>5} tiros, {nombre:>7}:",
              np.round(fracciones(tirar(n, probs, rng)), 3))

# Paso 3: la fracción de seises de un dado justo, 10 000 veces
def fraccion_seises(n, probs, rng, repeticiones=10_000):
    tiros = rng.choice(caras, size=(repeticiones, n), p=probs)
    return np.mean(tiros == 6, axis=1)

seises_justo = fraccion_seises(100, justo, rng)
bajo, alto = np.percentile(seises_justo, [2.5, 97.5])
print(f"Con 100 tiros, el 95 % de los dados justos da entre "
      f"{bajo:.2f} y {alto:.2f}")

# Paso 4: cuántas veces la regla acusa a un dado justo
print(f"Falsa alarma: {np.mean(seises_justo > alto):.3f}")

# Paso 5: cuántos tiros hacen falta
for n in [30, 100, 300, 1000]:
    limite = np.percentile(fraccion_seises(n, justo, rng), 97.5)
    descubre = np.mean(fraccion_seises(n, cargado, rng) > limite)
    print(f"n = {n:>4}: límite {limite:.3f}, "
          f"descubre al cargado el {descubre:.0%} de las veces")

   60 tiros,   justo: [0.15  0.15  0.167 0.15  0.2   0.183]
   60 tiros, cargado: [0.133 0.233 0.183 0.167 0.1   0.183]
 6000 tiros,   justo: [0.158 0.169 0.174 0.164 0.166 0.168]
 6000 tiros, cargado: [0.151 0.147 0.152 0.149 0.149 0.252]
Con 100 tiros, el 95 % de los dados justos da entre 0.10 y 0.24
Falsa alarma: 0.021
n =   30: límite 0.300, descubre al cargado el 19% de las veces
n =  100: límite 0.240, descubre al cargado el 54% de las veces
n =  300: límite 0.210, descubre al cargado el 94% de las veces
n = 1000: límite 0.190, descubre al cargado el 100% de las veces

• Con 60 tiros, las fracciones del dado justo se alejan bastante de 1/6 ≈ 0,167, y el cargado ni siquiera se delata: en esta simulación su cara más frecuente fue el 2, no el 6. Con 6000, el cargado muestra claramente su 0,25 en el 6 y el justo queda cerca de 0,167 en todas las caras: es la ley de los grandes números.
• Con 100 tiros, un dado justo puede dar entre unos 0,10 y 0,24 de seises solo por azar. Una fracción de 0,22 no prueba nada.
• La regla acusa a un dado justo pocas veces (alrededor del 2 %, el porcentaje que dejamos por encima del límite). Esa es la falsa alarma que aceptamos.
• Con 30 o 100 tiros, el dado cargado se escapa muchas veces. Con 300 tiros se lo descubre en la gran mayoría de los casos, y con 1000 casi siempre. Recomendaríamos unos 300 a 1000 tiros. Hay un equilibrio: cuanto más estricto el límite, menos falsas alarmas, pero hacen falta más tiros para descubrir al tramposo. Es el mismo equilibrio que vas a ver al evaluar modelos de clasificación.