Kilka słow o planowaniu eksperymentu

Kilka słów o planowaniu eksperymentu

Wojciech Myszka

25 maja 2010

Spis treści

1 Analiza regresji
1.1 Przypadek szczególny: dwa parametry
1.2 Przypadek ogólny
1.3 Nieliniowa funkcja regresji
2 Współczynnik korelacji wielowymiarowej
2.1 Istotność współczynnika korelacji
3 Rozkład prawdopodobieństwa współczynników regresji
4 Planowanie eksperymentu
4.1 Eksperyment dwupoziomowy
4.2 Plan ułamkowy
5 Plany eksperymentów ułamkowych

1 Analiza regresji

Niech nasz obiekt badań będzie opisany zależnością:

y = f(x ,x ,...,x ,z) 1 2 S

(1)

Mamy S wejść, jedno wyjście oraz dodatkową zmienną mającą wpływ na zachowanie obiektu, a mającą charakter „zaburzeń”.

Jeżeli przyjmiemy, że zaburzenia z mają charakter losowy nasza zależność z prostej zależności funkcyjnej przekształci się w zależność stochastyczną. Istotą tej zależności jest zależność rozkładu prawdopodobieństwa zmiennej losowej Y od wartości x przybranej przez zmienną losową X.¹

Będziemy mówić, że zmienna losowa Y jest zależna stochastycznie od zmiennej losowej X, jeśli dystrybuanta F zmiennej losowej przy warunku X = x jest nie tylko funkcją wartości y zmiennej losowej Y , lecz także wartości x

F (y|x) = P (Y < y|X = x) = f(y,x)

(2)

Budując opis matematyczny obiektu najbardziej interesuje nas wartość oczekiwana zmiennej losowej Y przy warunku X = x. Jeśli zachodzi:

E(Y |X = x) = f(x)

(3)

to mówimy o korelacji między zmiennymi losowymi Y i X.

Dokonujemy N pomiarów wartości x₁,x₂,…,x_S oraz odpowiadających im wartości y. Dane te możemy zgrupować w tablicy wejść i wyjść

x11, x12, ..., x1S, y1 x , x , ..., x , y 21 22 2S 2 ... ... ... ... ... xN1, xN2, ..., xNS, yN

(4)

Naszym celem jest wyznaczenie charakterystyki zastępczej o postaci

ˆy = f(x1,x2,...,xs;b0,b1,b2,...,bK )

(5)

zawierającej K + 1 nieznanych współczynników b_k;k = 0,1,2,…,K ≤ N

Dla n-tego pomiaru (w chwili n) funkcja przyjmuje wartość

ˆy = f(x ,x ,...,x ;b,b ,b ,...,b ) n n1 n2 nS 0 1 2 K

(6)

Charakterystyka 5 zwana jest zazwyczaj modelem matematycznym obiektu 1.

Współczynniki b_k,k = 0,1,…,K wyznaczamy w taki sposób aby wyznaczone z modelu wartości ŷ_n jak najmniej różniły się (były jak najbliżej) od wartości zmierzonych na obiekcie y_n. Najwygodniejsza jest odległość Euklidesa w przestrzeni N-wymiarowej.

N N Q = ∑ (y - ˆy )2 = ∑ [y - f(x ,x ,...,x ;b ,b ,b ,...,b )]2 n=1 n n n=1 n n1 n2 nS 0 1 2 K

(7)

Funkcja ŷ = f(x₁,x₂,…,x_s;b₀,b₁,b₂,…,b_K) gdzie współczynniki b_k wyznaczone zostały tak aby minimalizować 7 nazywana bywa funkcją regresji.

Bez większej straty ogólności możemy przyjąć, że funkcje regresji są liniowe względem parametrów

yˆ= ˆy(x;b ,b ,...,b ) = b + b f (x)+ ⋅⋅⋅+ b f (x) 0 1 K 0 1 1 k k

