import ipywidgets
import matplotlib.pyplot as plt
import numpy as np
from scipy import stats

plt.rc("font", size=14)

Variable aleatoria

Una variable aleatoria es una variable que representa las posibles resultados de una medición o experimento. Esta queda completamente caracterizada por su distribución de probabilidad, la probabilidad de que salga cada valor.

Un ejemplo simple es el resultado de tirar un dado (de 6 caras). ¿Cuál es su distribución de probabilidad? ¿Cuál es la probabilidad de que salga cada número?

prob_dado = 1 / 6

plt.bar([1, 2, 3, 4, 5, 6], prob_dado, width=0.2)
plt.xlabel("Número del dado")
plt.ylabel("Probabilidad")
plt.yticks((0, 1 / 6), (0, "1/6"));

En Python, podemos generar dado (una variable aleatoria) con la siguiente función:

dado = stats.randint(low=1, high=7)

dado

<scipy.stats._distn_infrastructure.rv_frozen at 0x7f3c82e31eb0>

Podemos tomar una muestra de esta variable aleatoria, es decir, tirar el dado, con la función .rvs:

dado.rvs()

2

dado.rvs(size=5)

array([1, 2, 6, 3, 3])

x = dado.rvs(size=60)

plt.hist(
    x,
    bins=np.arange(0.5, 7),
     rwidth=0.2,
     weights=np.full(x.shape, 1 / x.size),
    label="Experimental",
)
plt.bar(np.arange(1, 7) + 0.2, prob_dado,
        width=0.2, color="orange", label="Teórico")
plt.grid()
plt.legend(bbox_to_anchor=(1, 1))

<matplotlib.legend.Legend at 0x7f3c82dfd400>

Promedio "teórico"

¿Qué significa calcular el promedio "teórico" de una variable aleatoria? Podemos reinterpretar el cálculo de un promedio en términos de la frecuencia con la que salieron los datos:

Supongamos que tiramos el dado 5 veces y salen $\{x_1, \ldots, x_5\} = \{1, 3, 3, 4, 5\}$.

$$ \begin{align} \bar{x} &= \frac{1}{N} \sum_i^N x_i \\ \text{reemplazando} \\ &= \frac{1}{5} \; (1 + 3 + 3 + 4 + 5) \\ \text{sacando factor común} \\ &= \frac{1}{5} \; (1 \times 1 + 2 \times 3 + 1 \times 4 + 1 \times 5) \\ \text{distribuyendo el N=5} \\ &= \frac{1}{5} \times 1 + \frac{2}{5} \times 3 + \frac{1}{5} \times 4 + \frac{1}{5} \times 5 \\ \text{reinterpretando} \\ \bar{x} &= \sum_k x_k \, f_k \end{align} $$

donde $f_k$ es la fracción de veces que salió cada valor posible $x_k$.

Entonces, reemplazando la frecuencia experimental por la teórica, se puede calcular el promedio teórico $\mu$ como:

$$ \mu = \sum_k x_k \, p_k $$

donde $p_k$ es la probabilidad de que salga $x_k$.

Para el dado:

$$ \begin{align} \mu &= \sum_{i=1}^6 x_i \, \frac{1}{6} \\\\ &= \frac{1}{6} \times 1 + \frac{1}{6} \times 2 + \dots + \frac{1}{6} \times 6 \\\\ &= 3.5 \end{align} $$

Se puede calcular el promedio teórico de la variable aleatoria dado con la función .mean():

dado.mean()

3.5

El promedio "experimental" $\bar{x}$ estima el promedio "teórico" $\mu$:

x.mean()

3.55

donde x es la muestra de $N$ tiradas del dado (x = dado.rvs(size=N)).

Desviación estandar "teórica"

Analogamente, la desviación estándar "experimental" $S$ estima la desviación estandar "teórica" $\sigma$:

$$ \begin{align} S^2 &= \frac{1}{N} \sum_i (x_i - \bar{x})^2 \\ &= \sum_k f_k \; (x_k - \bar{x})^2 \\ \\ \sigma^2 &= \sum_k p_k \; (x_k - \mu)^2 \end{align} $$

dado.std(), x.std()

(1.707825127659933, 1.8113070786957504)

Ley de los grandes números

La ley (teorema) dice (parafraseado):

El promedio (experimental) de un gran número de experimentos se acerca cada vez más al valor esperado (promedio teórico) a medida que aumenta la cantidad de experimentos.

