Ajuste de mínimos cuadrados versus promedio

Question

Supongamos que quiero encontrar la constante elástica de un resorte a partir de los puntos de datos medidos de $F$ y $x$. Hay dos formas básicas de hacer esto.

1. Podría calcular la constante elástica para cada punto de datos y promediar los resultados. $k = \frac{1}{N} \sum_i \frac{F_i}{x_i}$
2. Podría encontrar el ajuste por mínimos cuadrados de una función lineal a los datos, lo cual me daría $k = \frac{\sum F_i x_i}{\sum x_i^2}$.

¿Hay alguna razón para preferir uno de estos métodos sobre el otro? He estado enseñando a mis estudiantes de secundaria a hacer el segundo como una introducción simple a las técnicas de ajuste de datos, pero no tengo una buena respuesta para explicarles por qué no deberían simplemente hacer el primero, que les parece más obvio.

Answer 1

Esa es una buena pregunta. La respuesta depende de cómo se recopilan los datos. Vamos a plantear dos modelos simples y proporcionar el mejor estimador no sesgado para cada modelo a través del Teorema de Gauss-Márkov.

Consideraremos que el desplazamiento $x$ está fijado de antemano y que la fuerza $F$ se mide con error. Cada modelo asume que el error es aditivo, con media cero e independiente. Sin embargo, cada modelo describirá de manera diferente la variabilidad del error.

---

Primero, supongamos que la fuerza $F$ se mide con una varianza de error independiente del desplazamiento $x$. Entonces, GM prueba que $\frac{\sum_{i=1}^n x_iF_i}{\sum_{i=1}^n x_i^2}$ es el estimador no sesgado de menor varianza. Por lo tanto, en este caso tienes razón en que tus estudiantes deben usar tu estimador.

Segundo, supongamos que la fuerza $F$ se mide con una desviación estándar de error proporcional a $x$. Por lo tanto, las fuerzas más grandes se miden con más variabilidad. Luego, GM prueba que $\frac{1}{n} \sum_{i=1}^n \frac{F_i}{x_i}$ es el estimador no sesgado de menor varianza. Por lo tanto, en este caso la intuición de tus estudiantes es acertada.

---

En conclusión, qué estimador es mejor depende de en qué medida la varianza del error depende del desplazamiento $x$. ¿Crees que una fuerza dada y una fuerza cien veces mayor se miden con la misma precisión? Lo dudo. En la práctica, creo que la desviación estándar es proporcional a $x$ (y por lo tanto a $F$). Por lo tanto, usaría el estimador que prefieren tus estudiantes.

Answer 2

Aquí está el enfoque cuantitativo (simulación) que mencionó @Dr. Momo: podemos simular algunos puntos para $x_i$ y $F_i$, asignarles algunos errores gaussianos, estimar $k$ con los dos enfoques y ver qué tan cerca están las estimaciones del verdadero parámetro.

Los supuestos subyacentes a esta simulación son probablemente una buena aproximación para experimentos simples de laboratorio, pero uno debe tener cuidado al elegir el método a aplicar a un problema específico para verificar cómo se distribuyen las incertidumbres; no todo es gaussiano en realidad, y descifrar estas distribuciones es un asunto complicado.

Dicho eso: estos dos estimadores, aunque ambos sin sesgo, difieren en su varianza (@usuario551504 ha discutido la razón por la cual).

El resultado se ve así: la normalización es irrelevante, pero se puede ver claramente que los errores son en promedio más pequeños con el enfoque de ajuste lineal.

Aquí está el código (en python - también requiere los paquetes numpy y matplotlib):

import numpy as np
import matplotlib.pyplot as plt

# desviaciones estándar de los errores de x y F
sigma_x = .1
sigma_F = .1

# parámetro elástico verdadero
k_true = 1.

# número de puntos de medida ficticios
n_measurement_points = 20

# número de pruebas a realizar
n_trials = 1000

def trial(seed: int):

    # generar puntos x igualmente espaciados entre 1 y 20, inclusive
    xi = np.linspace(1, 20, num=n_measurement_points)
    # generar las fuerzas "verdaderas", de acuerdo con la ley lineal
    # (sin el signo negativo por simplicidad, no cambia nada)
    Fi = k_true * xi

    # agregar ruido aleatorio gaussiano a x y F
    rng = np.random.default_rng(seed)    
    xi += rng.normal(scale=sigma_x, size=xi.shape)
    Fi += rng.normal(scale=sigma_F, size=Fi.shape)

    # estimar el parámetro elástico con los dos enfoques
    k_estimate_avg = 1 / n_measurement_points * np.sum(Fi / xi)
    k_estimate_fit =  np.sum(Fi * xi) / np.sum(xi**2)

    return k_estimate_avg - k_true, k_estimate_fit - k_true

# realizar la prueba varias veces para obtener estadísticas sobre la varianza
# del estimador
avg_errs = []
fit_errs = []
for n in range(n_trials):
    # sembramos el RNG con el índice del bucle for - no importa realmente,
    # es útil para la reproducibilidad
    err_avg, err_fit = trial(n)
    avg_errs.append(err_avg)
    fit_errs.append(err_fit)

# hacer un histograma de los resultados
plt.hist(avg_errs, alpha=.5, density=True, label='promedio')
plt.hist(fit_errs, alpha=.5, density=True, label='ajuste')
plt.xlabel('Estimación menos verdadero $k$')
plt.ylabel('Densidad de probabilidad')
plt.legend()