(8)

Przez x oznaczyłem [x₁,x₂,…,x_S]. O funkcjach f_k(x) zakłada się, że są znane i liniowo niezależne.

1.1 Przypadek szczególny: dwa parametry

Popatrzmy dokładniej na proces wyznaczania współczynników b_k. Załóżmy, że obiekt chcemy opisać zależnością liniową

y = b0 + b1x

(9)

Dokonaliśmy N(= 9) pomiarów wartości x i y, otrzymujemy zatem N równań o postaci

y = b + b x , n = 1,2,...,N n 0 1 n

(10)

o dwu niewiadomych b₀ i b₁.

Układu (w ogólnym przypadku) rozwiązać się nie da (jest sprzeczny?) chyba, że punkty (x_n,y_n) leżą na prostej.² Zatem, zamiast pierwotnego związku przyjmujemy związek zmodyfikowany:

yn = b0 + b1xn + en

(11)

(i mamy teraz N równań o N + 2 niewiadomych).

Kryterium Q dane równaniem 7 dba o to, by ∑ _n=1^Ne_n² → min!

Aby wyznaczyć b₀ i b₁ szukamy miejsc zerowych pochodnych cząstkowych ∂∂Qb0 i -∂Q
∂b1

∂Q N∑ ∂b--= - 2 (yn - b0 - b1xn) = 0 0 n=1

(12)

∂Q N∑ ∂b--= - 2 (yn - b0 - b1xn)xn = 0 1 n=1

(13)

Otrzymany układ (liniowych) równań normalnych rozwiązujemy otrzymując:

N N N N ∑ xnyn - ∑ xn ∑ yn --n=1-------n=1---n=1---- b1 = ∑N ( ∑N )2 N x2- xn n=1 n n=1 ( N N ) b = -1 ∑ y - b ∑ x 0 N n=1 n 1 n=1 n

(14)

Co ciekawsze ∑ _n=1^Ne_n = ∑ _n=1^N(y_n - b₀ - b₁x_n) = 0 Czemu?

1.2 Przypadek ogólny

Szukamy wielowymiarowej funkcji regresji w postaci:

yˆ= b + bx + b x + ⋅⋅⋅+ b x 0 1 1 2 2 S S

(15)

Obserwacje wejść i wyjść zapisane są w macierzy wejść X

⌊ ⌋ x10, x11, x12, ..., x1S || x20, x21, x22, ..., x2S || X = |⌈ ... ... ... ... ... |⌉ x , x , x , ... x N0 N1 N2 NS

(16)

W macierzy tej wprowadzono elementy „fikcyjne” x₁₀ = x₂₀ = ⋅⋅⋅ = x_N0 = 1

Wyniki obserwacji wyjść obiektu oraz wyjść modelu zapiszemy jako

′ y = [y1,y2,...,yN]

(17)

ˆy = [yˆ,ˆy ,...,ˆy ]′ 1 2 N

(18)

Nieznane parametry zapiszemy

b = [b0,b1,...,bS]′

(19)

W tym zapisie funkcja regresji przyjmuje postać

ˆy = Xb

(20)

Kryterium

∑N 2 ′ Q = (yn - ˆyn) = (y - ˆy) (y - ˆy ) n=1 = (y - Xb )′(y - Xb ) = y′y - y′Xb - b ′X ′y+ b ′X ′Xb ′ ′ ′ ′ ′ = y y - 2b X y + b X Xb → min!

(21)

Podobnie jak poprzednio wyznaczamy pochodną:

∂Q- = - 2X′y + 2X ′Xb ∂b

(22)

i otrzymujemy układ równań normalnych

X ′Xb = X ′y

(23)

który rozwiązujemy

b = (X ′X )-1X ′y

(24)

Łatwo sprawdzić (licząc drugą pochodną Q), że jest to istotnie minimum.

1.3 Nieliniowa funkcja regresji

