Maksymalna optymalizacja prawdopodobieństwa w przypadku niepewności w x i y została omówiona przez York (2004). Oto kod R jego funkcji.
„YorkFit”, napisany przez Ricka Wehra, 2011, przetłumaczony na R. przez Rachel Chang
Uniwersalna procedura znajdowania najlepszego dopasowania linii prostej do danych ze zmiennymi, skorelowanymi błędami, w tym błędami i oszacowaniami poprawności dopasowania, po równaniu (13) z York 2004, American Journal of Physics, który z kolei był oparty na York 1969, Earth and Planetary Sciences Letters
YorkFit <- funkcja (X, Y, Xstd, Ystd, Ri = 0, b0 = 0, printCoefs = 0, makeLine = 0, eps = 1e-7)
X, Y, Xstd, Ystd: fale zawierające punkty X, punkty Y i ich odchylenia standardowe
OSTRZEŻENIE: Xstd i Ystd nie mogą być zerowe, ponieważ spowoduje to, że Xw lub Yw będą NaN. Zamiast tego użyj bardzo małej wartości.
Ri: współczynniki korelacji dla błędów X i Y - długość 1 lub długość X i Y
b0: wstępne oszacowanie nachylenia (można uzyskać ze standardowego dopasowania najmniejszych kwadratów bez błędów)
printCoefs: ustaw wartość równą 1, aby wyświetlić wyniki w oknie poleceń
makeLine: zestaw równy 1, aby wygenerować falę Y dla linii dopasowania
Zwraca macierz ze znakiem przecięcia i nachylenia oraz ich niepewności
Jeśli nie podano wstępnego przypuszczenia dla b0, po prostu użyj OLS, jeśli (b0 == 0) {b0 = lm (Y ~ X) $ współczynniki [2]}
tol = abs(b0)*eps #the fit will stop iterating when the slope converges to within this value
a, b: końcowy punkt przecięcia i nachylenie a.err, b.err: oszacowane niepewności dotyczące punktu przecięcia i nachylenia
# 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)
}
lm