Dopasowanie funkcji do wynikw pomiarw.

Dopasowanie funkcji liniowej i wielomianu Dopasowanie funkcji liniowej

Zmienna

Y

jest liniową funkcją

X

X

A

A

Y

=

₀

+

₁

⋅

Wyniki pomiarów

Y

mają rozkład normalny o dyspersji

σ

_Y, która w

ogólnym przypadku również zależy od

X

. Wartości

X

nie są obarczone

błędami lub, w najgorszym wypadku, można je pominąć w porównaniu z błędami pomiaru wielkości

Y

. W wyniku wykonanych pomiarów

dysponujemy zestawem wartości

}

,

,

{

i y i i

y

x

σ

gdzie

i

=

1

,...,

n

, a i y

σ

są dyspersjami wartości mierzonych

y

_i.

Do znalezienia estymatorów

a

₀ i

a

₁ parametrów

A

₀ i

A

₁ możemy

wykorzystać Metodę Największej Wiarogodności. Odpowiednie prawdopodobieństwo będzie teraz wynosiło

(

)

∏

=

















₋

₋

−

=

n i y i i y y i i i i i

x

a

a

y

y

x

a

a

P

1 2 2 1 0 1 0

2

exp

2

1

})

,

,

{

;

,

(

σ

σ

π

σ

Maksymalizacja tego prawdopodobieństwa odpowiada minimalizacji funkcji

∑

=

















₋

₋

=

n i y i i y i i i i

x

a

a

y

y

x

a

a

1 2 1 0 1 0 2

₍

_,

_;

_{

_,

_,

_})

σ

σ

χ

ze względu na wartości

a

₀ i

a

₁.

Szukamy zatem rozwiązania układu równań

(

)

0

1

2

1 1 0 2 0 2

=

















−

−

−

=

∂

∂

_∑

= n i i i y

x

a

a

y

a

σ

_i

χ

(

)

0

2

1 1 0 2 1 2

=

















−

−

−

=

∂

∂

_∑

= n i i i y i

x

a

a

y

x

a

σ

_i

χ

Można je przekształcić w układ równań liniowych dla szukanych estymatorów

a

₀ i

a

₁

















=

















+

















∑

∑

∑

= = = n i y i n i y i n i yi i i

y

x

a

a

1 2 1 2 1 1 2 0

1

σ

σ

σ

















=

















+

















∑

∑

∑

= = = n i y i i n i y i n i y i i i i

y

x

x

a

x

a

1 2 1 2 2 1 1 2 0

_σ

_σ

_σ

Wprowadzamy następujące oznaczenia

2

1

i y i

w

=

σ

∑

=

i w

w

S

∑

=

i i x

w

x

S

∑

=

i i y

w

y

S

∑

=

2 i

x

w

S

_{xx i}

∑

=

i i i xy

w

x

y

S

i przy ich pomocy zapiszemy rozwiązanie układu

2 0 x xx w xy x y xx

S

S

S

S

S

S

S

a

−

−

=

2 1 x xx w y x xy w

S

S

S

S

S

S

S

a

−

−

=

Wariancje estymatorów 2 0 a

σ

i 2 1 a

σ

można ustalić podobnie jak poprzednio, korzystając z reguły propagacji wariancji

∑

=

























∂

∂

=

n i y i a i

y

a

1 2 2 0 2 0

σ

σ

∑

=

























∂

∂

=

n i y i a i

y

a

1 2 2 1 2 1

σ

σ

W tym celu należy znaleźć pochodne cząstkowe

a

₀ i

a

₁ względem wszystkich

y

_i. 2 0 x xx w x i i xx i i

S

S

S

S

x

w

S

w

y

a

−

−

=

∂

∂

2 1 x xx w x i w i i i

S

S

S

S

w

S

x

w

y

a

−

−

=

∂

∂

Po podstawieniu i wykonaniu przekształceń otrzymujemy

=

























∂

∂

=

∑

= n i y i a i

y

a

1 2 2 0 2 0

σ

σ

=

























−

−

∑

= n i w xx x i x i i xx i

w

S

S

S

S

x

w

S

w

1 2 2

1

(

)





=









−

+

−

∑

= n i i x i i xx x i i xx i x xx w

w

S

x

w

S

S

x

w

S

w

S

S

S

1 2 2 2 2 2 2 2 2

2

1

(

)

(

−

)

=

+

−

∑

= 2 2 1 2 2 2

₂

x xx w n i x i i xx x i i xx i

S

S

S

S

x

w

S

S

x

w

S

w

(

)

(

) (

)

(

−

)

=

+

−

∑

∑

∑

= = = 2 2 1 2 2 1 1 2

₂

x xx w n i x i i n i xx x i i n i xx i

S

S

S

S

x

w

S

S

x

w

S

w

(

−

)

=

+

−

2 2 2 2 2

₂

x xx w x xx xx x xx w

S

S

S

S

S

S

S

S

S

(

)

(

₂

)

2

(

2

)

2 2

2

x xx w xx x xx w x x xx w xx

S

S

S

S

S

S

S

S

S

S

S

S

−

=

−

+

−

=

























∂

∂

=

∑

= n i y i a i

y

a

1 2 2 1 2 1

σ

σ

=

























−

−

∑

= n i _{w xx x i} x i w i i

w

S

S

S

S

w

S

x

w

1 2 2

1

(

)



=











−

+

−

∑

= n i i x i x w i i w i i x xx w

w

S

w

S

S

x

w

S

x

w

S

S

S

1 2 2 2 2 2 2 2 2

2

1

(

)

(

−

)

=

+

−

∑

= 2 2 1 2 2 2

₂

x xx w n i x i x w i i w i i

S

S

S

S

w

S

S

x

w

S

x

w

(

)

(

) ( )

(

−

)

=

+

−

∑

∑

∑

= = = 2 2 1 2 1 1 2 2

₂

x xx w n i x i n i x w i i n i w i i

S

S

S

S

w

S

S

x

w

S

x

w

(

−

)

=

+

−

2 2 2 2 2

₂

x xx w x w x w xx w

S

S

S

S

S

S

S

S

S

(

)

(

₂

)

2

(

2

)

2 2

2

x xx w w x xx w x x xx w w

S

S

S

S

S

S

S

S

S

S

S

S

−

=

−

+

−

Czyli wariancje estymatorów

a

₀ i

a

₁ wynoszą

2 2 0 x xx w xx a

S

S

S

S

−

=

σ

2 2 1 x xx w w a

S

S

S

S

−

=

σ

Jeżeli parametry zależności

Y

od

X

są wielkościami, które mamy

wyznaczyć na podstawie wyników pomiarów

{

x

_i

,

y

_i

,

σ

_yi

}

, to możemy

zapisać: 2 0 x xx w xy x y xx

S

S

S

S

S

S

S

a

−

−

=

,

(

₀

)

₂ x xx w xx

S

S

S

S

a

u

−

=

2 1 x xx w y x xy w

S

S

S

S

S

S

S

a

−

−

=

,

(

₁

)

₂ x xx w w

S

S

S

S

a

u

−

=

Przykład. Rozkład Poisson'a - dopasowanie linii prostej

W pomiarach zależności intensywności promieniowania od odległości

r

od źródła rejestrowane szybkości zliczeń, o rozkładzie Poisson'a, zależą liniowo od odwrotności kwadratu odległości. Rejestracja przebiega w stałych odcinkach czasu i liczba zliczonych impulsów

y

_i ma rozkład

Poisson'a o średniej

ax

_i

+

b

, gdzie

x

_i

=

1

r

_i2 . Wynika z tego, że

dyspersja

i

-tego pomiaru wynosi _y

ax

_i

b

i

=

+

σ

, a waga

b

ax

w

i i

+

=

1

. Informację o zależności błędu (wagi) pomiaru od zmiennej niezależnej można wykorzystać na dwa sposoby. Można te wagi wstawić do

poprzednio ustalonego wyrażenia na

_χ

2_{, otrzymanego dla rozkładu}

normalnego (zakładając, że liczby zliczeń są na tyle duże, że rozkład Poisson’a jest bliski normalnemu),

(

)

∑

=

+

−

−

=

n i i i i

b

ax

b

ax

y

b

a

1 2 2

₍

_,

₎

χ

dostając po zróżniczkowaniu układ dwóch równań

(

)

0

1 2 2 1 2

=

+

−

=

∂

∂

_∑

_∑

= = n i _i i i n i i

b

ax

y

x

x

a

χ

(

)

0

1 2 2 2

=

+

−

=

∂

∂

∑

= n i _i i

b

ax

y

n

b

χ

które są nieliniowe ze względu na parametry

a

i

b, i rozwiązać je

Drugi sposób polega na zastąpieniu rozkładu normalnego rozkładem Poisson'a o wartości średniej

µ

_i

=

y

(

x

_i

)

=

a

⋅

x

_i

+

b

i ponownym

skonstruowaniu funkcji rozkładu gęstości prawdopodobieństwa

})