Podobne rozważania można przeprowadzić w przypadku gdy

ˆy = b0f0(x)+ b1f1(x)+ ⋅⋅⋅+ bK fK(x)

(25)

gdzie f₀(X) = 1 jest funkcją fikcyjną.

Macierz pomiarów ma teraz postać

⌊ f0(x1), f1(x1), ..., fK (x1) ⌋ | f (x ), f (x ), ..., f (x ) | X = || 0 2 1 2 K 2 || ⌈ ..., ..., ..., ... ⌉ f0(xN ), f1(xN), ..., fK (xN )

(26)

gdzie x_n = [x_n1,x_n2,…,x_nS]; n = 1,2,…,N.

Pozostałe rozważania są identyczne jak w poprzednim przypadku.

W praktyce funkcje f mogą przyjmować postacie

f(x) = 1 f(x) = xs, s = 1,2,...,S 2 f(x) = x s, s = 1,2,...,S f(x) = xixj, i,j = 1,2,...,S, lecz i < j

(27)

co daje funkcję w postaci

∑S ∑S 2 ∑S ˆy = b0 + bsxs + bssx s + bijxixj s=1 s=1 i,j=1i<j

(28)

2 Współczynnik korelacji wielowymiarowej

Aby ocenić „jakość” (natężenie związku) związku między dwiema wielkościami (wyjściem obiektu y i wyjściem modelu ) wyznacza się współczynnik korelacji wielowymiarowej określony wzorem

N ∑ (y - Ży)(yˆ - Ży) n=1 n n R = R(y ˆy) = ┌│------------------------- │∘ N∑ 2∑N 2 (yn - Ży) (ˆyn - Ży) n=1 n=1

(29)

gdzie

∑N ∑N yŻ= 1- yn = ˆyn N n=1 n=1

(30)

Do praktycznych obliczeń używa się powyższego wzoru w wersji przekształconej, macierzowej (nie wymagającej znajomości wyjść modelu ŷ).

∘ ------------- b′X ′y - N Ży2 R = --′-------2- y y - N Ży

(31)

(Przekształcenie do tej postaci wymaga trochę rachunków i „zauważenia” paru rzeczy.)

2.1 Istotność współczynnika korelacji

Współczynnik korelacji R to liczba (0 ≤ R ≤ 1) zależna (patrz 31) od danych pomiarowych (wejść i wyjść obiektu), o których zakładaliśmy, że są zmiennymi losowymi. Zatem i R jest zmienną losową. Zatem ma jakiś rozkład. Pozwala to szacować czy wyznaczona wartość współczynnika korelacji jest bardzo prawdopodobne czy nie.

DO ZROBIENIA!

3 Rozkład prawdopodobieństwa współczynników regresji

Podobnie jak wcześniej zakładać będziemy, że obiekt opisany jest równaniem

y = X β + e

(32)

Przy czym teraz rozróżniamy prawdziwe wartości współczynników (β) od ich estymat (b). e to „wektor zaburzeń” normalnych o zerowej wartości średniej

E(e) = 0

(33)

i niezależnych składowych, o macierzy kowariancji

⌊ ⌋ σ2 0 ... 0 ′ || 0 σ2 ... 0 || cov(e) = E(ee ) = |⌈ ... ... ... ... |⌉ 0 0 ... σ2

(34)

Wartość oczekiwana zmiennej y

E(y ) = E ((X )β + e) = Xβ + E (e) = X β

(35)

Kowariancja y

