Régression lorsque chaque point a sa propre incertitude en

J'ai fait mesures de deux variables et . Ils ont tous deux des incertitudes connues et qui leur sont associées. Je veux trouver la relation entre et . Comment puis-je le faire? $n$ $x$ $y$ $\sigma_x$ $\sigma_y$ $x$ $y$

EDIT : chaque a un différent associé, et le même avec le . $x_i$ $\sigma_{x,i}$ $y_i$

Exemple R reproductible:

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

Le problème avec cet exemple est que je pense qu'il suppose qu'il n'y a pas d'incertitudes dans . Comment puis-je réparer cela? $x$

r regression deming-regression rhombidodécaèdre
la source

lm

Y

$Y$

P (Y | X)

$P(Y | X)$

Y

$Y$

X

$X$

X

$X$

Pour votre cas plutôt spécial (univarié avec un rapport connu de niveaux de bruit pour X et Y), la régression de Deming fera l'affaire, par exemple la Demingfonction dans le package R MethComp .

conjugateprior

@conjugateprior Merci, cela semble prometteur. Je me demande: la régression de Deming fonctionne-t-elle toujours si j'ai une variance différente (mais toujours connue) sur chaque individu x et y? c'est-à-dire si les x sont des longueurs, et j'ai utilisé des règles avec des précisions différentes pour obtenir chaque x

rhombidodécaèdre

Je pense que la façon de le résoudre quand il y a des variances différentes pour chaque mesure est d'utiliser la méthode de York. Quelqu'un sait-il s'il existe une implémentation R de cette méthode?

rhombidodécaèdre

@rhombidodecahedron Voir l'ajustement "avec des erreurs mesurées" dans ma réponse: stats.stackexchange.com/questions/174533/… (qui est tiré de la documentation du démodage de paquet).

Roland

Réponses:

$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 vraie ligne est représentée en bleu pointillé. Le long de celui-ci, les points d'origine sont tracés sous forme de cercles creux. Des flèches grises les relient aux points observés, tracées sous forme de disques noirs solides. La solution est dessinée sous forme d'une ligne rouge continue. Malgré la présence de grands écarts entre les valeurs observées et réelles, la solution est remarquablement proche de la ligne correcte dans cette région.

#
# 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
la source

+1. Pour autant que je sache , cela répond aussi à cet ancien Q: stats.stackexchange.com/questions/178727 ? Nous devrions alors le fermer en double.

amibe dit Réintégrer Monica

En outre, selon mon commentaire à la réponse dans ce fil, il semble que la demingfonction puisse également gérer les erreurs variables. Il devrait probablement donner un ajustement très similaire au vôtre.

amibe dit Réintégrer Monica

Je me demande si le déroulement de la discussion a plus de sens si vous changez la place des 2 paragraphes au-dessus et en dessous de la figure?

gung - Réintégrer Monica

Ce matin, un électeur m'a rappelé que cette question avait été posée et répondue de plusieurs manières, avec un code de travail, il y a plusieurs années sur le site de Mathematica SE .

whuber

Cette solution a-t-elle un nom? et éventuellement une ressource pour une lecture plus approfondie (en dehors du site Mathematica SE, je veux dire)?

JustGettinStarted

L'optimisation du maximum de vraisemblance pour le cas d'incertitudes en x et y a été traitée par York (2004). Voici le code R pour sa fonction.

"YorkFit", écrit par Rick Wehr, 2011, traduit en R par Rachel Chang

Routine universelle pour trouver le meilleur ajustement en ligne droite aux données avec des erreurs corrélées variables, y compris l'erreur et la qualité des estimations de l'ajustement, suivant l'équation. (13) de York 2004, American Journal of Physics, qui était basé à son tour sur York 1969, Earth and Planetary Sciences Letters

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

X, Y, Xstd, Ystd: ondes contenant des points X, des points Y et leurs écarts-types

AVERTISSEMENT: Xstd et Ystd ne peuvent pas être nuls car cela fera que Xw ou Yw soit NaN. Utilisez plutôt une très petite valeur.

Ri: coefficients de corrélation pour les erreurs X et Y - longueur 1 ou longueur de X et Y

b0: estimation initiale approximative de la pente (peut être obtenue à partir d'un ajustement des moindres carrés standard sans erreurs)

printCoefs: mis égal à 1 pour afficher les résultats dans la fenêtre de commande

makeLine: mis égal à 1 pour générer une onde Y pour la ligne d'ajustement

Renvoie une matrice avec l'ordonnée à l'origine et la pente plus leurs incertitudes

Si aucune supposition initiale pour b0 n'est fournie, utilisez simplement OLS si (b0 == 0) {b0 = lm (Y ~ X) $ coefficients [2]}

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

a, b: interception finale et pente a.err, b.err: incertitudes estimées en interception et en pente

# 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
la source

Notez également que le package R "IsoplotR" inclut la fonction york (), donnant les mêmes résultats que le code YorkFit ici.

Steven Wofsy