Regresión cuando cada punto tiene su propia incertidumbre tanto en

Hice mediciones de dos variables e . Ambos tienen incertidumbres conocidas y asociadas con ellos. Quiero encontrar la relación entre e . ¿Cómo puedo hacerlo? $n$ $x$ $y$ $\sigma_x$ $\sigma_y$ $x$ $y$

EDITAR : cada tiene una asociada a ella, y lo mismo con la . $x_i$ $\sigma_{x,i}$ $y_i$

Ejemplo R reproducible:

## pick some real x and y values 
true_x <- 1:100
true_y <- 2*true_x+1

## pick the uncertainty on them
sigma_x <- runif(length(true_x), 1, 10) # 10
sigma_y <- runif(length(true_y), 1, 15) # 15

## perturb both x and y with noise 
noisy_x <- rnorm(length(true_x), true_x, sigma_x)
noisy_y <- rnorm(length(true_y), true_y, sigma_y)

## make a plot 
plot(NA, xlab="x", ylab="y",
    xlim=range(noisy_x-sigma_x, noisy_x+sigma_x), 
    ylim=range(noisy_y-sigma_y, noisy_y+sigma_y))
arrows(noisy_x, noisy_y-sigma_y, 
       noisy_x, noisy_y+sigma_y, 
       length=0, angle=90, code=3, col="darkgray")
arrows(noisy_x-sigma_x, noisy_y,
       noisy_x+sigma_x, noisy_y,
       length=0, angle=90, code=3, col="darkgray")
points(noisy_y ~ noisy_x)

## fit a line 
mdl <- lm(noisy_y ~ noisy_x)
abline(mdl)

## show confidence interval around line 
newXs <- seq(-100, 200, 1)
prd <- predict(mdl, newdata=data.frame(noisy_x=newXs), 
    interval=c('confidence'), level=0.99, type='response')
lines(newXs, prd[,2], col='black', lty=3)
lines(newXs, prd[,3], col='black', lty=3)

El problema con este ejemplo es que creo que supone que no hay incertidumbres en . ¿Cómo puedo arreglar esto? $x$

r regression deming-regression romdodecaedro
fuente

Es cierto que se lmajusta a un modelo de regresión lineal, es decir: un modelo de la expectativa de con respecto a , en el que claramente es tan aleatorio y se considera conocido. Para lidiar con la incertidumbre en , necesitará un modelo diferente.

Y

$Y$

P (Y | X)

$P(Y | X)$

Y

$Y$

X

$X$

X

$X$

conjugateprior

Para su caso bastante especial (univariante con una relación conocida de niveles de ruido para X e Y), la regresión de Deming será suficiente, por ejemplo, la Demingfunción en el paquete R MethComp .

conjugateprior

@conjugateprior Gracias, esto parece prometedor. Me pregunto: ¿la regresión de Deming todavía funciona si tengo una varianza diferente (pero aún conocida) en cada x e y individual? es decir, si las x son longitudes, y usé reglas con diferentes precisiones para obtener cada x

rombodecaedro

Creo que quizás la forma de resolverlo cuando existen diferentes variaciones para cada medición es usar el método de York. ¿Alguien sabe si hay una implementación R de este método?

rombodecaedro

@rhombidodecahedron Vea el ajuste "con errores medidos" en mi respuesta allí: stats.stackexchange.com/questions/174533/… (que se tomó de la documentación del paquete de Deming).

Roland

Respuestas:

$L$ $\theta$ $\gamma$

(x, y) : \cos (θ) x + \sin (θ) y = γ .

$(x,y): \cos(\theta) x + \sin(\theta) y = \gamma.$

$(x,y)$

d (x, y; L) = \cos (θ) x + \sin (θ) y - γ .

$d(x,y;L) = \cos(\theta) x + \sin(\theta) y - \gamma.$

$x_i$ $\sigma_i^2$ $y_i$ $\tau_i^2$ $x_i$ $y_i$

Var (d (x_{i}, y_{i}; L)) = \cos^{2} (θ) σ_{i}^{2} + \sin^{2} (θ) τ_{i}^{2} .