{ } cov(y) = E [y- E((y)][y - E((y)]′ = E {[X β + e - X β][X β + e- Xβ ]′} = = E (ee ′) = Iσ2

(36)

Ocena parametrów β wyznaczona metodą najmniejszych kwadratów to zmienna losowa

b = (X ′X )-1X ′y

(37)

Wartość średnia b to

E(b ) = E [(X ′X)-1X ′y] = (X ′X)-1X ′E(y) = (X ′X )-1X ′X β = β

(38)

Kowariancja b

cov(b ) = E {[b - E (b][b - E (b]′} = E [(b - β)(b - β)′]

(39)

Wyliczmy b - β

b - β = (X ′X )-1X ′y- β = (X ′X)- 1X ′(Xβ + e) - β = (X′X )-1X ′e

(40)

Zatem

cov(b) = E {[(X ′X )-1X ′e][(X′X )-1X ′e]′} = E {(X′X )-1X ′ee′X [(X ′X )]′}

(41)

Wzór jest dosyć symetryczny, i można go przekształcać dalej

cov(b) = (X ′X )- 1X′E (ee ′)X (X′X )-1 = (X ′X)- 1X ′IX(X ′X)- 1σ2 = (X ′X)-1σ2

(42)

Macierz kowariancji mówi jakości wyznaczenia współczynników wektora b. na głównej przekątnej są wariancje składowych wektora, pozostałe wyrazy macierzy mówią jak bardzo jedna składową zależy od drugiej. Macierz jest symetryczna.

⌊ ⌋ ⌊ ⌋ var(b0) cov (b0b1) ... cov(b0bK) c00 c01 ... c0K || cov(b1b0) var(b1) ... cov(b1bK) || || c10 c11 ... c1K || cov(b) = |⌈ ... ... ... ... |⌉ = |⌈ ... ... ... ... |⌉σ2 cov(bK b0) ... ... var(bK) cK0 cK1 ... cKK

(43)

Współczynniki c_jk zwane są mnożnikami Gaussa.

var(bk) = σ2bk = ckkσ2 2 cov(bkbj) = ckjσ R (bkbj) = √-ckj--- ckkcjj

(44)

R(b_kb_j) to współczynnik korelacji między współczynnikami b_k i b_j.

To co jest najciekawszą cechą macierzy kowariancji cov(b) jest to, iż (w tym przypadku) zależy ona tylko od macierzy pomiarów. Oznacza to tyle, że dobierając odpowiednio wymuszenia przekazywane na obiekt możemy mieć wpływ na dokładność (i niezależność) wyznaczania parametrów funkcji regresji b.

Jeżeli uda się nam skonstruować tak macierz wymuszeń abu X′X była macierzą diagonalną, wówczas poszczególne składowe wektora b będą wyznaczane niezależnie (w sensie statystycznym) od siebie.

DOKOŃCZYĆ

4 Planowanie eksperymentu

4.1 Eksperyment dwupoziomowy

Bardzo często zdarza się, że zmienne wejściowe obiektu opisanego charakterystyką

y = f (x1,x2,...,xs,z)

(45)

mają ograniczony zakres

x0- Δx ≤ x ≤ x0- Δx , s = 1,2,...,S s s s s s

(46)

Dosyć naturalnym postępowaniem w takiej sytuacje jest przeprowadzenie eksperymentu polegającego na wygenerowaniu wszystkich możliwych 2^S kombinacji minimalnych i maksymalnych wartości zmiennych sterujących.

Jeżeli nasz model jest liniowy

S y = b0 + ∑ bsxs s=1

(47)

oraz jeżeli dokonamy zamiany zmiennych

xs - x0s ts = -Δx----, s = 1,2,...,S s

(48)

czyli

0 xs = xs + tsΔxs

(49)

to nasz model ulegnie przekształceniu

∑S ˆy = b0 + bs(x0s + tsΔxs ) = s=1 ∑S ∑S = b0 + bsx0s + bsΔxs s=1 s=1

(50)

Powyższe możemy zapisać jako

S∑ ˆy = k0 + ksts s=1

(51)

identyczny z pierwotnym model liniowy (o zmodyfikowanych współczynnikach).

Teraz jednak, zmienna t_s przyjmuje (przy eksperymencie dwupoziomowym) wartości +1 i -1.

Gdy S = 3 macierz wejść T (przez analogie do X) przyjmie postać

⌊ ⌋ 1 - 1 - 1 - 1 || 1 +1 - 1 - 1 || || 1 - 1 +1 - 1 || || 1 +1 +1 - 1 || T = || || || 1 - 1 - 1 +1 || || 1 +1 - 1 +1 || ⌈ 1 - 1 +1 +1 ⌉ 1 +1 +1 +1

(52)

Wyliczmy

⌊ ⌋ 8 0 0 0 || 0 8 0 0 || T′T = |⌈ 0 0 8 0 |⌉ = 8I 0 0 0 8

(53)

W przypadku ogólnym

N = 2S

(54)

oraz

1 T ′T = N I (T ′T )-1 = --I N

(55)

Natomiast współczynniki k wyliczamy z wzoru

′ -1 ′ 1 ′ k = (T T ) T Ty = N- T y

(56)

4.2 Plan ułamkowy

Gdy zmiennych wejściowych jest dużo, plan dwupoziomowy wymaga wykonania ogromnej liczby eksperymentów (gdy S = 10 wówczas 2¹⁰ = 1024, gdy S = 30 trzeba wykonać 1 Gi (Giga?) pomiarów. Zakładając, że każdy eksperyment (pomiar) zajmuje tylko 1s pomiary w przypadku eksperymentu całkowitego, dwupoziomowego zajmą ponad 30 lat.

Opracowano metody wykonywania eksperymentów ułamkowych zawierających tylko pewną liczbę doświadczeń wybranych z eksperymentu całkowitego.

Jeżeli przypatrzeć się dokładniej istocie elementu całkowitego, to spełnia on trzy następujące warunki:

Symetria doświadczeń względem środka eksperymentu:
$∑N tns = 0, s = 1,2,...,S. n=1$ (57)
Ortogonalność polegającą na zerowaniu się wszystkich iloczynów skalarnych wektorów kolumnowych macierzy eksperymentu T
$∑N tnitnj = 0, i,j = 1,2,...,S, i ⁄= j. n=1$ (58)
Równość sum kwadratów we wszystkich kolumnach macierzy eksperymentu T
$N∑ t2ns = N, s = 1,2,...,S n=1$ (59)

Najistotniejszą cechą jest ortogonalność, dzięki czemu macierz T′T jest diagonalna (co ogromnie ułatwia obliczenia i pozwala wyznaczać parametry modelu niezależnie).
Plany eksperymentów ułamkowych tworzy się tak aby spełniały warunki eksperymentu całkowitego, ale zawierały mniej punktów.

5 Plany eksperymentów ułamkowych

Plan całkowity eksperymentu dwupoziomowego dla trzech zmiennych (we współrzędnych standaryzowanych:

n	t₁	t₂	t₃

1	-	-	-
2	+	-	-
3	-	+	-
4	+	+	-
5	-	-	+
6	+	-	+
7	-	+	+
8	+	+	+

Mając do dyspozycji 8 pomiarów można wyznaczyć co najwyżej 8 parametrów funkcji regresji (czemu?). Jeżeli nasz model ma postać:

ˆy = k0 + k1ti + k2t2 + k3t3

(60)

wszystko jest w porządku.

Pełna funkcja regresji „wielomianowej” dla trzech zmiennych będzie zawierać:

wyraz zerowy — k₀ (1)
wszystkie zmienne t_i (3),
wszystkie iloczyny t_it_j (6),
wszystkie iloczyny postacji t_it_jt_k (9).

Zatem maksymalny wielomian bedzie zawierał 18 nieznanych parametrów. Żaden eksperymet dwupoziomowy nie pozwoli wyznaczyć wszystkich tych parametrów.

t₀

t₁

t₂

t₃

t₁²

t₂²

t₃²

t₁t₂

t₁t₃

t₂t₃

t₁t₂t₃

t₁²t₂

t₁³