,

{

;

,

(

a

b

x

_i

y

_i

p

otrzymania zestawu wyników

{

x

_i

,

y

_i

}

pod warunkiem,

że parametry określające wartość średnią

µ

mają wartości

a

i

b

∏

= −









=

n i x y i y i i i

e

y

x

y

b

a

p

1 ) (

!

)]

(

[

)

,

(

i bezpośrednim zastosowaniu MNW. Zamiast maksymalizować prawdopodobieństwo można równie dobrze maksymalizować jego logarytm, co znacznie upraszcza wzory

=

−

−

=

∑

∑

∑

= = = n i i n i n i i i i

y

x

y

x

y

y

b

a

p

1 1 1

!

)

(

)]

(

ln

[

)

,

(

ln

.

)

(

)]

(

ln

[

1 1

const

x

y

x

y

y

n i n i i i i

−

+

∑

∑

= = Po zróżniczkowaniu względem

a

otrzymujemy

∑

∑

= =









∂

∂

−













∂

∂

=

∂

∂

n i n i i i i i

a

x

y

a

x

y

x

y

y

a

p

p

1 1

)

(

)

(

)

(

1

1

_.

Warunkiem zerowania się pochodnej

a

p

∂

∂

_{jest zerowanie się prawej}

strony. Biorąc pod uwagę, że i _i

x

a

x

y

₌

∂

∂

(

)

otrzymujemy pierwsze z układu równań

0

1 1

=

+

−

∑

∑

= = n i i i i n i i

b

ax

y

x

x

Różniczkowanie logarytmu

ln

p

(

a

,

b

)

względem

b

daje:

∑

∑

= =









∂

∂

−













∂

∂

=

∂

∂

n i n i i i i i

b

x

y

b

x

y

x

y

y

b

p

p

1 1

)

(

)

(

)

(

1

1

a uwzględniając, że

(

)

=

1

∂

∂

b

x

y

_{i otrzymujemy drugie równanie}

0

1

=

+

−

∑

= n i _i i

b

ax

y

n

Ostatecznie układ równań przyjmuje w tym przypadku postać

0

1 1

=

+

−

∑

∑

= = n i _i i i n i i

b

ax

y

x

x

0

1

=

+

−

∑

= n i i i

b

ax

y

n

bardzo podobną do poprzedniego układu.

Różnicy między tymi dwoma sposobami można się spodziewać dla mniejszych liczb zliczonych impulsów, kiedy różnice między

y

_i a

b

x

a

⋅

_i

+

stają się relatywnie większe. Wykorzystanie funkcji

χ

2,

uwarunkowane normalnym rozkładem gęstości prawdopodobieństwa, wymaga większych liczb impulsów lub grupowania pomiarów. Uniknięcie problemów związanych z grupowaniem wyników i tym samym

zmniejszeniem liczby punktów pomiarowych oraz ustalaniem właściwych wag jest możliwe przez bezpośrednie wykorzystanie Metody Największej Wiarogodności, jak w powyższym przykładzie.

Dopasowanie wielomianu

Jeżeli linia prosta nie daje się dobrze dopasować do danych

pomiarowych to można spróbować dopasować bardziej złożoną funkcję, to jest wielomian stopnia

m

∑

=

=

m k k k

x

a

x

y

0

)

(

lub w bardziej ogólnej formie

∑

=

=

m k k k

f

x

a

x

y

0

)

(

)

(

gdzie

f

_k

(x

)

są dowolnymi funkcjami nie zawierającymi parametrów

a

_k.

Jeżeli wyniki pomiarów

y

_i mają rozkład normalny o dyspersji i y

σ

, to maksymalizacja funkcji gęstości prawdopodobieństwa

∏

∑

∑

= = =



























 −

−

⋅

















=

n i n i m k i k k i y y y i i m

x

y

y

a

f

x

a

a

P

i i i 1 1 2 0 2 0

(

)

1

2

1

exp

2

1

})