x = dado.rvs(size=10000)

plt.hist(x, bins=np.arange(0.5, 7), rwidth=0.2, weights=np.full(x.shape, 1 / x.size))
plt.bar(np.arange(1, 7) + 0.2, prob_dado, width=0.2, color="orange")

print("Promedio", "\n       Teórico:", dado.mean(), "\n  Experimental:", x.mean())
print("\nDesv. est.", "\n       Teórico:", dado.std(), "\n  Experimental:", x.std())

Promedio 
       Teórico: 3.5 
  Experimental: 3.4994

Desv. est. 
       Teórico: 1.707825127659933 
  Experimental: 1.705109861563178

Con más muestras, cada vez más se va a parecer la distribución experimental a la teórica. Por lo tanto, cualquier calculo que hagamos sobre las muestras, como el promedio o la desviación estandar.

Suma de variables aleatorias

Supongamos que el experimento es tirar dos datos y sumar sus resultados.

dado_1, dado_2 = dado.rvs(), dado.rvs()

print("Dado 1:", dado_1)
print("Dado 2:", dado_2)
print("\nSuma:", dado_1 + dado_2)

Dado 1: 5
Dado 2: 1

Suma: 6

Nuestra variable aleatoria ahora es $Y$, la suma de dos dados:

$$ Y = X_1 + X_2 $$

donde $X_1$ y $X_2$ son dados.

¿Cómo es esta distribución?

1. ¿Cuál es el valor mínimo que puede salir? ¿Y el máximo?

2. ¿Cuál es la probabilidad de cada valor?

N = 100000
dado_1 = dado.rvs(size=N)
dado_2 = dado.rvs(size=N)
x = dado_1 + dado_2

plt.hist(x, bins=np.arange(0.5, 14), rwidth=0.2, weights=np.full(x.shape, 1 / x.size))

# Predicción
plt.bar(np.arange(1, 13) + 0.2,
        1/12,
        width=0.2, color="orange")

plt.ylabel("Probabilidad")

Text(0, 0.5, 'Probabilidad')

Para que salga $Y=2$, tanto $X_1$ como $X_2$ tienen que ser $1$.

En cambio, para $Y=7$, hay múltiples $(x_1, x_2)$ que dan ese resultado: $(1, 6), (2, 5), (3, 4), (4, 3), \ldots$

¿Qué pasará con el promedio y la desviación estandar?

x.mean(), x.std()

(7.00069, 2.424885466140618)

$$ Y = X_1 + X_2 $$

Se puede ver que el promedio de la suma es la suma de los promedios:

$$ \langle Y \rangle = \langle X_1 + X_2 \rangle = \langle X_1 \rangle + \langle X_2 \rangle = 3.5 + 3.5 = 7 $$

Para la desviación estandar, podemos usar la fórmula de propagación*:

$$ \begin{align} \sigma_Y^2 &= \left( \frac{\partial Y}{\partial X_1} \sigma_{X_1} \right)^2 + \left( \frac{\partial Y}{\partial X_2} \sigma_{X_2} \right)^2 \\\\ &= \left( 1 \times \sigma_{X_1} \right)^2 + \left( 1 \times \sigma_{X_2} \right)^2 \\\\ &= 2 \; \sigma_{X}^2 \approx 2 \times (1.70)^2 \approx (2.41)^2 \end{align} $$

Extra

*En general, es una aproximación, pero en el caso de funciones lineales es exacta.

*Solo vale para variables independientes.

Teorema central del límite

La suma de $N$ variables independientes $X_i$ (con varianza finita $\sigma^2 < \infty$):

$$ Y = \sum_i^N X_i $$

tiende a una distribución normal o gaussiana cuando $N$ tiende a infinito.

---

Recordatorio:

Una distribución gaussiana tiene dos parámetros: el centro $\mu$, y el ancho $\sigma$:

$$ G(x) \propto \exp \left( -\frac{(x - \mu)^2}{2 \sigma^2} \right) $$

llamados así porque son el promedio y la desviación estandar "teóricos". Estos se pueden estimar como el promedio y la desviación estándar (experimentales), respectivamente.

---

¿Cuántas variables hay que sumar para que se parezca a una gaussiana?