# para guardar la figura, se puede hacer
# plt.savefig('ajuste_lineal_vs_avg.png', dpi=150)

# para mostrarla de forma interactiva, en cambio, hacer
# plt.show()

También señalaré que el resultado es cualitativamente el mismo incluso si los errores solo se incluyen para $x$ o para $F$.

Answer 3

Como en todas las cosas estadísticas, todo depende completamente de tus suposiciones.

Por qué los promedios son útiles: Si tienes una variable x de la cual puedes tomar muestras, y proviene de una distribución normal con alguna media $\mu$ y desviación estándar $\sigma$, resulta que el estimador de máxima verosimilitud para la verdadera media $\mu$ es el promedio de la muestra $\frac{1}{N} \sum x_{i}$: el promedio de la muestra es tu mejor estimación de la verdadera media. Un ejemplo de este proceso es un tablero de gacha muy grande (o tabla de Galton). Si no sabes de dónde caen las pelotas, puedes estimarlo midiendo las ubicaciones donde terminan las pelotas. Para introducir brevemente las matemáticas relevantes, cuando intentas maximizar el producto de las Gaussianas que representan la probabilidad de que $\mu$ sea la verdadera media: $\Pi_{i} \exp{(x_{i}-\mu)^{2}}$ y tomas una derivada igual a 0, obtienes que $\sum (x_{i} - \mu) = 0$ entonces $\mu = \frac{1}{N}\sum x_{i}$ es la suposición de mayor verosimilitud para la media.

Por qué los mínimos cuadrados son útiles: Esto está en el nombre. Minimiza el error cuadrático mínimo. Este se convierte en un estimador de máxima verosimilitud cuando tus residuos $\epsilon = F - kx$ están distribuidos de forma normal: ocurriendo con una probabilidad relativa $P(\epsilon) \propto \exp{(-\epsilon^{2}/\sigma^{2})}$. Esta vez, al intentar estimar $k$ para maximizar la verosimilitud, tomamos derivadas con respecto a $k$ que está dentro de $\epsilon = F - kx$. Cada término del producto contribuirá con un $2\epsilon \frac{d\epsilon}{dk} = 2\epsilon x = 2 Fx - 2 kxx$. Nuevamente, la suma de estos nos dará la mejor suposición para k: $\sum F_{i}x_{i} - \sum k x_{i}^{2} = 0$ entonces $k = \frac{\sum F_{i}x_{i}}{\sum x_{i}^{2}}$

Entonces la pregunta es, ¿crees que cada vez que tomas una muestra de $F_{i}, x_{i}$ estás tomando una muestra de $k_{i} = F_{i}/x_{i}$ que proviene de una distribución normal centrada alrededor de la verdadera constante k, o crees que el resorte tiene una constante $k$ que te dará residuos distribuidos de forma normal en $F_{i}$, $x_{i}$?

Creo que una gran lección para tus estudiantes es no solo tomar promedios o hacer mínimos cuadrados por instinto, sino comprender que estas son maquinarias matemáticas comúnmente utilizadas por una razón, derivadas por la razón, y vienen con suposiciones.

Answer 4

El enfoque de linealidad tiene algunas ventajas:

1. Da menos peso a los puntos cercanos al origen. Esto suele ser algo positivo, ya que en general esos puntos tienen barras de error relativamente más grandes.

2. Suaviza el ruido de manera más efectiva y se ve menos afectado por los valores atípicos que el primer enfoque.

Puedes experimentar con esto en Excel y ver cómo difieren los dos enfoques. Estoy seguro de que alguien puede ofrecer más información cuantitativa, pero esa es mi experiencia con los dos enfoques.

Como advertencia, si estás utilizando el enfoque de linealidad, debes recopilar datos en intervalos iguales a lo largo de tu dominio. Incluso profesionalmente, a menudo veo gráficos de linealidad con 4-5 puntos agrupados en x = 1 y un solo punto en x = 10. Esto demuestra muy poco sobre la linealidad del sistema.

Answer 5

Creo que la respuesta depende de lo que quieras hacer con la constante elástica que extraigas, pero en general creo que el segundo enfoque es más útil.

Presumiblemente deseas un modelo simple de un resorte que produzca una fuerza (a través de la ley de Hooke) con un pequeño error, dado el desplazamiento. A menos que te importe un rango específico de desplazamientos, esto es precisamente lo que te ofrece el segundo enfoque, minimizando el error en el sentido de los mínimos cuadrados.

El primer enfoque todavía minimiza algún error, pero uno que puede que no sea tan útil. Minimiza el error en $F/x$ en el sentido de los mínimos cuadrados, lo que, como señala el Dr. Momo, enfatiza los errores en $F$ para desplazamientos pequeños. Esto significa que la constante elástica extraída de esta manera es un peor modelo para desplazamientos grandes.

Si hay un rango específico de desplazamientos que te importan más, también podrías usar una función de ponderación $w(x)$ que priorice los desplazamientos más importantes, obteniendo $$k = \frac{\sum w^2(x_i) F_i x_i}{\sum w^2(x_i) x_i^2}$$ que minimiza el error en $w(x)F$ en el sentido de los mínimos cuadrados. $w(x) = 1/x$ para tu primer enfoque y $w(x) = 1$ para el segundo.