,

,

{

;

,...,

(

σ

σ

π

σ

sprowadza się do minimalizacji funkcji

∑

∑

= =



























 −

=

n i m k i k k i y y i i m

x

y

y

a

f

x

a

a

i i 1 2 0 2 0 2

₍

_,...,

_;

_{

_,

_,

_})

1

₍

₎

σ

σ

χ

względem parametrów

a

_k.

Różniczkowanie

_χ

2 _{względem parametru}

l

a

daje

=



























 −

−

=

∂

∂

_∑

_∑

= = n i i l m k i k k i y l

x

f

x

f

a

y

a

i 1 2 0 2

)

(

)

(

1

2

σ

χ

[

]

=





























+

















−

∑

∑

∑

= = = n i m k i k i l k y n i _y i l i

a

f

x

f

x

x

f

y

i i 1 0 2 1 2

)

(

)

(

1

2

)

(

2

σ

σ

∑ ∑

∑

= = =

































+

















−

m k n i y i k i l k n i y i l i i i

x

f

x

f

a

x

f

y

0 1 2 1 2

)

(

)

(

2

)

(

2

σ

σ

W rezultacie otrzymujemy układ

m + 1

równań liniowych (dla wskaźnika

l

przebiegającego wartości od

0

do

m

).

∑ ∑

∑

= = =

































=

















m k n i y i k i l k n i y i l i i i

x

f

x

f

a

x

f

y

0 1 2 1 2

)