$\operatorname{Var}(d(x_i,y_i;L)) = \cos^2(\theta)\sigma_i^2 + \sin^2(\theta)\tau_i^2.$

$\theta$ $\gamma$

$\sigma_i$ $\tau_i$ $0$

$\tau_i$ $\sigma_i$ $x$ $n=8$

La línea verdadera se muestra en azul punteado. A lo largo, los puntos originales se trazan como círculos huecos. Las flechas grises los conectan a los puntos observados, trazados como discos negros sólidos. La solución se dibuja como una línea roja continua. A pesar de la presencia de grandes desviaciones entre los valores observados y los reales, la solución está notablemente cerca de la línea correcta dentro de esta región.

#
# Generate data.
#
theta <- c(1, -2, 3) # The line is theta %*% c(x,y,-1) == 0
theta[-3] <- theta[-3]/sqrt(crossprod(theta[-3]))
n <- 8
set.seed(17)
sigma <- rexp(n, 1/2)
tau <- rexp(n, 1)
u <- 1:n
xy.0 <- t(outer(c(-theta[2], theta[1]), 0:(n-1)) + c(theta[3]/theta[1], 0))
xy <- xy.0 + cbind(rnorm(n, sd=sigma), rnorm(n, sd=tau))
#
# Fit a line.
#
x <- xy[, 1]
y <- xy[, 2]
f <- function(phi) { # Negative log likelihood, up to an additive constant
  a <- phi[1]
  gamma <- phi[2]
  sum((x*cos(a) + y*sin(a) - gamma)^2 / ((sigma*cos(a))^2 + (tau*sin(a))^2))/2
}
fit <- lm(y ~ x) # Yields starting estimates
slope <- coef(fit)[2]
theta.0 <- atan2(1, -slope)
gamma.0 <- coef(fit)[1] / sqrt(1 + slope^2)
sol <- nlm(f,c(theta.0, gamma.0))
#
# Plot the data and the fit.
#
theta.hat <- sol$estimate[1] %% (2*pi)
gamma.hat <- sol$estimate[2]
plot(rbind(xy.0, xy), type="n", xlab="x", ylab="y")
invisible(sapply(1:n, function(i) 
  arrows(xy.0[i,1], xy.0[i,2], xy[i,1], xy[i,2], 
         length=0.15, angle=20, col="Gray")))
points(xy.0)
points(xy, pch=16)
abline(c(theta[3] / theta[2], -theta[1]/theta[2]), col="Blue", lwd=2, lty=3)
abline(c(gamma.hat / sin(theta.hat), -1/tan(theta.hat)), col="Red", lwd=2)

whuber
fuente

+1. Por lo que yo entiendo, esto también responde a esta pregunta anterior: stats.stackexchange.com/questions/178727 ? Deberíamos cerrarlo como un duplicado entonces.

ameba dice Reinstate Monica

Además, según mi comentario a la respuesta en ese hilo, parece que la demingfunción también puede manejar errores variables. Probablemente debería producir un ajuste muy similar al tuyo.

ameba dice Reinstate Monica

Me pregunto si el flujo de la discusión tiene más sentido si cambia los lugares de los 2 párrafos arriba y abajo de la figura.

gung - Restablecer Monica

Esta mañana, un votante me recordó que esta pregunta había sido formulada y respondida de varias maneras, con un código de trabajo, hace varios años en el sitio de Mathematica SE .

whuber

¿Esta solución tiene un nombre? y posiblemente un recurso para lecturas adicionales (además del sitio de Mathematica SE, quiero decir)?

JustGettinStarted

York (2004) ha abordado la optimización de máxima probabilidad para el caso de incertidumbres en x e y. Aquí está el código R para su función.

"YorkFit", escrito por Rick Wehr, 2011, traducido a R por Rachel Chang

Rutina universal para encontrar el mejor ajuste en línea recta a los datos con errores correlacionados y variables, incluyendo errores y estimaciones de bondad de ajuste, siguiendo la ecuación (13) de York 2004, American Journal of Physics, que se basó a su vez en York 1969, Earth and Planetary Sciences Letters

YorkFit <- función (X, Y, Xstd, Ystd, Ri = 0, b0 = 0, printCoefs = 0, makeLine = 0, eps = 1e-7)