@ipywidgets.interact(N=(1, 10))
def suma_de_N_dados(N=1):
    # Sumo N dados (10 mil veces)
    x = 0
    for _ in range(N):
        x += dado.rvs(size=10000)

    # Histograma
    plt.hist(
        x,
        bins=np.arange(0.5, 6 * N + 1),
        rwidth=0.5,
        weights=np.full(x.shape, 1 / x.size),
    )

    # Comparación con gaussiana
    gaussiana = stats.norm(loc=x.mean(), scale=x.std())
    t = np.linspace(0, 6 * N + 1, 1000)
    plt.plot(t, gaussiana.pdf(t), lw=3)

Entonces, la gaussiana puede ser una buena aproximación en cuando sumamos múltiples variables aleatorias.

Multiples fuentes de error

Cuando en nuestra medición $M$ de una magnitud cuyo valor real es $M_0$ intervienen multiples fuentes de error $X, Y, Z, \dots$:

$$ M = M_0 + X + Y + Z + \dots $$

es razonable suponer que nuestra medición va a tener una distribucion gaussiana.

El promedio

Otro caso donde es razonable asumir una distribución gaussiana es cuando se calcula un promedio:

$$ \bar{X} = \frac{X_1 + X_2 + \dots + X_N}{N} $$

donde $X_1, X_2, \dots, X_N$ son variables aleatorias de la misma distribución (un dado, por ejemplo)

Por lo tanto, tienen la misma desviación estandar $\sigma_X$.

Si propagamos el error de $\bar{X}$:

$$ \begin{align} \sigma^2_\bar{x} &= \left( \frac{\partial \bar{x}}{\partial X_1} \sigma_{X_1} \right)^2 + \dots + \left( \frac{\partial \bar{x}}{\partial X_N} \sigma_{X_N} \right)^2 \\ \text{reemplazamos derivadas} \\ &= \left( \frac{1}{N} \, \sigma_{X_1} \right)^2 + \dots + \left( \frac{1}{N} \, \sigma_{X_N} \right)^2 \\ \text{usamos } \sigma_{X_i} = \sigma_X \\ &= \left( \frac{1}{N} \, \sigma_X \right)^2 + \dots + \left( \frac{1}{N} \, \sigma_X \right)^2 \\ \text{N términos iguales} \\ &= N \left( \frac{1}{N} \, \sigma_X \right)^2 \\ &= \left( \frac{\sigma_X}{\sqrt{N}} \right)^2 \end{align}$$

Simplificando,

$$ \sigma_\bar{x} = \frac{\sigma_X}{\sqrt{N}} $$

también conocido como error estadístico $\sigma_e$ (del promedio).

Ejemplo de las naranjas

Caso 1: tomamos una única naranja y medimos su diámetro varias veces. No es tan simple medir el diametro de una naranja con una regla, a veces uno se pasa, o mide de menos. ¿Cuál es el diametro de dicha naranja?

Caso 2: tomamos muchas naranjas, y medimos una única vez el diámetro de cada una. ¿Cuál es el diámetro de las naranjas?

En ambos casos, el valor medido $Y$ se puede modelar como la suma de dos variables aleatorias:

$$ Y = D + X $$

donde $D$ sería el diámetro y $X$ el efecto de la medición.

Nosotros queremos caracterizar $D$, pero obtenemos $Y$ cuando medimos.

Asumamos que todas tienen distribución gaussiana (aunque no hace falta):

$ D \sim N(\mu_d, \; \sigma_d^2) $

$ X \sim N(\mu_c, \; \sigma_x^2) $

y, por lo tanto,

$ Y \sim N(\mu_d + \mu_x, \; \sigma_d^2 + \sigma_x^2) $

np.random.seed(42)

D = stats.norm(loc=10, scale=1)
X = stats.norm(loc=0, scale=0.3)

d = D.rvs(size=1000)
x = X.rvs(size=d.size)

td = np.linspace(-1, 1, 100)
td = 5 * td + 10

tx = 3 * np.linspace(-1, 1, 100)

@ipywidgets.interact(i=(0, d.size-1), j=(0, 4))
def animacion(i=0, j=0):
    fig, ax = plt.subplots(1, 2, sharex=True, figsize=(8, 4))
    
    ax[0].plot(td, D.pdf(td), color="C1")
    
    di = d[i]
    yd = D.pdf(di)
    
    if j == 0: return
    ax[0].vlines(di, 0, yd, color="C1")
    ax[0].scatter(di, yd, color="C1")
    
    if j == 1: return
    ax[0].plot(tx + di, yd * X.pdf(tx) / X.pdf(0), color="C0")
    
    if j == 2: return
    xi = x[i]
    ax[0].vlines(di + xi, 0, yd * X.pdf(xi) / X.pdf(0), color="C0")
    ax[0].scatter(di + xi, yd * X.pdf(xi) / X.pdf(0), color="C0")
    
    if j == 3: return
    ax[1].hist(d[:i+1] + x[:i+1], bins=td[::5])

