Cálculo de la intercepción en estimador de Theil-Sen

Question

Acabo de escribir algo de código C para el Theil-Sen, después de buscar un poco en Google (yo no tengo ninguna documentación definitiva). Mi comprensión de la intersección de cálculo es que por primera vez me calcular la mediana pendiente, y luego construir una línea a través de cada punto de datos con esta pendiente, encontrar la intersección de cada línea y, a continuación, tomar la mediana de interceptar.

La única manera en que puedo encontrar para probar el código es comparar los resultados con los de Kendall-Theil Robusta Línea de programa, a partir de la USGS. En un conjunto de datos de los 237 puntos (datos de salud, con una correlación de Pearson de ~0.55), estamos de acuerdo exactamente en la mediana de la pendiente, pero no está de acuerdo en el intercepto (1,4%). Según mis datos, la KTRL interceptar no es la mediana de interceptar, pero es en lugar de 46% de la forma a través de la gama.

Después de algunos escarbando en la KTRL código, parece que calcular la intersección mediante la creación de una sola "línea media", en lugar de la mediana de todos los intercepta. Su intercepción es medianY - medianX * median slope.

Cualquier comentario sobre cual es la forma "correcta" de hacer esto, si hay uno, o cómo se maneja este asunto en R/etc?

Gracias.

Answer 1

El Theil-Sen estimador es esencialmente un estimador de la pendiente solo; la línea que se ha construido en una variedad de maneras diferentes.

De hecho, hay una gran variedad de formas para calcular la intersección.

Usted dijo:

Mi comprensión de la intersección de cálculo es que por primera vez me calcular la mediana pendiente, y luego construir una línea a través de cada punto de datos con esta pendiente, encontrar la intersección de cada línea y, a continuación, tomar la mediana de interceptar.

Esa es una manera en que se hace.

Hay un par de métodos que calculan los tipos de ponderado medianas, pero en base a ese mismo enfoque (poner más peso en los puntos más separados en el x-espacio).

Otro común (probablemente la más común) es calcular la mediana($y-bx$). Esto es lo que Sen miró, por ejemplo.

Otra es tratar de obtener un estimador con la mayor eficiencia posible en el normal (semejante a la del estimador de la pendiente en situaciones típicas) y desglose similar punto de la inclinación de la estimación (probablemente hay poco punto en tener un mejor desglose a expensas de la eficiencia), como el uso de la Hodges-Lehmann estimador (mediana de pares de medias) en $y-bx$. Este tiene una especie de simetría en la forma de las laderas y las intersecciones son definidos ... y generalmente se da algo muy cerca de la LS de línea cuando la normal supuestos cerca de retención, mientras que la Sen-interceptar puede ser - relativamente hablando - muy diferentes.

Y algunas personas simplemente calcular la media.

Todavía hay otras sugerencias que se han mirado. Realmente, no hay " una " interceptar a ir con la pendiente de la estimación.

Dietz listas de varias posibilidades, posiblemente incluso de todos los que he mencionado, pero que de ninguna manera exhaustiva.

Answer 2

Kendall–Theil Robusto programa de Línea de la USGS tiene un compañero PDF.

En la página 8 (PDF, en la página 15) indica el método utilizado y la fórmula que encontró, pero da la referencia como Conover.

Interceptar

La estimación de la intersección se calcula mediante el uso de la Conover (1980) ecuación

$$b = Y_\text{median} -m\times X_\text{median}\ \ ,\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ (6)$$

donde
$\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ b\ \ \ \ \ $ es la estimación de la intersección,
$\ \ \ \ \ \ \ Y_\text{median}\ \ \ \ $ es la media de las variables de respuesta,
$\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ m\ \ \ \ $ es la estimación de la pendiente,
y
$\ \ \ \ \ \ X_\text{median}\ \ \ \ $ es la media de las variables explicativas.

Puedo confirmar esto no produce el mismo resultado que el programa. Si hay métodos superiores y así sucesivamente, como siempre, una cuestión de opinión y de sus circunstancias particulares.

El M-algoritmo de estimación es probablemente errónea.

for i = 1,   # dat-1 do
for j = i+1, # dat do

Cambio j de indización de i+1 y no de instancias de proceso de i == j.

Luego clasificar (sort) el resultado de la elección de índice como se describe, o podría decirse que tomar la mediana, la cual le dará un poco diferente resultado. Si usted trama, los datos se verá como una CDF de la parcela.