Regresja liniowa wielokrotna
x = (x, z)
E(Y | X = x) =�(x) = a0 + a1x + a2z
Funkcja kryterialna:
n
S(a0, a1, a2) = �(xi ) - yi 2
( )
"
i=1
n
S(a0, a1, a2) = a0 + a1xi + a2zi - yi 2
( )
"
i=1
Kryterium wyznaczenia
rozwiązania a0, a1, a2:
S(a0,a1,a2)
min
a0 ,a1,a2
n n
"S " �ł 2 �ł "
2
= a0 + a1xi + a2zi - yi �ł = a0 + a1xi + a2zi - yi
( ) ( )
" " ( )
"a0 "a0 �ł "a0
�ł i=1 łł i=1
n n
"
=
( ) ( ) ( )
"2 a0 + a1xi + a2zi - yi "a0 a0 + a1xi + a2zi - yi = "2 a0 + a1xi + a2zi - yi 1
i=1 i=1
n n n n n
= yi =
( )
"2 a0 + a1xi + a2zi - yi = 2a0"1+ 2a1"x + 2a2"z - 2"
i i
i=1 i=1 i=1 i=1 i=1
n n n
= 2na0 + 2a1 i + 2a2 i - 2 yi = 0
"x "z "
i=1 i=1 i=1
n n n
na0 + a1 i + a2 i = yi
"x "z "
i=1 i=1 i=1
n n
"S " �ł 2 �ł "
2
= a0 + a1xi + a2zi - yi �ł = a0 + a1xi + a2zi - yi
( ) ( )
" " ( )
"a1 "a1 �ł "a1
�ł i=1 łł i=1
n n
"
=
( ) ( ) ( )
"2 a0 + a1xi + a2zi - yi "a1 a0 + a1xi + a2zi - yi = "2 a0 + a1xi + a2zi - yi xi
i=1 i=1
n n n n
2
= 2a0 yixi = 0
"x + 2a1"x + 2a2"z xi - 2"
i i i
i=1 i=1 i=1 i=1
n n n n
2
a0 yixi
"x + a1"x + a2"z xi = "
i i i
i=1 i=1 i=1 i=1
n n
"S " �ł 2 �ł "
2
= a0 + a1xi + a2zi - yi �ł = a0 + a1xi + a2zi - yi
( ) ( )
" " ( )
"a2 "a2 �ł "a2
�ł i=1 łł i=1
n n
"
=
( ) ( ) ( )
"2 a0 + a1xi + a2zi - yi "a2 a0 + a1xi + a2zi - yi = "2 a0 + a1xi + a2zi - yi zi
i=1 i=1
n n n n
= 2a0 zi + 2a1 ixi + 2a2 i2 - 2 yizi = 0
" "z "z "
i=1 i=1 i=1 i=1
n n n n
a0 i + a1 zixi + a2 zi2 = yizi
"z " " "
i=1 i=1 i=1 i=1
n n n
na0 + a1 yi
"x + a2"z = "
i i
i=1 i=1 i=1
n n n n
2
a0 yixi
"x + a1"x + a2"z xi = "
i i i
i=1 i=1 i=1 i=1
n n n n
2
a0 yizi
"z + a1"z xi + a2"z = "
i i i
i=1 i=1 i=1 i=1
n n n
�ł �ł łł
n yi śł
"x "z łł "
i i
�ł śł �ł
i=1 i=1 i=1
�ł śł �ł śł
a0
�ł łł
n n n n
�ł śł
2
� yixi śł
"x "x "z xi śł �ła śł = �ł"
i i i 1
�ł śł �ł
�ł śł
i=1 i=1 i=1 i=1
�ł śł �ł śł
�ł śł
�ła2 �ł
n n n n
�ł śł �ł śł
2
yizi śł
"z "z xi "z �ł "
�ł i i i śł �ł
�ł i=1 i=1 i=1 �ł i=1 �ł
Met. Najmniejszych Kwadratów dla
pomiarów bezpośrednich i pośrednich
Postulat:
Wynik kolejnego pomiaru yj można uważać za sumę wielkości
x  nieznanej (poszukiwanej wielkości podlegającej pomiarowi)
oraz
�j - błędu pomiarowego
czyli yj = x + �j dla każdego j=1,...,n
Metoda Najmniejszych Kwadratów postuluje dobieranie tak
wielkości �j , aby suma kwadratów błędów �j była najmniejsza:
n n
2
2
x - yj minimum
( )
"� = "
j
j=1 j=1
Pomiary bezpośrednie
Przypadek równych błędów pomiarowych (pierwiastek z wariancji �2 )
Wykonujemy n pomiarów nieznanej wielkości x. Wyniki pomiarów obarczone są
błędem pomiarowym �j , o którym zakładamy że ma rozkład normalny z wartością
średnią równą zeru:
yj = x + �j E(�j)=0 (wart. oczekiwana) E(�j2)=�2 (wariancja)
(założenie to  o rozkładzie  jest często uzasadnione przez centralne twierdzenie
graniczne).
Najlepszym estymatorem dla x jest średnia arytmetyczną:
n
1
v
x = y = yj
"
n
j=1
którego wariancja wyraża się
v
�2(x) = �2 / n
a identyfikując błąd z odchyleniem standardowym otrzymujemy
�x = �2 / n
Przypadek różnych błędów pomiarowych (pierwiastków z wariancji �j2 )
Wykonujemy n pomiarów nieznanej wielkości x. Wyniki pomiarów obarczone są
błędem pomiarowym �j , o którym zakładamy że ma rozkład normalny z wartością
średnią równą zeru:
yj = x + �j E(�j)=0 (wart. oczekiwana) E(�j2)=�j2=1/gj (wariancja)
Najlepszym estymatorem dla x jest (średnia ważona pomiarów):
n n
v
x =
"g yj /"g
j j
j=1 j=1
którego wariancja wyraża się
-1
-1
n n
�ł �ł
�ł �ł
1
v
�2(x) = =
�ł �ł
�ł �ł
" "g łł
j
�ł
�2 �ł
j=1 j=1
j
�ł
�ł łł
a estymatory błędów pomiarowych
v
� = yj - x
j
mają rozkład normalny z wartością średnią równą zeru i wariancją �j2 , a to oznacza
v
� / �
że wielkości pochodzą ze standardowego rozkładu Gaussa N(0,1), co
j j
powoduje że suma ich kwadratów
n n n
2 2
v v
! = � / � = (yj - x) / � = g
(v ) ( )
" " "(y - x)2
j j j j j
j=1 j=1 j=1
ma rozkład �2 o n-1 stopniach swobody. Własności te można wykorzystać do
weryfikacji słuszności założeń przy użyciu testu �2 .
Przykład. Najlepsze wartości stałych fizycznych uzyskuje się przez obliczanie średniej
ważonej wszystkich pomiarów przeprowadzonych przez różne ośrodki naukowe-
doświadczalne. Dla cząstek elementarnych takie średnie są wyznaczane regularnie,
bowiem coraz bardziej doskonalone są metody pomiarowe. Np. średnia masa
neutralnego mezonu K uśredniono w oparciu o wyniki 4 doświadczeń wykonanych
różnymi technikami doświadczalnymi (tabela).
�j 1/(�j)^2=Gj
j Yj Yj*Gj Yj-x Gj(Yj-x)^2
1 498.10 0.40 6.30 3038.0 0.20 0.3
2 497.44 0.33 10.00 4974.4 -0.46 2.1
3 498.90 0.50 4.00 1995.6 1.00 4.0
4 497.44 0.50 4.00 1989.8 -0.46 0.8
Ł 24.30 11997.8 7.2
estymata x = ŁYjGj/ŁGj = 497.9 bład �x = 1/(ŁGj)=0.20
Wartość !min wynosi w tym przykładzie 7.2 (tabela) . Dla poziomu istotności 5%, oraz
stopni swobody n-1=4-1=3 wartość krytyczna �2=7.95. Wniosek: uważamy ze
uzyskany wynik ( masa mezonu K = 497.9 ą0.2 MeV ) stanowi najlepszą wartość z
poziomem ufności 95%. Tak pozostało do czasu (upłynęło 20 lat) gdy przeprowadzono
nowe doświadczenia i uzyskano wynik 497.671 ą0.030 MeV .
UWAGA: jeśli test �2 da rezultat negatywny (odrzucamy hipotezę o dobrym
dopasowaniu bowiem �2 > �2 ) powodem tego może być fakt, że co najmniej jeden z
kryt
pomiarów jest obarczony błędem systematycznym (nie o rozkładzie normalnym z
wartością średnią 0). Czasami niektóre pomiary odbiegają od innych. Usunięcie ich z
analizy pozwala na przyjęcie hipotezy zerowej, czyli test �2 daje wynik pozytywny.
Pomiary pośrednie
Często, interesujące nas w pracy doświadczalnej wielkości, nie podlegają
bezpośredniemu pomiarowi. Np. w kinetyce chemicznej stała szybkości reakcji
opisywana jest zwykle równaniem Arrheniusa
Q
ńł �ł
K(T ) = k0 exp
�ł- żł
RT
ół �ł
gdzie: Q jest energią aktywacji procesu chemicznego i nie podlega bezpośredniemu
pomiarowi
K(T)  stała szybkości reakcji jako funkcja wyłącznie temperatury
(bezwzględnej) procesu,
R  uniwersalna stała gazowa
Badania kinetyczne prowadzi się zwykle dla szeregu dobranych warunków
izotermicznych (przy różnych wartościach temperatury bezwzględnej) i dla
wyznaczonych stałych szybkości reakcji poszukuje się stałych k0 i Q w równaniu
Arrheniusa. Dla ułatwienia wyznaczenia tych stałych, równanie to po obustronnym
zlogarytmowaniu uzyskuje liniową postać
Y = A + Bt
1
gdzie:
Y = ln(K(T )) ; A = ln k0 ; t = ; B = -Q / R
T
Uogólniamy zagadnienie pomiarów pośrednich w następujący sposób.
Załóżmy, że interesuje nas wyznaczenie zespołu r nieznanych (a priori) wielkości
(x1,...,xr), które nie podlegają bezpośredniemu pomiarowi. Wiemy zaś (lub
postulujemy), że istnieje liniowy (bądz
nieliniowy) związek między nimi a pewną
mierzalną zmienną �. Na początek, rozpatrywać będziemy związek liniowy:
�j = pj0 + pj1x1 +...+ pjr xr lub f = �j + aj0 + aj1x1 + ...+ ajr xr = 0
j
T
�ł
a =
Jeśli wektor wartości wsp. ajr zapisać w notacji wektora kolumnowego
j
�łaj1,aj 2,..., ajr łł
�ł
f = �j + aj0 + aT x = 0 j = 1,..., n
to powyższe relacje można zapisać jako:
j j
a po zdefiniowaniu macierzy i wektorów:
T T T
x = x1,..., xr ; � = �1,...,�n ; a0 = a01,...,a0n
[ ] [ ] [ ]
a11 a12 ... a1r a1 f1
�ł łł �ł łł �ł łł
�ła a22 ... a2r śł �ła śł �ł śł
f2
21 2
�ł śł �ł śł �ł śł
A = = f =
�ł śł �ł śł �ł śł
�" �" ń" �" �" �"
�ł �ł śł
ar ... arr śł �ł śł fn
�łar1 2 �ł �łar �ł �ł �ł
powyższy układ równań można zapisać (dla n dokonanych pomiarów):
f = � + a0 + Ax
Zakładamy, zgodnie z założeniami ogólnymi, że każdy pomiar zmiennej � o wartości
yj jest obarczony błędem �j o rozkładzie normalnym ze średnią zero:
yj = �j + � ; E � = 0 ; E �2 = �2 = 1/ g
( ) ( )
j j j j j
Zmienne yj są niezależne, stad związek między nimi jest wyrażony zerowymi
kowariancjami (cov(yi,yj)=0), a wtedy macierz kowariancji zmierzonych wartości jest
diagonalna, i także jej odwrotność zwana macierzą wag pomiarów:
2 g1 0 ... 0
�ł łł
�ł łł
�1 0 ... 0
�ł śł
�ł śł
-1 0 g2 ... 0
-1
0 �2 ... 0 �ł śł
2 G = Cy = = G� = C�
�ł śł ( )
( )
y
Cy = = C�
�ł śł
�" �" ń" �"
�ł śł
�" �" ń" �"
�ł śł
�ł śł
0 0 ... gn �ł
�ł
0 0 ... �2 �ł
�ł śł
�ł n
Wprowadzając wektor pomiarów i błędów:
T T
y = � + � ; y = y1,..., yn ; � = �1,...,�n
[ ] [ ]
Otrzymujemy układ równań
y - � + a0 + Ax = 0
Poszukujemy rozwiązania tego układu przy warunku, że suma kwadratów błędów �j
ważonych przez czynniki gj jest minimalna (jednocześnie to statystyka �2(n-r) dobroci
dopasowania model do danych):
n n
2
! = = =
"� / �2 "g �2 �TG �
j j j j y
j=1 j=1
Stosując tzw. metodę mnożników Lagrangea można uzyskać rozwiązanie tego
zagadnienia ze względu na x, otrzymując wektor estymat
-1
v
x = - ATGyA ATG (y + a0)
( )
y
z macierzą kowariancji wyznaczonych niewiadomych
-1
-1
Cx = Gx = ATGyA
( )
( )
v v
Pierwiastki kwadratowe z elementów przekątni głównej można utożsamiać z  błędami
pomiarowymi estymat zmiennych xi, mimo ze nie podlegały bezpośredniemu
pomiarowi.
Wektor estymat błędów pomiarowych przyjmuje wartość
-1
v v
� = Ax + y + a0 = -A ATG A ATG (y + a0) + (y + a0)
( )
y y
a wektor dopasowanych wartości pomiarów ( estymaty również)
-1
v v
� = y - � = A ATG A ATG (y + a0) - a0
( )
y y
-1 -1
-1
z macierzą kowariancji
C� = G� = A ATGyA AT = A Gx AT
( )
( ) ( )
v v v
W ramach pracy własnej, proszę sobie udowodnić, że gdy r=1 ( pomiary bezpośrednie
o indywidualnej dokładności �j2)
-1
n n n n n
�ł �ł �ł �ł
v
ATGyA = yj / �2 �ł =
�ł �ł �ł
"1/ �2 ; x = "1/ �2 " j "g yj /"g j
j j j
j=1 j=1 j=1 j=1 j=1
�ł łł �ł łł
Przykład. Dopasowanie linii prostej.
Niech zbiór wartości yj, pochodzący z pomiaru bezpośredniego, odpowiada zbiorowi
wartości tj (znanych dokładnie  tzw. zmienna kontrolowana t). Jeżeli ti, co się zdarza
nader często, jest również obarczone błędami,, to dopasowanie prostej staje się
problemem nieliniowych, osobno omawianym.
j t y �
1 0.0 1.4 0.5
�j = yj - � = x1 + x2t ; j = 1,...,n
j j
2 1.0 1.5 0.2
3 2.0 3.7 1.0
� - x1 + x2t = 0 ; a0 = 0 macierzowo
4 3.0 4.1 0.5
Obliczenia macierzowe dają następujące wyniki:
�ł-1 -t1 4 0 0 0
łł �ł-1 0
łł �ł łł
�ł śł �ł �ł0 25 0 0śł
�ł-1 -1 -1 -1 34 39 -1 0.094 -0.057
łł �ł łł �ł łł
�ł-1 -t2 śł �ł-1 -1śł �ł śł
śł
A = = AT = G = ATG A =
( )
y v y
�ł �ł39 65śł Cx = ATG A = �ł śł
�ł-1 -t3 śł �ł śł
śł �ł-1 -2 0 -1 -2 -3śł y 0 0 1 0
�ł �ł �ł �ł �ł-0.057 0.049 �ł
�ł śł �ł śł �ł śł
�ł-1 -t4 �ł �ł-1 -3�ł �ł0 0 0 4�ł
0.636 0.094 0.038 -0.019 -0.076 0.31
�ł łł �ł łł �ł łł
�ł1.702śł �ł śł �ł0.17śł
0.636 0.307 0.038 0.031 0.023 0.016
�ł łł �ł łł
�ł śł �ł śł �ł śł
v v v v v
x = �x = � = -Ax = C� = �� =
�ł1.066śł �ł0.222śł
�ł śł �ł-0.019 0.023 0.065 0.107 0.26
śł �ł śł
2.768
�ł �ł �ł �ł
�ł śł �ł śł �ł śł
�ł3.834�ł �ł-0.076 0.016 0.107 0.199 �ł �ł0.45�ł
Minimum funkcji kryterialnej ! = 4.507, dobroć dopasowania dla l.st.sw.=4-2=2
testujemy względem �2 (0.05,2)=5.99. Wniosek: nie ma podstaw do podważenia
kryt
dobrego dopasowania linii prostej do danych empirycznych (poziom ufności
dopasowania 95%, poziom istotności 5%).
Należy zwrócić uwagę, na fakt że błędy po dokonaniu dopasowania (błędy a posteriori)
są mniejsze niż błędy wyjściowe. Poniżej porównanie 2 linii regresji, czerwona dla
jednakowych błędów wielkości mierzonych ( błędów a priori), zielona dla
zróżnicowanych błędów.
Linia czerwona - regresja z jednakowymi błędami y=1.130+1.030*t
linia zielona - regresja z róznymi błędami y=0.636+1.066*t
5
4
3
2
1
0
0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5
t
Wykonanie obliczeń przedstawione w pliku: Dopasowanie_liniowe.xls
y