(

)

(

)

(

σ

σ

Układ równań można zastąpić równoważnym mu równaniem macierzowym

α

a

β

=

gdzie

a

i

β

są macierzami jednowymiarowymi o

m + 1

elementach

(wektorami), a

α

jest kwadratową i symetryczną macierzą

(m + 1)×(m + 1)

. Elementami macierzy

a

są szukane parametry dopasowania

a

_k, a elementy macierzy

β

i

α

wynoszą odpowiednio

∑

=

















=

n i y i k i k i

x

f

y

1 2

)

(

σ

β

∑

=

















=

=

n i y i k i l kl lk i

x

f

x

f

1 2

)

(

)

(

σ

α

α

Jeżeli macierz

α

jest nieosobliwa (ma wyznacznik różny od zera), to możemy obie strony równania macierzowego pomnożyć prawostronnie przez macierz

ε

odwrotną do macierzy

α

(

_ε

₌

_α

−1_,

_α

_α

−1

₌

₁

_).

Otrzymujemy wtedy

a

α

α

a

ε

α

a

ε

β

=

=

−1

=

co daje 1 −

=

=

β

ε

β

α

a

Rozwiązanie w formie macierzowej można przedstawić inaczej w formie skalarnej

∑

=

=

m k kl k l

a

0

ε

β

Macierz

α

nazywana jest macierzą krzywizny funkcji

_χ

2 _{w przestrzeni}

pochodnymi funkcji

_χ

2_. k l n i y i k i l lk

a

a

x

f

x

f

i

∂

∂

∂

=

















=

∑

= 2 2 1 2

2

1

)

(

)

