Regression, wenn jeder Punkt sowohl in

Ich habe $n$ Messungen mit zwei Variablen und . Sie haben beide bekannte Unsicherheiten und die mit ihnen verbunden sind. Ich möchte die Beziehung zwischen und . Wie kann ich es tun? $x$ $y$ $\sigma_x$ $\sigma_y$ $x$ $y$

BEARBEITEN : Mit jedem ist ein anderes und mit dem dasselbe . $x_i$ $\sigma_{x,i}$ $y_i$

Reproduzierbares R-Beispiel:

## 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)

Das Problem bei diesem Beispiel ist, dass ich davon ausgehe, dass es in keine Unsicherheiten gibt . Wie kann ich das beheben? $x$

r regression deming-regression rhombidodecahedron
quelle

Es ist wahr, lmpasst ein lineares Regressionsmodell, das heißt: ein Modell der Erwartung von

Y

$Y$ in Bezug auf , in dem eindeutig zufällig ist und als bekannt betrachtet wird. Um mit Unsicherheiten in umzugehen, benötigen Sie ein anderes Modell.

P (Y | X)

$P(Y | X)$

Y

$Y$

X

$X$

X

$X$

Conjugateprior

Für Ihren eher speziellen Fall (univariat mit einem bekannten Rauschpegelverhältnis für X und Y) reicht die Deming-Regression aus, z. B. die DemingFunktion im R-Paket MethComp .

Conjugateprior

@conjugateprior Danke, das sieht vielversprechend aus. Ich frage mich: Funktioniert Deming-Regression immer noch, wenn ich für jedes einzelne x und y eine andere (aber immer noch bekannte) Varianz habe? dh wenn die x Längen sind und ich Lineale mit unterschiedlichen Genauigkeiten verwendet habe, um jedes x zu erhalten

rhombidodecahedron

Ich denke, der Weg, es zu lösen, wenn es für jede Messung unterschiedliche Varianzen gibt, ist die Verwendung der York-Methode. Weiß zufällig jemand, ob es eine R-Implementierung dieser Methode gibt?

Rhombidodecahedron

@rhombidodecahedron Siehe "mit gemessenen Fehlern", die in meine Antwort passen: stats.stackexchange.com/questions/174533/… (die aus der Dokumentation der Paketdeming entnommen wurde).

Roland

Antworten:

Lassen Sie die wahre Linie , die durch einen Winkel und einen Wert ist, die Menge $L$ $\theta$ $\gamma$

(x, y) : \cos (θ) x + Sünde (θ) y = γ .

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

$(x,y)$

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

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

Die Varianz von sei und die von sei , Unabhängigkeit von $x_i$ $\sigma_i^2$ $y_i$ $\tau_i^2$ $x_i$ $y_i$

Var (d (x_{ich}, y_{ich}; L)) = \cos^{2} (θ) σ_{ich}^{2} + {Sünde}^{2} (θ) τ_{ich}^{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$

Die wahre Linie ist blau gepunktet dargestellt. Entlang sind die ursprünglichen Punkte als hohle Kreise eingezeichnet. Graue Pfeile verbinden sie mit den beobachteten Punkten, die als durchgezogene schwarze Scheiben dargestellt sind. Die Lösung wird als durchgezogene rote Linie gezeichnet. Trotz der großen Abweichungen zwischen den beobachteten und den tatsächlichen Werten ist die Lösung in diesem Bereich bemerkenswert nahe an der korrekten Linie.

#
# 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
quelle

+1. Soweit ich weiß , beantwortet dies auch dieses ältere Q: stats.stackexchange.com/questions/178727 ? Wir sollten es dann als Duplikat schließen.

Amöbe sagt Reinstate Monica

Laut meinem Kommentar zur Antwort in diesem Thread sieht es auch so aus, als ob die demingFunktion auch mit variablen Fehlern umgehen kann. Es sollte wahrscheinlich eine sehr ähnliche Passform ergeben.

Amöbe sagt Reinstate Monica

Ich frage mich, ob der Ablauf der Diskussion sinnvoller ist, wenn Sie die Stellen der beiden Absätze über und unter der Abbildung vertauschen.

gung - Wiedereinsetzung von Monica

Ich wurde heute Morgen (von einem Wähler) daran erinnert, dass diese Frage vor einigen Jahren auf der Mathematica SE-Website auf verschiedene Weise mit Arbeitscode gestellt und beantwortet wurde .

whuber

Hat diese Lösung einen Namen? und möglicherweise eine Quelle für weiteres Lesen (neben der Mathematica SE-Site, meine ich)?

JustGettinStarted

Die Optimierung der maximalen Wahrscheinlichkeit für den Fall von Unsicherheiten in x und y wurde von York (2004) angesprochen. Hier ist R-Code für seine Funktion.

"YorkFit", geschrieben von Rick Wehr, 2011, übersetzt in R von Rachel Chang

Universelle Routine zum Ermitteln der besten geraden Anpassung an Daten mit variablen, korrelierten Fehlern, einschließlich Fehler- und Anpassungsgüte-Schätzungen, nach Gl. (13) von York 2004, American Journal of Physics, das sich wiederum auf York 1969, Earth and Planetary Sciences Letters, stützte

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

X, Y, Xstd, Ystd: Wellen mit X-Punkten, Y-Punkten und deren Standardabweichungen

WARNUNG: Xstd und Ystd dürfen nicht Null sein, da dies dazu führt, dass Xw oder Yw NaN sind. Verwenden Sie stattdessen einen sehr kleinen Wert.

Ri: Korrelationskoeffizienten für X- und Y-Fehler - Länge 1 oder Länge von X und Y

b0: grobe anfängliche Schätzung für die Steigung (kann aus einer fehlerfreien Standardanpassung der kleinsten Quadrate erhalten werden)

printCoefs: Setzen Sie den Wert auf 1, um die Ergebnisse im Befehlsfenster anzuzeigen

makeLine: Gleich 1 setzen, um eine Y-Welle für die Anpassungslinie zu erzeugen

Gibt eine Matrix mit dem Achsenabschnitt und der Steigung sowie deren Unsicherheiten zurück

Wenn keine anfängliche Schätzung für b0 gegeben ist, dann verwende einfach OLS, wenn (b0 == 0) {b0 = 1m (Y ~ X) $ Koeffizienten [2]}

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

a, b: Endschnitt und Steigung a.err, b.err: geschätzte Unsicherheiten in Schnitt und Steigung

# 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
quelle

Beachten Sie auch, dass das R-Paket "IsoplotR" die Funktion york () enthält, die die gleichen Ergebnisse wie der YorkFit-Code liefert.

Steven Wofsy