X, Y, Xstd, Ystd: ondas que contienen puntos X, puntos Y y sus desviaciones estándar

ADVERTENCIA: Xstd e Ystd no pueden ser cero ya que esto hará que Xw o Yw sean NaN. Use un valor muy pequeño en su lugar.

Ri: coeficientes de correlación para errores X e Y - longitud 1 o longitud de X e Y

b0: aproximación inicial aproximada para la pendiente (se puede obtener de un ajuste estándar de mínimos cuadrados sin errores)

printCoefs: establecido igual a 1 para mostrar resultados en la ventana de comandos

makeLine: establece igual a 1 para generar una onda Y para la línea de ajuste

Devuelve una matriz con la intersección y la pendiente más sus incertidumbres

Si no se proporciona una suposición inicial para b0, simplemente use OLS if (b0 == 0) {b0 = lm (Y ~ X) $ coeficientes [2]}

tol = abs(b0)*eps #the fit will stop iterating when the slope converges to within this value

a, b: intersección final y pendiente a.err, b.err: incertidumbres estimadas en intersección y pendiente

# WAVE DEFINITIONS #

Xw = 1/(Xstd^2) #X weights
Yw = 1/(Ystd^2) #Y weights


# ITERATIVE CALCULATION OF SLOPE AND INTERCEPT #

b = b0
b.diff = tol + 1
while(b.diff>tol)
{
    b.old = b
    alpha.i = sqrt(Xw*Yw)
    Wi = (Xw*Yw)/((b^2)*Yw + Xw - 2*b*Ri*alpha.i)
    WiX = Wi*X
    WiY = Wi*Y
    sumWiX = sum(WiX, na.rm = TRUE)
    sumWiY = sum(WiY, na.rm = TRUE)
    sumWi = sum(Wi, na.rm = TRUE)
    Xbar = sumWiX/sumWi
    Ybar = sumWiY/sumWi
    Ui = X - Xbar
    Vi = Y - Ybar

    Bi = Wi*((Ui/Yw) + (b*Vi/Xw) - (b*Ui+Vi)*Ri/alpha.i)
    wTOPint = Bi*Wi*Vi
    wBOTint = Bi*Wi*Ui
    sumTOP = sum(wTOPint, na.rm=TRUE)
    sumBOT = sum(wBOTint, na.rm=TRUE)
    b = sumTOP/sumBOT

    b.diff = abs(b-b.old)
  }     

   a = Ybar - b*Xbar
   wYorkFitCoefs = c(a,b)

# ERROR CALCULATION #

Xadj = Xbar + Bi
WiXadj = Wi*Xadj
sumWiXadj = sum(WiXadj, na.rm=TRUE)
Xadjbar = sumWiXadj/sumWi
Uadj = Xadj - Xadjbar
wErrorTerm = Wi*Uadj*Uadj
errorSum = sum(wErrorTerm, na.rm=TRUE)
b.err = sqrt(1/errorSum)
a.err = sqrt((1/sumWi) + (Xadjbar^2)*(b.err^2))
wYorkFitErrors = c(a.err,b.err)

# GOODNESS OF FIT CALCULATION #
lgth = length(X)
wSint = Wi*(Y - b*X - a)^2
sumSint = sum(wSint, na.rm=TRUE)
wYorkGOF = c(sumSint/(lgth-2),sqrt(2/(lgth-2))) #GOF (should equal 1 if assumptions are valid), #standard error in GOF

# OPTIONAL OUTPUTS #

if(printCoefs==1)
 {
    print(paste("intercept = ", a, " +/- ", a.err, sep=""))
    print(paste("slope = ", b, " +/- ", b.err, sep=""))
  }
if(makeLine==1)
 {
    wYorkFitLine = a + b*X
  }
 ans=rbind(c(a,a.err),c(b, b.err)); dimnames(ans)=list(c("Int","Slope"),c("Value","Sigma"))
return(ans)
 }

Steven Wofsy
fuente

También tenga en cuenta que el paquete R "IsoplotR" incluye la función york (), dando los mismos resultados que el código YorkFit aquí.

Steven Wofsy