(

χ

σ

α

Macierz

ε

jest nazywana macierzą (wariancji-)kowariancji. Macierz

ε

jest również symetryczna. Jej elementy diagonalne są równe wariancjom estymatorów

a ...

₀

a

_m, a pozostałe kowariancjom między odpowiednimi

parami estymatorów. jj aj

ε

σ

2

=

jl a a_{j l}

ε

σ

=

Jeżeli w szczególności dopasowywane funkcje

f

_k

(x

)

tworzą układ

ortogonalny, to w macierzach

α

i

ε

znikają elementy poza główną przekątną. Oznacza to brak korelacji między parametrami i ułatwia

rozszerzanie modelu opisanego równaniem

y

(

x

)

=

∑

a

_k

f

_k

(

x

)

o kolejne

składniki.

Uwaga o linearyzacji modelu zależności między wielkościami mierzonymi Przedstawione powyżej techniki dopasowania funkcji do wyników

pomiarów były ograniczone tylko do funkcji, które są liniowe względem szukanych parametrów. W przypadku kiedy funkcja nie jest liniowa względem swoich parametrów stosuje się często przekształcenie zmiennych w taki sposób, żeby funkcja opisująca zależność

przekształconych zmiennych była już liniowa względem parametrów. Jeżeli jednak przekształcana jest zmienna zależna (przekształceniem nieliniowym), to zniekształceniu ulegają różnice między pomiarami i modelem i rozwiązanie problemu dla modelu zlinearyzowanego nie odpowiada minimum

_χ

2 _{dla modelu oryginalnego.}

Przykłady rachunkowe

Pomiary spadku napięcia wzdłuż drutu oporowego z prądem

b ax U = + Nr pomiaru położenie cm różnica potencjału V niepewność pomiaru napięcia V xi2 xiUi dopasowane napięcie axi + b kwadraty różnic 1 10,0 0,37 0,05 100 3,70 0,33 0,001328 2 20,0 0,58 0,05 400 11,60 0,60 0,000247 3 30,0 0,83 0,05 900 24,90 0,86 0,000778 4 40,0 1,15 0,05 1600 46,00 1,12 0,000897 5 50,0 1,36 0,05 2500 68,00 1,38 0,000494 6 60,0 1,62 0,05 3600 97,20 1,64 0,000595 7 70,0 1,90 0,05 4900 133,00 1,91 0,000043 8 80,0 2,18 0,05 6400 174,40 2,17 0,000127 9 90,0 2,45 0,05 8100 220,50 2,43 0,000365 N = 9 Sx = 450,0 SU = 12,44 Sxx = 28500 SxU = 779,30 Σ=0,004874

( )

2 2

∑

∑

− = ∆ N xi xi = 54000

(

−

)

∆ =

∑

2

∑

∑

∑

/ i i i i i U x xU x b = 7,14·10-2 σU2

∑

xi2 /∆= 1,319·10 -3 _{u(b)= 3,63·10}-2 _V

(

−

)

∆ = N

∑

xiUi

∑

xi

∑

Ui / a = 262·10-4 Nσ_U2 /∆= 4,167·10-7 u(a)= 6,45·10-4 V/cm 0,00 0,50 1,00 1,50 2,00 2,50 3,00 0,0 20,0 40,0 60,0 80,0 100,0 po³o¿e nie , cm ró¿nica potencja³u, V

Liczby zliczeń impulsów licznika GM w jednakowych odcinkach czasu w funkcji odległości źródła cząstek od licznika b d a N = 1₂ + Nr pomiaru i odległość d (cm) x = 1/d 2 zliczenia N u(N) wi wixi wixi 2 wiNi wixiNi dopasowana liczba zliczeń axi + b kwadrat różnicy 1 20 25,00 901 30,0 0,00111 0,0277 0,6937 1,0 25,00 886,94 0,22 2 25 16,00 652 25,5 0,00153 0,0245 0,3926 1,0 16,00 610,66 2,62 3 30 11,11 443 21,0 0,00226 0,0251 0,2787 1,0 11,11 460,58 0,70 4 35 8,16 339 18,4 0,00295 0,0241 0,1966 1,0 8,16 370,09 2,85 5 40 6,25 283 16,8 0,00353 0,0221 0,1380 1,0 6,25 311,36 2,84 6 45 4,94 281 16,8 0,00356 0,0176 0,0868 1,0 4,94 271,09 0,35 7 50 4,00 240 15,5 0,00417 0,0167 0,0667 1,0 4,00 242,29 0,02 8 60 2,78 220 14,8 0,00455 0,0126 0,0351 1,0 2,78 204,77 1,05 9 75 1,78 180 13,4 0,00556 0,0099 0,0176 1,0 1,78 174,07 0,20 10 100 1,00 154 12,4 0,00649 0,0065 0,0065 1,0 1,00 150,19 0,09 n=10 Sumy Sw=0,03570 Sx=0,1868 Sxx=1,9122 Sy=10,0 Sxy=81,02 Ss=10,95 = − = ∆ _{( )}2 x xx wS S S 0,03339 = − =SySxx SxSxy b 119,50 u(b)= Sxx/∆ = 7,57 = − =SwSxy SxSy a 30,70 u(a)= S_w/∆ = 1,03

0 100 200 300 400 500 600 700 800 900 1000 0 20 40 60 80 100 120 o d le g ³o œæ, cm liczba zlicze ñ 0 100 200 300 400 500 600 700 800 900 1000 0,0 5,0 10,0 15,0 20,0 25,0 30,0

kwadrat odwrotnoœci odle g³oœci, 1/cm^2

liczba zlicze