Para describir a las naranjas, queremos describir $D$:

$ D \sim N(\mu_d, \; \sigma_d^2) $

pero tenemos el efecto de la medición:

$ X \sim N(\mu_x, \; \sigma_x^2) $

y lo que medimos es:

$ Y \sim N(\mu_d + \mu_x, \; \sigma_d^2 + \sigma_x^2) $

Una hipotesis razonable es que $\mu_x = 0$, es decir, que cuando medimos, "medimos de menos o de más por igual". Entonces, el promedio las mediciones es el promedio de los diámetros:

$$ \langle Y \rangle = \langle D \rangle $$

Y, si la dispersión de los diámetros $\sigma_d$ es mucho mayor al error de medición $\sigma_x$:

$$ \sigma_y^2 = \sigma_d^2 + \sigma_x^2 \approx \sigma_d^2 $$

sigma_d = 1
sigma_x = 0.1  # 10 veces más chico

np.sqrt(sigma_d**2 + sigma_x**2)

1.004987562112089

@ipywidgets.interact(m_d=(0, 5), m_x=(0, 5),
                     s_d=(0.01, 2), s_x=(0.01, 2))
def naranjas(m_d=10, s_d=1, m_x=0, s_x=1):
    t = np.linspace(-5, 10, 100)

    parameters = {
        "D ~ $N(\mu_d, \sigma_d)$": (m_d, s_d),
        "X ~ $N(\mu_x, \sigma_x)$": (m_x, s_x),
        "Suma": (m_d + m_x, np.sqrt(s_d ** 2 + s_x ** 2)),
    }

    for label, (m, s) in parameters.items():
        gaussian = stats.norm(loc=m, scale=s)
        plt.plot(t, gaussian.pdf(t) / gaussian.pdf(m), label=label)

    plt.legend(bbox_to_anchor=(1, 1))

Extras

Máximo - Mínimo

¿Por qué el $\max(X_i) - \min(X_i)$ es un mal estimador para la dispersión de los datos (de una gaussiana)?

A diferencia de la desviación estándar, el valor de máximo menos mínimo depende de la cantidad de datos que tomemos, y empeora a medida que tomamos más datos.

En la figura de abajo, tomamos $N$ datos de una distribución gaussiana, con valores teóricos $\mu=0, \sigma=1$, y calculamos la desviación estándar y máximo - mínimo. Repetimos este procedimiento 1000 veces, es decir, calculamos 1000 desviaciones estándares de $N$ datos, y gráficamos un histograma para ese $N$. Y volvemos a repetir lo mismo pero para otro valor de $N$.

fig, ax = plt.subplots(1, 2, figsize=(12, 4))

Ns = (10, 30, 100, 300, 1000)

cmap = plt.get_cmap("viridis")
colors = cmap(np.linspace(0, 1, len(Ns)))

for N, color in zip(Ns, colors):
    x = np.random.normal(size=(1000, N))
    ax[0].hist(x.std(1), bins="fd", histtype="step", color=color, label=N)
    ax[1].hist(x.ptp(1), bins="fd", histtype="step", color=color, label=N)

ax[0].set(title="Desviación estándar")
ax[1].set(title="Máximo - Mínimo")

ax[1].legend(title="Cantidad de muestras N", bbox_to_anchor=(1, 1))

<matplotlib.legend.Legend at 0x7f3c82b53250>

Lo que se ve es que la desviación estándar $S$ da siempre alrededor de $1$, pero con menor dispersión a medida que aumenta la cantidad de datos. En cambio, al aumentar la cantidad de datos, el valor de máximo - mínimo aumenta.

La razón de esto es que los valores extremos (alejados del centro) tienen muy baja probabilidad de salir, y empiezan a salir cuando tomamos muchas muestras. En particular, para la gaussiana, el máximo - mínimo tiende a infinto a medida que tomamos más datos.

Ojo, eso no significa que el máximo menos mínimo no sirva nunca. Si tenemos una distribucion uniforme (por ejemplo, un dado), podemos usarlo para estimar su ancho. Pero, obviamente, va a tener una interpretación distinta a la que pensamos con el $\sigma$ de una gaussiana.