Texty ke cviˇ cen´ım
Matematick´ a Statistika
Ivan Nagy, Pavla Pecherkov´ a, Jitka Homolov´ a
Obsah 1 Popisn´ a statistika
4
1.1
Sezn´amen´ı s programem MATLAB . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4
1.2
Popisn´e charakteristiky v MATLABu . . . . . . . . . . . . . . . . . . . . . . . . . . .
4
1.3
Pˇr´ıklady . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4
2 N´ ahodn´ y v´ ybˇ er
4
3 Pˇ r´ıklady
5
4 Bodov´ e odhady a jejich vlastnosti
5
4.1
Exponenci´aln´ı tˇr´ıda rozdˇelen´ı . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5
4.2
Nˇekter´a rozdˇelen´ı ve tvaru exponenci´aln´ı tˇr´ıdy . . . . . . . . . . . . . . . . . . . . .
5
4.3
Metoda maxim´aln´ı vˇerohodnosti . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
6
4.4
Metoda moment˚ u. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
6
4.5
Pˇr´ıklad (exponenci´aln´ı rozdˇelen´ı) . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
6
4.6
Pˇr´ıklad (alternativn´ı rozdˇelen´ı) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
7
4.7
Vlastnosti bodov´ ych odhad˚ u. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
7
4.8
Pˇr´ıklad (nestrannost a konzistence pro Alt(π)) . . . . . . . . . . . . . . . . . . . . .
8
4.9
Pˇr´ıklad (vydatnost)
8
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5 Pˇ r´ıklady
8
6 Intervaly spolehlivosti
9
6.1
Funkce MATLAB . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
9
6.2
Pˇr´ıklady . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
9
7 Testy parametrick´ ych hypot´ ez
9
7.1
Funkce MATLAB . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
9
7.2
Pˇr´ıklady . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
9
8 Pˇ r´ıklady
9
9 Chi2 testy hypot´ ez
10
9.1
Funkce MATLAB . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
9.2
Pˇr´ıklady . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
10 Dalˇ s´ı neparametrick´ e testy
10
10.1 Funkce MATLAB . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 10.2 Pˇr´ıklady . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 2
11 Regresn´ı anal´ yza
10
11.1 Funkce MATLAB . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 11.2 Pˇr´ıklady . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 12 Korelaˇ cn´ı anal´ yza
11
12.1 Funkce MATLAB . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 12.2 Pˇr´ıklady . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 13 Anal´ yza rozptylu (ANOVA)
11
13.1 Funkce MATLAB . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 13.2 Pˇr´ıklady . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 14 Anal´ yza hlavn´ıch komponent
11
14.1 Pˇr´ıklady . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 14.2 Funkce MATLAB . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 15 Pˇ r´ıklady ˇ reˇ sen´ e pomoc´ı MATLAB
12
15.1 Popisn´a statistika . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 15.2 Intervaly spolehlivosti . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 15.3 Testy hypot´ez . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 15.4 Chi2 testy hypot´ez . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 15.5 Dalˇs´ı neparametrick´e testy hypot´ez . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 15.6 Regresn´ı anal´ yza . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 15.7 Korelaˇcn´ı anal´ yza . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 15.8 ANOVA . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 15.9 Anal´ yza hlavn´ıch komponent . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43
3
Cviˇ cen´ı 1 1
Popisn´ a statistika
1.1
Sezn´ amen´ı s programem MATLAB
1.2
Popisn´ e charakteristiky v MATLABu
bas− stat bas− stat− 2 bas− stat− sort mean std var cov cor mode median iqr moment meansq sumsq sorted unsorted table
1.3
popisn´a statistika pro netˇr´ıdˇen´a data popisn´a statistika -”- dva soubory popisn´a statistika pro tˇr´ıdˇen´a data pr˚ umˇer smˇerodatn´a odchylka v´ ybˇerov´ y rozptyl v´ ybˇerov´a kovariance korelaˇcn´ı koeficient modus median mezikvartilov´e rozpˇet´ı v´ ybˇerov´e momenty pr˚ umˇer ˇctverc˚ u souˇcet ˇctverc˚ u tˇr´ıdˇen´a data netˇr´ıdˇen´a data kontingenˇcn´ı tabulka
Pˇ r´ıklady
• Udˇelat pr˚ umˇer a porovnat s mean. A dalˇs´ı. • Sestavit nˇeco jako bas− stat, ale podle sv´eho. • Generovat kostku a ovˇeˇrit 1/6 pro jedno ˇc´ıslo. A dalˇs´ı (sud´e,...). – Kostka: fix(6rand+1). – V´ıce hod˚ u: k=fix(6rand(1,1000)+1), udˇelat min, max, hist. – Souˇcet na dvou kostk´ach: k2=fix(6rand(1,1000)+1)+fix(6rand(1,1000)+1). MATLAB: p011, p012
2
N´ ahodn´ y v´ ybˇ er • Popsat v´ ybˇer na pˇr´ıkladech. Zd˚ uraznit n´ahodnost charakteristik v´ ybˇeru. V´ ybˇerov´ y pr˚ umˇer a jeho charakteristiky. • Rozdˇelen´ı v´ ybˇerov´ ych charakteristik. • Normov´an´ı a odnormov´an´ı v´ ybˇerov´eho pr˚ umˇeru. Udˇelat funkci. 4
3
Pˇ r´ıklady
cv01 a cv02.
Cviˇ cen´ı 2 4
Bodov´ e odhady a jejich vlastnosti • Vlastnosti bodov´ ych odhad˚ u na grafu hp statistiky. • Metoda moment˚ u. • Metoda maxim´aln´ı vˇerohodnosti. • Pˇ r´ıklady: – pro exponenci´aln´ı rozdˇelen´ı, – pro alternativn´ı rozdˇelen´ı, – MATLAB: pr˚ umˇer z pr˚ umˇer˚ u se bl´ıˇz´ı m´ı.
4.1
Exponenci´ aln´ı tˇ r´ıda rozdˇ elen´ı
Hp f (x, θ) je exponenci´aln´ı tˇr´ıdy jestliˇze f (x, θ) = exp{Q(θ)U (x) + R(θ) + V (x)} mnoˇzina {x|f (x, θ) 6= 0} nez´avis´ı na θ.
4.2
Nˇ ekter´ a rozdˇ elen´ı ve tvaru exponenci´ aln´ı tˇ r´ıdy
Alternativn´ı x
1−x
f (x, π) = π (1 − π) Poissonovo
f (x, λ) = e−λ
π = exp ln x + ln(1 − π) 1−π
λx = exp {ln(λ)x − λ − ln(x!)} x!
Norm´aln´ı (µ i σ 2 nezn´am´e) (
1 1 f (x, [µ, σ ]) = √ exp − 2 2πσ 2
(
= exp
[− 2σ1 2 ,
µ ] σ2
"
x2 x
#
1 − 2
x−µ σ
2 )
= )
!
µ2 1 + ln(σ 2 ) − ln(2π) σ2 2
Exponenci´aln´ı 1 x f (x, δ) = exp − δ δ
5
1 = exp − x − ln(δ) δ
4.3
Metoda maxim´ aln´ı vˇ erohodnosti
Exponenci´aln´ı tˇr´ıda f (x, θ) = exp{Q(θ)U (x) + R(θ) + V (x)} Podm´ınka extr´emu Q0 S + nR0 = 0
kde
S=
n X
U (xi )
i=1
Ovˇeˇren´ı maxima Q00 S + nR00 < 0
4.4
Metoda moment˚ u
Porovn´ame obecn´e momenty souboru E[X k ] =
Z
xk f (x, θ)dx
θ
a v´ ybˇeru
n 1X Mk = xk , n i=1
pro k = 1, 2, . . . podle toho, kolik parametr˚ u odhadujeme. Vˇetˇsinou odhadujeme jeden parametr, potom pouˇzijeme prv´e momenty, tj. E[X] = x.
4.5
Pˇ r´ıklad (exponenci´ aln´ı rozdˇ elen´ı)
Zkonstruujte odhadovou statistiku pro odhad parametru δ exponenci´aln´ıho rozdˇelen´ı. ˇ e ˇs e n´ı: R Maxim´ aln´ı vˇ erohodnost Exponenci´aln´ı tˇr´ıda 1 f (x, δ) = exp − x − ln(δ) δ
Odtud je U =x Q = − 1δ R = − ln(δ) Podm´ınka
S = ni=1 xi Q0 = δ12 R0 = − 1δ
P
Q00 = − δ23 R00 = δ12 n 1 X 1 xi − n = 0 2 δ i=1 δ
Ovˇeˇren´ı −
n 2 X 1 xi + n 2 < 0 3 δ i=1 δ
Momenty E[X] = δ,
n 1X xi = x, n i=1
Pozn´ amka: Udˇelat sp´ıˇse pro a = 1/δ. 6
δˆ = x.
⇒
⇒
1 < 1. 2
⇒
δˆ = x
4.6
Pˇ r´ıklad (alternativn´ı rozdˇ elen´ı)
Proved’te konstrukci odhadu parametru π n´ahodn´e veliˇciny s alternativn´ım rozdˇelen´ım. ˇ e ˇs e n´ı: R Maxim´ aln´ı vˇ erohodnost: Exponenci´aln´ı tˇr´ıda π f (x, θ) = exp ln x + ln(1 − π) 1−π
Odtud je U =x π Q = ln 1−π R = ln(1 − π)
Pn
S= Q0 = R0 =
Podm´ınka
i=1 xi 1 π(1−π) −1 1−π
Q00 = R00 =
2π−1 π 2 (1−π)2 −1 (1−π)2
n X 1 1 −n =0 π(1 − π) i=1 1−π
P
Ovˇeˇren´ı (dosad´ıme
⇒
π ˆ=x=p
xi = nπ) 2π − 1 1 nπ − n <0 π 2 (1 − π)2 (1 − π)2
⇒
π < 1.
Momenty E[X] =
1 X
xπ x (1 − π)1−x = 0π 0 (1 − π)1 + 1π 1 (1 − π)0 = π,
i=0 n 1X xi = x = p n i=1
Odtud porovn´an´ım
4.7
π ˆ = p.
Vlastnosti bodov´ ych odhad˚ u
Obecnˇe oznaˇc´ıme: T je odhadov´a statistika, θ je odhadovan´ y parametr. Nestrannost E[T ] = θ Konzistence lim P (|Tn − θ| > ) = 0
n→∞
Kriterium: - asymptotick´a nestrannost E[Tn ] = θ, pro n → ∞, - limita rozptylu jde k nule
D(Tn ) = 0, pro n → ∞.
Vydatnost Vyjadˇruje pˇresnost bodov´eho odhadu. Pro nestrann´e odhady se vydatnost mˇeˇr´ı rozptylem
D[T ].
Pro odhady, kter´e nejsou nestrann´e je mˇeˇr´ıtkem vydatnosti M SE M SE = E[(T − θ)2 ] = D[T ] + B(θ)2 . 7
4.8
Pˇ r´ıklad (nestrannost a konzistence pro Alt(π))
Je d´an v´ ybˇer o rozsahu n z n´ahodn´e veliˇciny X s pravdˇepodobnostn´ı funkc´ı f (x) = π x (1−π)1−x , x = P 0, 1. Ukaˇzte, ˇze statistika T = p = n1 ni=1 Xi je nestrann´ ym a konzistentn´ım odhadem parametru π. ˇ e ˇs e n´ı: R D[X] = π(1 − π).
E[X] = π, Nestrannost
n n 1X 1X E[T ] = E Xi = Xi = π. n i=1 n i=1
"
#
Konzistence lim P (|p − π| > ) ≤ lim
n→∞
Cviˇcen´ı: tot´eˇz, ale pro T = n+ =
4.9
n→∞
P
D[X] π(1 − π) = lim =0 n→∞ n2 n2
Xi . (Je nestran´ y, nen´ı konzistentn´ı).
Pˇ r´ıklad (vydatnost)
Z n´ahodn´e veliˇciny byl uˇcinˇen v´ ybˇer o rozsahu 3000 a vypoˇcteny pr˚ umˇery xi z prvn´ı poloviny dat, x2 ze dvou tˇretin druh´e poloviny a x3 ze zbytku. Pro odhad stˇredn´ı hodnoty byly pouˇzity dvˇe statistiky X1 + X2 X2 + X3 T1 = a T2 = . 2 2 a) Ovˇeˇrte nestrannost obou statistik. b) Zjistˇete, kter´a ze statistik ve vydatnˇejˇs´ı. ˇ e ˇs e n´ı: R a) E[T1 ] =
µ1 + µ2 = µ, 2
E[T2 ] =
µ2 + µ3 = µ. 2
Jsou nestrann´e b) 1 1 D[T1 ] = D X 1 + X 2 2 = (D[X 1 ] + D[X 2 ]) = 4 4
σ2 σ2 + 1500 1000
1 1 D[T2 ] = D X 2 + X 3 2 = (D[X 2 ] + D[X 3 ]) = 4 4
σ2 σ2 + 1000 500
h
h
i
i
T1 je vydatnˇejˇs´ı. Pozn´amka: Rozptyl souˇctu je zde souˇcet rozptyl˚ u. Proˇc?
5
Pˇ r´ıklady
cv03 a cv04
8
!
=
!
=
σ2 . 2400
σ2 . 1333.3
Cviˇ cen´ı 3 6 6.1
Intervaly spolehlivosti Funkce MATLAB
z− int t− int t− int− 2s t− int− 2n t− int− 2p prop− int prop− int− 2 var− int
6.2
interval interval interval interval interval interval interval interval
pro pro pro pro pro pro pro pro
stˇredn´ı hodnotu (zn´am´ y rozptyl) stˇredn´ı hodnotu (nezn´am´ y rozptyl) dvˇe stˇredn´ı hodnoty (sdruˇzen´ y) dvˇe stˇredn´ı hodnoty (nesdruˇzen´ y) dvˇe stˇredn´ı hodnoty (p´arov´e v´ ybˇery) pod´ıl dva pod´ıly rozptyl
Pˇ r´ıklady
p041, p042, p043
7 7.1
Testy parametrick´ ych hypot´ ez Funkce MATLAB
z− test z− test− 2 t− test t− test− 2s t− test− 2n t− test− 2p prop− test prop− test− 2 var− test var− test− 2
7.2
test test test test test test test test test test
stˇredn´ı hodnoty (zn´am´e sig) dvou stˇredn´ıch hodnot (zn´am´e sig) stˇredn´ı hodnoty (nezn´am´e sig) dvou stˇredn´ıch hodnot (sdruˇzen´ y) dvou stˇredn´ıch hodnot (nesdruˇzen´ y) dvou stˇredn´ıch hodnot (p´arov´ y) pod´ılu dvou pod´ıl˚ u rozptylu dvou rozptyl˚ u
Pˇ r´ıklady
p051, p052, p053, p054, p055, p056, p057
8
Pˇ r´ıklady
cv05, cv06 a cv07
9
Cviˇ cen´ı 4 9 9.1
Chi2 testy hypot´ ez Funkce MATLAB
chisquare− test chisquare− test− homogeneity chisquare− test− independence
9.2
vzd´alenost O a E dobr´a shoda nez´avislost
Pˇ r´ıklady
p061, p062, p063
10 10.1
Dalˇ s´ı neparametrick´ e testy Funkce MATLAB
sign− test wztest cor− test kolmogorov− smirnov− test
10.2
test test test test
dvou medi´an˚ u nez´avislosti prvk˚ u v´ ybˇeru nez´avislosti v´ ybˇer˚ u typu rozdˇelen´ı
Pˇ r´ıklady
p071, p072, p073, p074
Cviˇ cen´ı 5 11 11.1
Regresn´ı anal´ yza Funkce MATLAB
reg− desc reg− info lin− reg exp− reg pol− reg
popisn´a lin. regrese lin. regrese (soubor charakteristik) line´arn´ı regrese exponenci´aln´ı regrese polynomick´a regrese
10
11.2
12 12.1
Pˇ r´ıklady
Korelaˇ cn´ı anal´ yza Funkce MATLAB
reg− infe t− test− reg f− test− reg
12.2
inferenˇcn´ı lin. regrese t-test lin. regrese f-test lin. regrese
Pˇ r´ıklady
Cviˇ cen´ı 6 13 13.1
Anal´ yza rozptylu (ANOVA) Funkce MATLAB
anova anova− 2
13.2
anal´ yza rozptylu - jednoduch´a anal´ yza rozptylu - dvojn´a
Pˇ r´ıklady
Cviˇ cen´ı 7 14 14.1
Anal´ yza hlavn´ıch komponent Pˇ r´ıklady
p101, p102
14.2
Funkce MATLAB
pca− eig pca− svd
rozklad kovarianˇcn´ı matice rozklad datov´e matice
11
15 15.1
Pˇ r´ıklady ˇ reˇ sen´ e pomoc´ı MATLAB Popisn´ a statistika
****************** p011 ****************** % P011 - Popisn´ a statistika % ------------------------% - char. polohy a variability pro prosta data % - char. polohy a variability pro tˇ r´ ıdˇ en´ a data intro % ZAD´ AN´ I Pˇ R´ IKLADU % Mˇ eˇ ren´ ım jsme z´ ıskali dva n´ ahodn´ e v´ ybˇ ery "x" a "y" a % uloˇ zili je do datov´ eho souboru "data01". % % Urˇ cete: % a) stˇ redn´ ı hodnotu, rozptyl, smˇ erodatnou odchylku, modus, % medi´ an, variaˇ cn´ ı rozpˇ et´ ı, a mezikvartilov´ e rozpˇ et´ ı pro "x" % b) kovarianci a korelaˇ cn´ ı koeficient pro "x" a "y". % % Soubor "x" roztˇ ridte do pˇ eti interval˚ u a % c) urˇ cete v´ aˇ zen´ y pr˚ umˇ er a rozptyl, % d) urˇ cete modus a medi´ an tˇ r´ ıdˇ en´ ych dat. % ˇ REˇ SEN´ I load data01; disp(’Priklad P011’); disp(’ a) prosta data’); mx = mean(x); vx = var(x); sd = std(x); mo = mode(x,200); me = median(x); rr = max(x)-min(x); kr = iqr(x); fprintf(’stredni hodnota: %g\n’,mx) fprintf(’rozptyl: %g\n’,vx) fprintf(’smerodatna odchylka: %g\n’,sd) fprintf(’modus: %g\n’,mo) fprintf(’median: %g\n’,me) fprintf(’variacni rozpeti: %g\n’,rr) fprintf(’mezikvartilove rozpeti: %g\n’,kr) disp(’ b)’); cxy = cov(x,y); rxy = cor(x,y); fprintf(’kovariance: %g\n’,cxy) fprintf(’korelacni koeficient: %g\n’,rxy) 12
disp(’ c) tridena data’); [ff xx]=hist(x,5); xn.d=xx’; xn.n=ff’; mxs = mean(xn); vxs = var(xn); fprintf(’stredni hodnota: %g\n’,mxs) fprintf(’rozptyl: %g\n’,vxs)
disp(’ d)’); [a k]=max(xn.n); mos = xn.d(k); mes = median(unsorted(xn)); fprintf(’modus: %g\n’,mos) fprintf(’median: %g\n’,mes) ****************** p012 ****************** % P012 - Popisn´ a statistika % ------------------------% - tr´ ıdeni dat % - char. tridenych a netridenych dat intro % ZAD´ AN´ I Pˇ R´ IKLADU % Mˇ eˇ ren´ ım jsme z´ ıskali dva n´ ahodn´ e v´ ybˇ ery "x" a "y" a % uloˇ zili je do datov´ eho souboru "data01". % Soubor "x" roztˇ ridte interval˚ u s hranicemi h = [-2 0 2 6 10]; % a % a) urˇ cete v´ aˇ zen´ y pr˚ umˇ er a rozptyl takto tˇ r´ ıdˇ en´ ych dat, % b) porovnejte s pr˚ umˇ erem a rozptylem stejn´ ych dat, % tˇ r´ ıdˇ en´ ych rovnomˇ ernˇ e do pˇ eti skupin. % ˇ REˇ SEN´ I load data01; disp(’Priklad P012’); disp(’ a) trideni do predapsanych intervalu’); intervaly=h f = zeros(5,1); g = f; for j=1:4 f(j)=sum(x>h(j) & x<=h(j+1)); g(j)=(h(j+1)+h(j))/2; end%for f(1)=f(1)+1; xn.d=g; xn.n=f; 13
mx1 = mean(xn) vx1 = var(xn) disp(’ b) trideni do 5 ekvidistantnich intervalu’); [f g]=hist(x,5); intervaly=g xn.d=g; xn.n=f; mx2 = mean(xn) vx2 = var(xn)
15.2
Intervaly spolehlivosti
****************** p041 ****************** % P041 - Intervaly spolehlivosti (pravdˇ epodobnost ˇ ze ..) % -----------------------------------------------------intro % % % % % % % % % %
ZAD´ AN´ I Pˇ R´ IKLADU Pˇ redpokl´ ad´ ame, ˇ ze velk´ y roˇ cn´ ık va Vˇ S m´ a v´ yskedky testu rozdˇ eleny norm´ alnˇ e, se stˇ redn´ ı hodnotou 72 bod˚ u a smˇ erodatnou odchylkou 9 bod˚ u. a) Urˇ cete pravˇ epodobnost, ˇ ze n´ ahodnˇ e vybran´ y student bude m´ ıt v´ ysledek nad 80 bod˚ u. b) Provedeme v´ ybˇ er o rozsahu n=4 a n=9. Jak´ a je pravdˇ epodobnost, ˇ ze v´ ybˇ erov´ y pr˚ umˇ er z tˇ ehto v´ ybˇ er˚ u bude vˇ etˇ s´ ı neˇ z 80?
% ˇ REˇ SEN´ I % zadan´ e hodnoty m = 72; s = 9; h = 80; n = 4;
% % % %
% a) z=(h-m)/s; pv1=1-normal_cdf(z);
% normov´ an´ ı pro veliˇ cinu % pravdˇ epodobnost
% b) z=(h-m)/s*sqrt(n); pv2=1-normal_cdf(z);
% normov´ an´ ı pro v´ ybˇ er % pravdˇ epodobnost
stˇ redn´ ı hodnota smˇ erodatn´ a odchylka hranice bod˚ u velikost v´ ybˇ eru
fprintf(’Test velkeho rocniku na VS (1)\n\n’); fprintf(’Pr. nad %g pro studenta je: %g\n’,h,pv1); fprintf(’Pr. nad %g pro vyber n=%g je: %g\n’,h,n,pv2);
14
% % % %
POZN´ AMKY: ˇ Reˇ seni b) je pro n=4. Pro n=9 se hodnota n pˇ rep´ ıˇ se. Tak je moˇ zno mˇ enit i dakˇ s´ ı hodnoty, pˇ r´ ıpadnˇ e upravit cel´ y m-file.
% Funkce ’normal_cdf(x)’ je tot´ ez ˇ jako % ’normal_cdf(x,0,1)’ coˇ z je ’stdnormal_cdf(x)’ % => std na zaˇ c´ atku funkce je tedy moˇ zno vynechat! ****************** p042 ****************** % P042 - Intervaly spolehlivosti (v´ ypoˇ cet hranice intervalu) % --------------------------------------------------------intro % % % % % % % % % %
ZAD´ AN´ I Pˇ R´ IKLADU Velk´ y roˇ cn´ ık na Vˇ S psal test. Ze zkuˇ senosti v´ ıme, ˇ ze v´ ysledky testu jsou rozdˇ eleny norm´ alnˇ e, se smˇ erodatnou odchylkou 9 bod˚ u. Stˇ redn´ ı hodnotu, danou konkr´ etn´ ı obt´ ıˇ znost´ ı testu, nezn´ ame. Provedli jsme v´ ybˇ er x = [64, 66, 89, 77]; a) Urˇ cete interval, ve kter´ em leˇ z´ ı stˇ redn´ ı hodnota m´ ı s pravdˇ epodobnost´ ı 0.95 (nebo, ve kter´ em leˇ z´ ı 95% vˇ sech odhad˚ u stˇ redn´ ı hodnoty - coˇ z je tot´ eˇ z). b) Tot´ eˇ z, ale pro nezn´ am´ y rozptyl. c) Tot´ eˇ z, ale pro rozptyl.
% ˇ REˇ SEN´ I % zadan´ e hodnoty sig = 9; al = .05;
% smˇ erodatn´ a odchylka % pravdˇ epodobnost
% a) n=length(x); mx=mean(x); z_kr=-normal_inv(al/2); del=sig/sqrt(n)*z_kr; is1=[mx-del, mx+del];
% % % % %
rozsah v´ ybˇ eru pr˚ umˇ er x krit. hodnota polomˇ er IS IS
% b) s=std(x); t_kr=-t_inv(al/2,n-1); del=s/sqrt(n)*t_kr; is2=[mx-del, mx+del];
% % % %
smˇ erodatn´ a odchylka krit. hodnota polomˇ er IS IS
% c) ch_krL=chisquare_inv(1-al/2,n-1); % kr. hod. lev´ a ch_krR=chisquare_inv(al/2,n-1); % kr. hod. prav´ a s2=var(x); % rozptyl 15
is3(1)=(n-1)*s2/ch_krL; is3(2)=(n-1)*s2/ch_krR;
% IS lev´ a hranice % IS prav´ a hranice
fprintf(’Test velkeho rocniku na VS (2)\n\n’); fprintf(’IS pro mi (sig zname) je: (%g, %g)\n’,is1); fprintf(’IS pro mi (sig nezname) je: (%g, %g)\n’,is2); fprintf(’IS pro sig2 je: (%g, %g)\n’,is3); % % % % %
POZN´ AMKY: lze pouˇ z´ ıt funkce pro IS z_int(x,sig**2,.95); t_int(x,.95); var_int(x,.95);
****************** p043 ****************** % P043 - Intervaly spolehlivosti (pro m´ ı, rozsah IS) % -------------------------------------------------intro % % % % % % % %
ZAD´ AN´ I Pˇ R´ IKLADU N´ ahodn´ a veliˇ cina X m´ a norm´ aln´ ı rozdˇ elen´ ı se stˇ redn´ ı hodnotou m´ ı a rozptylem sig2. Provedli jsme n´ ahodn´ y v´ ybˇ er x = [25, 17, -5, -4, 10, 22, 11, 2, 5, 12]; Urˇ cete a) 95% IS pro stˇ redn´ ı hodnotu, b) maxim´ aln´ ı chybu oboustrann´ eho intervalu, c) rozsah v´ ybˇ eru pro pro maxim´ aln´ ı chybu 0.2.
% ˇ REˇ SEN´ I % zadan´ e hodnoty al = .05; del0 = .2; % a) is1=t_int(x, al, ’<>’); is2=t_int(x, al, ’>’); is3=t_int(x, al, ’<’);
% IS obou.. % IS pravo.. % IS levo..
% b) n=length(x); s=std(x); t_kr=-t_inv(al/2,n-1); del=s/sqrt(n)*t_kr;
% % % %
% c) nn=s/del0*t_kr; n0=fix(nn)+1;
% rozsah v´ ybˇ eru (zlomek) % -’(zaokrouhleno)
rozsah v´ ybˇ eru smˇ erodatn´ a odchylka kritick´ a hodnota polomˇ er intervalu
16
disp(’ a)’); fprintf(’IS pro mi a rozsah vyberu\n\n’); fprintf(’IS obou.. je: (%g, %g)\n’,is1); fprintf(’IS pravo.. je: (%g, %g)\n’,is2); fprintf(’IS levo.. je: (%g, %g)\n’,is3); disp(’ b)’); fprintf(’maximalni chyba: %g\n’,del); disp(’ c)’); fprintf(’rozsah vyberu: %g\n’,n0); % POZN´ AMKY: ****************** p044 ****************** % PO44 - Interval spolehlivosti pro rozd´ ıl dvou m´ ı % -----------------------------------------------intro % % % % %
% % % %
´N´ ZADA I Pˇ R´ IKLADU Sledovali jsme ´ uˇ cinek dvou protikorozn´ ıch l´ atek. Prvn´ ı jsme pouˇ zili ve 20 pˇ r´ ıpadech, druhou ve 25 pˇ r´ ıpadech. Po stanoven´ e dobˇ e jsme zjistili stupeˇ n poˇ skozen´ ı s v´ ysledky mx1=82.4; s2x1=12; mx2=80; s2x2=10; Urˇ cete 95% IS pro shodu stˇ redn´ ıch hodnot obou test˚ u za pˇ redpokladu normality obou rozdˇ elen´ ı. Rozptyly povaˇ zujte a) za stejn´ e, b) za r˚ uzn´ e. V´ ysledky porovnejte.
% ˇ REˇ SEN´ I % zadan´ e hodnoty al = .05; n1 = 20; n2 = 25; alt = ’<>’; % konstrukce dat (jsou zad´ any jiˇ z vypoˇ cten´ e momenty) x1.n=n1; x1.m=mx1; x1.v=s2x1; x2.n=n2; x2.m=mx2; x2.v=s2x2; % a) sdruˇ zen´ y test is_s = t_int_2s(x1,x2,al,alt);
17
% b) nesdruˇ zen´ y test is_n = t_int_2n(x1,x2,al,alt); fprintf(’Dve stredni hodnoty\n\n’); fprintf(’interval spolehlivosti: (%g, %g)\n’, is_s); fprintf(’interval spolehlivosti: (%g, %g)\n’, is_n); % POZN´ AMKY: % Protoˇ ze v´ ybˇ erov´ e rozptyly jsou bl´ ızk´ e, jsou bl´ ızk´ e i % IS pro sdruˇ zen´ y a nesdruˇ zen´ y v´ ybˇ er % Protoˇ ze 0 neleˇ z´ ı v intervalech, lze uˇ cinit z´ avˇ er (s % pravdˇ epodobnost´ ı 0.95), ˇ ze metody nejsou shodn´ e. ****************** p045 ****************** % P045 - Interval pro pod´ ıl % ------------------------intro % % % %
ZAD´ AN´ I Pˇ R´ IKLADU Na d´ alniˇ cn´ ım ´ useku s rychlost´ ı omezenou na 80 km/h byla n´ ahodn´ ym v´ ybˇ erem u 20 automobil˚ u zmˇ eˇ rena rychlost s v´ ysledky
x=[78 89 75 78 82 98 77 73 76 78 68 79 81 79 80 ... 77 73 99 80 76]; % Urˇ cete levo, pravo i oboustrann´ y 0.95-IS pro pod´ ıl % automobil˚ u, pˇ rekraˇ cuj´ ıc´ ıch povolenou rychlost. % N´ avod: poˇ cet pˇ rekraˇ cuj´ ıc´ ıch je sum(x>80) % ˇ REˇ SEN´ I % zadan´ e hodnoty al = .05; n = 20; alt_l = ’>’; alt_p = ’<’; alt_o = ’<>’; n1=sum(x>80); p=n1/n; is_l= prop_int(p,n,al,alt_l); is_p= prop_int(p,n,al,alt_p); is_o= prop_int(p,n,al,alt_o); fprintf(’Dva podily\n\n’); fprintf(’levostrann´ y IS: fprintf(’pravostrann´ y IS:
(%6.4g, %g)\n’, is_l); (%g, %6.4g)\n’, is_p); 18
fprintf(’oboustrann´ y IS:
(%6.4g, %6.4g)\n’, is_o);
% POZN´ AMKY: % V testu lze zadat bud pod´ ıl p nebo poˇ cet n1. V´ ysledek % je stejn´ y.
15.3
Testy hypot´ ez
****************** p051 ****************** % P051 - Test hypot´ ezy pro stˇ redn´ ı hodnotu % ---------------------------------------intro % % % % %
ZAD´ AN´ I Pˇ R´ IKLADU Byl proveden v´ ybˇ er o rozsahu 16 z norm´ aln´ ıho rozdˇ elen´ ı a byl vypoˇ cten v´ ybˇ erov´ y pr˚ umˇ er 4518 a v´ ybˇ erov´ a smˇ erodatn´ a odchylka 115. Testujte hypot´ ezu H0: mi0=4500 proti alternativˇ e HA: mi>mi0 na hladinˇ e v´ yznamnosti 0.05.
% ˇ REˇ SEN´ I % zadan´ e hodnoty n = 16; mi0 = 4500; mx = 4518; sx = 115; al = .05;
% % % % %
rozsah v´ ybˇ eru H0 pr˚ umˇ er smˇ erodatn´ a odchylka pravdˇ epodobnost
t_kr=-t_inv(al,n-1); t_r=(mx-mi0)/sx*sqrt(n); W=[t_kr, inf]; pval=1-t_cdf(t_r,n-1);
% % % %
kritick´ a hodnota realizovan´ a statistika kritick´ y obor p-hodnota
fprintf(’Test stredni hodnoty\n\n’); fprintf(’statistika: %g\n’,t_r); fprintf(’krit. obor: (%g, %g)\n’,W); fprintf(’p-hodnota: %g\n’,pval); if t_r>t_kr disp(’=> H0 zam´ ıt´ ame’); else disp(’=> H0 nezam´ ıt´ ame’); end%if % POZN´ AMKY: % m-file je pˇ ripraven pro zmˇ enu zad´ an´ ı. Proto je z´ avˇ er % formulov´ an obecnˇ e (i pro zam´ ıtnut´ ı H0) % Ovˇ eˇ ren´ ı pomoc´ ı funkce: disp(’ ’); 19
disp(’Overeni’); % konstrukce v´ ybˇ eru (momenty) x.n = n; x.m = mx; x.v = sx^2; [pval, t, df]=t_test(x, mi0, ’>’)
****************** p052 ****************** % P052 - Test pro stˇ redn´ ı hodnotu % ------------------------------intro % % % % % %
´N´ ZADA I Pˇ R´ IKLADU Provedli jsme n´ ahodn´ y v´ ybˇ er o rozsahu 10 z norm´ aln´ ıho rozdˇ elen´ ı se smˇ erodatnou odchylkou 3 a vypoˇ cetli pr˚ umˇ er 25.25. Na hladinˇ e 0.05 testujte nulovou hypot´ ezu H0: mi0=24 proti alternativˇ e HA: a) mi != mi0; b) mi > mi0; c) mi < mi0.
% ˇ REˇ SEN´ I % zadan´ e hodnoty n = 10; sx = 3; mx = 25.25; mi0 = 24; al = .05; % konstrukce v´ ybˇ eru (momenty) x.n = n; x.m = mx; x.v = sx^2; % a) disp(’Oboustranny test’); pval=t_test(x, mi0) disp(’ ’); % b) disp(’Levostranny test’); pval=t_test(x, mi0, ’<’) disp(’ ’); % c) disp(’Pravostranny test’); pval=t_test(x, mi0, ’>’) % POZN´ AMKY: 20
% Lze tak´ e vyuˇ z´ ıt vnitˇ rn´ ıch tisk˚ u (zad´ an´ ı bez v´ ystupu). ****************** p053 ****************** % P053 - Test pro stˇ redn´ ı hodnotu % ------------------------------intro % % % % %
´N´ ZADA I Pˇ R´ IKLADU V´ yrobce tvrd´ ı, ˇ ze j´ ım vyr´ abˇ en´ e kuliˇ cky do loˇ zisek maj´ ı stˇ redn´ ı hodnotu pr˚ umˇ eru 10 mm a ne vˇ etˇ s´ ı. Ne hladinˇ e v´ yznamnosti 0.05 testujte uvedenou hypot´ ezu, jestliˇ ze ve v´ ybˇ eru 16 kuliˇ cek byl pr˚ umˇ er 10.3 a a) sig2=1, b) s2=1.21.
% ˇ REˇ SEN´ I % zadan´ e hodnoty alt = ’>’; mi0 = 10; mx = 10.3; n = 16; sig2 = 1; s2 = 1.21; % konstrukce v´ ybˇ eru (momenty) x.n = n; x.m = mx; x.v = s2; % a) pvala=z_test(x, mi0, sig2, alt); % b) pvalb=t_test(x, mi0, alt); fprintf(’Test stredni hodnoty\n\n’); fprintf(’pval pro sig2: %g\n’, pvala); fprintf(’pval pro s2: %g\n’, pvalb); % POZN´ AMKY: ****************** p054 ****************** % P054 - Test pro rozptyl % ----------------------intro % % % %
ZAD´ AN´ I Pˇ R´ IKLADU Nov´ a metoda mˇ eˇ ren´ ı d´ elky souˇ c´ astek byla ovˇ eˇ rov´ ana na etalonu. Disperze urˇ cen´ a z 10 mˇ eˇ ren´ ı byla 100. Je tento v´ ysledek ve shodˇ e s tvrzen´ ım, z ˇe disperze v´ ysledk˚ u nov´ e 21
% metody nen´ ı vˇ etˇ s´ ı neˇ z 50? Testujte na hladinˇ e 0.05. % ˇ REˇ SEN´ I % zadan´ e hodnoty n = 10; s2 = 100; s20 = 50; alt = ’>’; % konstrukce v´ ybˇ eru (momenty) x.n = n; x.v = s2; var_test(x, s20, alt); % POZN´ AMKY: ****************** p055 ****************** % P055 - Test pro dvˇ e stˇ redn´ ı hodnoty % ----------------------------------intro % ZAD´ AN´ I Pˇ R´ IKLADU % Zjiˇ stoval se vliv dvou kataliz´ ator˚ u na konverzi plynu. % V´ ysledky (v procentech) jsou x = [24.1 22.8 24.4 24.7 23.9]; y = [23.2 24.5 24.1 25.1 22.6 24.2 24.0]; % Ovˇ eˇ rte hypot´ ezu, ˇ ze vysledky obou katalyz´ ator˚ u jsou na % hladinˇ e v´ yznamnosti 0.05 stejn´ e. Testujte pro % a) stejn´ e rozptyly, b) r˚ uzn´ e rozptyly. % ˇ REˇ SEN´ I % zadan´ e hodnoty al = .05; alt = ’<>’; % a) stejn´ e rozptyly (sdruˇ zen´ y test) [pvala, ta, dfa] = t_test_2s (x, y, alt); % b) r˚ uzn´ e rozptyly (nesdruˇ zen´ y test) [pvalb, tb, dfb] = t_test_2n (x, y, alt); disp(’Sdruzeny t-test’) fprintf(’pval: %g\n’, pvala); if pvala
disp(’ H0 nezamitame’); end%if disp(’ ’); disp(’Nesdruzeny t-test’) fprintf(’pval: %g\n’, pvalb); if pvalb
´N´ ZADA I Pˇ R´ IKLADU Dvˇ ema laboraton´ ımi metodami se zjiˇ stoval obsah chemick´ e l´ atky v roztoku (v %). Pˇ redpokl´ ad´ ame. ˇ ze v´ ysledky maj´ ı norm´ aln´ ı rozdˇ elen´ ı. Bylo provedeno mˇ eˇ ren´ ı 5 vzork˚ u (kaˇ zd´ y vzorek byl mˇ eˇ ren obˇ ema metodami) s v´ ysledky x = [2.3 1.9 2.1 2.4 2.6]; y = [2.4 2.0 2.0 2.3 2.5];
% Ovˇ eˇ rte, zda existuje na hladinˇ e v´ yznamnosti 0.05 podstatn´ y % rozd´ ıl mezi metodami. % ˇ REˇ SEN´ I % zadan´ e hodnoty al = .05; n=length(x); if length(y)~=n, disp(’Delky x a y musi byt stejne’); return end%if [pval, t, df] = t_test_2p (x, y); fprintf(’Parovy t-test\n\n’); fprintf(’p-hodnota: %g\n’, pval); fprintf(’statistika: %g\n’, t); fprintf(’stupne vol.: %g\n’, df); % POZN´ AMKY: ****************** p057 ****************** 23
% P057 - IS a TH pro dva pod´ ıly % ----------------------------intro % % % % % % % %
´N´ ZADA I Pˇ R´ IKLADU Na dvou pracoviˇ st´ ıch sledujeme pracovn´ ı prostoje. Pracoviˇ stˇ e A bylo sledov´ ano v 800 ˇ casov´ ych okamˇ zic´ ıch a bylo zaznamen´ ano 116 prostoj˚ u. Kontrola pracoviˇ ste B prob´ ıhala v 1200 okamˇ zic´ ıch a bylo zjiˇ stˇ eno 78 prostoj˚ u. a) Stanovte 95% IS pro rozd´ ıl stˇ redn´ ıch pod´ ıl˚ u. b) Na hladinˇ e v´ yznamnosti 0.05 testujte hypot´ ezu o rovnosti stˇ redn´ ıch pod´ ıl˚ u.
% ˇ REˇ SEN´ I % zadan´ e hodnoty x1 = 116; n1 = 800; x2 = 78; n2 = 1200; al = .05; % a) is=prop_int_2(x1,n1,x2,n2,al); % b) [pval, z] = prop_test_2 (x1, n1, x2, n2); fprintf(’Dva podily\n\n’); fprintf(’interval spolehlivosti: fprintf(’p-hodnota testu:
(%g, %g)\n’, is); %g\n’, pval);
% POZN´ AMKY: % Souvislost mezi IS a TH: nula nen´ ı v IS => H0 se zam´ ıt´ a.
15.4
Chi2 testy hypot´ ez
****************** p061 ****************** % P061 - Test dobr´ e shody pro rovnomˇ ern´ e rozdˇ elen´ ı % -----------------------------------------------intro ZAD´ AN´ I Pˇ R´ IKLADU Testujte tvrzen´ ı, ˇ ze ve Sv´ edsku se dˇ eti rod´ ı rovnomˇ ernˇ e po cel´ y rok. K dispozici jsou ´ udaje z n´ ahodn´ eho v´ ybˇ eru 88 porod˚ u, sdruˇ zen´ e do ˇ ctyˇ r (r˚ uznˇ e dlohuh´ ych) obdobi. D´ elky obdob´ ı ve dnech jsou: d = [91 62 61 151]; % Poˇ cty porod˚ u za jednotliv´ a obdob´ ı: % % % % %
24
x = [27 20
8
33];
% ˇ REˇ SEN´ I n = length(x); % Pozorovan´ e ˇ cetnosti jsou zad´ any o = x; % % % % e
Teoretick´ e ˇ cetnosti urˇ c´ ıme tak, aby 1) byl stejn´ y poˇ cet porod˚ u (tj. 88), 2) porody byly rozdˇ eleny rovnomˇ ernˇ e, tj. ´ umˇ ernˇ e d´ elce obdob´ ı. Tedy = d/sum(d)*sum(x);
% Statistika: ch2 = sum((o-e).^2./e); % p-hodnota pval = 1-chisquare_cdf(ch2,n-1); % Ovˇ eˇ ren´ ı pomoc´ ı funkce ’chisquare_test’ [pv ch2]=chisquare_test(o,e); fprintf(’Test dobre shody rovnomerneho rozdeleni\n\n’); fprintf(’p-hodnota: %g\n’, pval); fprintf(’p-hodnota (funkce): %g\n’, pv); % POZN´ AMKY: ****************** p062 ****************** % P062 - Test dobr´ e shody pro norm´ aln´ ı rozdˇ elen´ ı % ---------------------------------------------intro % % % % %
ZAD´ AN´ I Pˇ R´ IKLADU V souboru ’data1’ je uloˇ zen datov´ y vzorek x. a) Urˇ cete bodov´ y odhad jeho stˇ redn´ ı hodnoty a rozptylu. b) Testujte, zda poch´ az´ ı z norm´ aln´ ıho rozdˇ elen´ ı s odhadnut´ ymi charakteristikami.
% ˇ REˇ SEN´ I load data01; n = length(x); m = 30;
% data % rozsah v´ ybˇ eru % poˇ cet interval˚ u pro ˇ cetnosti
% a) odhady mx = mean(x); vx = var(x); 25
% b) test dobr´ e shody % Pozorovan´ e ˇ cetnosti (na ’m’ intervalech) [o xx]=hist(x, m); s=xx+(xx(2)-xx(1))/2; % Teoretick´ e ˇ cetnosti f=normal_cdf(s, mx, vx); ff=[0 f]; ee=ff(2:end)-ff(1:end-1); e=sum(o)*ee; % Test [pval ch2] = chisquare_test(o, e); fig plot(xx,o,’bo’,xx,e,’rx’); title(’Dobra shoda pro norm. rozd. (blue=o, red=e)’) fprintf(’Test dobre shody pro normalni rozdeleni\n\n’); fprintf(’p-hodnota testu: %g\n’, pval); fprintf(’ (pozri tiez graf)\n’); % POZN´ AMKY: ****************** p063 ****************** % P063 - Test nez´ avislosti % -----------------------intro ´N´ ZADA I Pˇ R´ IKLADU Bylo dot´ az´ ano 400 osob na bydliˇ stˇ e a platovou skupinu. Z odpovˇ ed´ ı byla sestavena kontingenˇ cn´ ı tabulka = [28 42 30 24; 44 78 78 76]; % kde ˇ r´ adky odpov´ ıdaj´ ı bydliˇ sti a sloupce platov´ ym tˇ r´ ıd´ am. % Testujte hypot´ ezu H0: ’bydliˇ stˇ e a platy jsou nez´ avisl´ e’. % % % T
% ˇ REˇ SEN´ I % a) podle vzorc˚ u [nr nc]=size(T); ss=sum(sum(T)); t=T/ss; sc=sum(t); sr=sum(t’); tt=sr’*sc; TT=tt*ss;
% % % % % % %
rozmˇ ery tabulky celkov´ y souˇ cet normov´ an´ ı norm. souˇ cet pˇ res ˇ r´ adky (svisle) norm. souˇ cet pˇ res sloupce (vodorovnˇ e) nez´ avisl´ a normovan´ a tabulka naz´ avisl´ a tabulka absolutn´ ı 26
o=T(:); e=TT(:); n=(nr-1)*(nc-1);
% pozorovan´ e ˇ cetnosti % teoretick´ e c ˇetnosti % stupnˇ e volnosti
[pval ch2]=chisquare_test(o,e,n); % b) ovˇ eˇ ren´ ı funkc´ ı [pv ch2 df] = chisquare_test_i (T);
fprintf(’Test nez´ avislosti\n\n’); fprintf(’stat. (vzorce): %g\n’, ch2); fprintf(’pval (vzorce): %g\n\n’, pval); fprintf(’stat. (funkce): %g\n’, ch2); fprintf(’pval (funkce): %g\n’, pval); % % % %
´MKY: POZNA Kontingenˇ cn´ ı tabulka obsahuje absolutn´ ı ˇ cetnosti dvou znak˚ u: bydliˇ ste 1,2; plat 1,2,3,4. Potom napˇ r. 30 lid´ ı odpovˇ edˇ elo, ˇ ze bydl´ ı v 1 a maj´ ı plat 3.
15.5
Dalˇ s´ı neparametrick´ e testy hypot´ ez
****************** p071 ****************** % P071 - Znam´ enkov´ y test medi´ an˚ u % -----------------------------intro % % % %
ZAD´ AN´ I Pˇ R´ IKLADU Generujte dva vzorky dat d´ elky n. Oba s rozdˇ elen´ ım chi2, prvn´ ı s 5 a druh´ y s 8 stupni volnosti. Pomoc´ ı zanam´ enlov´ eho tesu ovˇ eˇ rte shodu medi´ an˚ u tˇ echto vzork˚ u.
% ˇ REˇ SEN´ I n = 100; % generov´ an´ ı vzork˚ u x=chisquare_rnd(5, 1, n); y=chisquare_rnd(8, 1, n);
% chi2 s 5 stupni volnosti % chi2 s 8 stupni volnosti
% v´ ypoˇ cet medi´ an˚ u med_x=median(x); med_y=median(y); % znam´ enkov´ y test pval=sign_test(x, y, ’<>’); fprintf(’Znamenkovy test medianu\n’); 27
fprintf(’(H0: jsou stejn´ e)\n\n’); fprintf(’p-hodnota testu:
%g\n’, pval);
% POZN´ AMKY: % vykreslen´ ı pr˚ ubˇ eh˚ u datov´ ych vzork˚ u fig plot(1:n,x,1:n,y); title(’Testovana data’) ****************** p072 ****************** % P072 - Test nez´ avislosti prvk˚ u v´ ybˇ eru (rezidu´ ı) % ----------------------------------------------intro % % % % % % %
ZAD´ AN´ I Pˇ R´ IKLADU V n´ asleduj´ ıc´ ım pˇ r´ ıkladˇ e generujeme data pomoc´ ı diferenˇ cn´ ı rovnice 2. ˇ r´ adu. Data modelujeme pomoc´ ı diferenˇ cn´ ı rovnice 1. ˇ r´ adu jej´ ıˇ z koeficienty odhadujeme metodou nejmenˇ s´ ıch ctverc˚ ˇ u. Odeˇ cten´ ım zmˇ eˇ ren´ ych dat a dat generovan´ ych z odhadnut´ e diferenˇ cn´ ı rovnice z´ ısk´ ame rezidua. Pokud je odhad dobr´ y, mus´ ı b´ yt rezidua nez´ avisl´ a. Testujte!
% ˇ REˇ SEN´ I k = 1; % ˇ rad odhadu (1 nebo 2) n = 200; % poˇ cet dat ni = 3; % zaˇ c´ atek generov´ an´ ı dat if k==1, nv = 5; % poˇ cet odhadovan´ ych parametr˚ u + 1 else nv = 6; end%if y=zeros(1,n); e=randn(1,n); u=randn(1,n); v=eye(nv)*1e-5; % statistika pro odhad for t=ni:n % simulace regr. modelu druh´ eho ˇ r´ adu y(t)=.9*y(t-1)-.5*y(t-2)+1.3*u(t)+.2*u(t-1)+.1*e(t); % identifikace regr. modelu prvn´ ıho ˇ r´ adu if k==1, d=[y(t) y(t-1) u(t) u(t-1) 1]; else d=[y(t) y(t-1) y(t-2) u(t) u(t-1) 1]; end%if v=v+d’*d; end vy=v(1,1); vfy=v(2:nv,1); vf=v(2:nv,2:nv); p=inv(vf)*vfy; % v´ ypoˇ cet odhadu ze statistiky for t=ni:n % predikce identifikovan´ ym modelem prvn´ ıho ˇ r´ adu if k==1, yp(t)=p(1)*y(t-1)+p(2)*u(t)+p(3)*u(t-1)+p(4); else yp(t)=p(1)*y(t-1)+p(2)*y(t-2)+p(3)*u(t)+p(4)*u(t-1)+p(5); end%if end 28
% chyby predikce ep=y(ni:n)-yp(ni:n); % test na bˇ elost (nez´ avislost) chyb predikce pval=wztest(ep); fprintf(’Test nezavislosti rezidui\n’); fprintf(’(H0: jsou nezavisla)\n\n’); fprintf(’p-hodnota testu: %g\n’, pval); % POZN´ AMKY: ****************** p073 ****************** % P073 - Test nez´ avislosti v´ ybˇ er˚ u (Pearson, Spearman, Kendal) % ----------------------------------------------------------intro % ZAD´ AN´ I Pˇ R´ IKLADU % Mˇ eˇ ren´ ım jsme z´ ıskali dva v´ ybˇ ery ’x’ a ’y’ x=[ 2.5 3.4 1.3 5.8 3.6 2.7 4.3 5.1 2.9 4.5]; y=[13.4 15.2 11.8 13.1 14.5 14.1 12.7 13.6 11.2 13.4]; % % % % % %
Testujte, zda poch´ azej´ ı z nez´ avisl´ ych rozdˇ elen´ ı. Pro test pouˇ zijte a) Pearson˚ uv test, b) Spearman˚ uv test, c) Kendal˚ uv test. V´ ysledky srovnejte!
% ˇ REˇ SEN´ I disp(’a) Pearsonuv test’) sP=cor_test(x,y,’<>’,’p’) disp(’a) Spearmanuv test’) sS=cor_test(x,y,’<>’,’s’) disp(’a) Kendalv test’) sK=cor_test(x,y,’<>’,’k’) % POZN´ AMKY: disp(’ *** H0: jsou nezavisle ***’); ****************** p074 ****************** % P074 - Kolmogorov-Smirnov˚ uv test typu rozdˇ elen´ ı % ----------------------------------------------intro % ZAD´ AN´ I Pˇ R´ IKLADU % Generujte data z a) norm´ aln´ ıho, b) rovnomˇ ern´ eho rozdˇ elen´ ı 29
% a testujte, zda gener´ ator skuteˇ cnˇ e pˇ redstavuje tato % rozdˇ elen´ ı % ˇ REˇ SEN´ I n = 1000;
% poˇ cet dat
% a) test rovnomˇ ern´ eho rozdˇ elen´ ı na intervalu (2,4) x=uniform_rnd(2, 4, 1, n); pv1=kolmogorov_smirnov_test(x, ’uniform’, 2, 4); % test norm´ aln´ ıho rozdˇ elen´ ı se stˇ r. h. 5 a rozpt. 9 y=normal_rnd(5, 9, 1, n); pv2=kolmogorov_smirnov_test(y, ’normal’, 5, 9);
fprintf(’Kolmogorov-Smirnov test\n\n’); fprintf(’p-hodnota pro rovnomerna data: fprintf(’p-hodnota pro normalni data:
%g\n’, pv1); %g\n’, pv2);
% POZN´ AMKY: % Dalˇ s´ ı typy rozdˇ elen´ ı. % Zad´ an´ ı je stejn´ e % gener´ ator: Rozdˇ elen´ ı_rnd(Parametry, ˇ r´ adky, sloupce) % test: kolmogorov_smirnov_test(data, Rozdˇ elen´ ı, Parametry) % % % % % % % % % % % % % %
Rozdˇ elen´ ı Parametry -----------------------------binomial n,p poisson lambda geometric p hypergeometric m,t,n uniform a,b exponential lambda lognormal a,v stdnormal normal m,v t df chisquare df f df_num,df_den
% Pˇ r´ ıklady dalˇ s´ ıch rozdˇ elen´ ı % x1=exponential_rnd(2,1,n); % pv_e=kolmogorov_smirnov_test(x1,’exponential’,2) % x2=binomial_rnd(10,.6,1,n); % pv_b=kolmogorov_smirnov_test(x2,’binomial’,10,.6) 30
% x3=poisson_rnd(4,1,n); % pv_p=kolmogorov_smirnov_test(x3,’poisson’,4) % x4=lognormal_rnd(10,9); % pv_l=kolmogorov_smirnov_test(x4,’lognormal’,10,9)
15.6
Regresn´ı anal´ yza
****************** p081 ****************** % P081 - Regresn´ ı pˇ r´ ımka % ---------------------intro % % % % % a b s x
´N´ ZADA I Pˇ R´ IKLADU Generujte datov´ e dvojice ’[x,y]’, pro ’x=1,2,...,10’ tak, aby splnovali rovnici ’y = ax + b + e’, kde ’e’ je norm´ aln´ ı n´ ahodn´ a veliˇ cina s nulovou stˇ redn´ ı hodnotou a smˇ erodatnou odchylkou ’s’. Volte: = 2; = 1; = 1; = 1:10;
% % % % % % % % % %
a) b) c) d) e) f) g)
Urˇ cete koeficienty regresn´ ı pˇ r´ ımky. Urˇ cete korelaˇ cn´ ı koeficient. Vypoˇ ctˇ ete predikce ve vˇ sech bodech ’x’. Body ’[x,y]’ i regresn´ ı pˇ r´ ımku zobrazte. Urˇ cete pˇ redpov´ ıdanou hodnotu ’ya’ pro ’xa=2.5’. Odhadnˇ ete, pro kter´ e ’xb’ bude platit ’yb=16’. Doporuˇ cen´ e experimenty: 1. Mˇ ente ’a’ na -1, -.1, .1, .5, 1, 5, 15 2. Mˇ ente ’s’ na 0, .1, 10, 100 3. Provedte pˇ redpovˇ ed ’yc’ pro ’xc=25, 125’
% ˇ REˇ SEN´ I n=length(x); % poˇ cet dat e=s*randn(1,n); % ˇ sum y=a*x+b+e; % hodnoty ’y’ mx=mean(x); % pr˚ umˇ er ’x’ my=mean(y); % pr˚ umˇ er ’y’ s2x=var(x)*(n-1); % souˇ cet ˇ ctverc˚ u pro ’x’ sxy=cov(x,y)*(n-1); % souˇ cet ˇ ctverc˚ u pro ’x,y’ s2y=var(y)*(n-1); % souˇ cet ˇ ctverc˚ u pro ’y’ % a) disp(’Koeficienty regresni primky’) b1=sxy/s2x % smˇ ernice b0=my-b1*mx % absolutn´ ı ˇ clen 31
% b) disp(’Korelacni koeficient’) r=sxy/sqrt(s2x*s2y) % korelaˇ cn´ ı koeficient % c) yp=b1*x+b0; % predikce % d) if 1, axis([0,11,0,30]); % osy else axis([min(x)-1,max(x)+1,min(y)-1,max(y)+1]); end fig plot(x,y,’o’,x,yp,’+-g’); % graf title(’Regesni analyza pro testovana data’) % e) disp(’Predikce y pro dane x’) xa=2.5; ya=b0+b1*xa % ’y’ pro dan´ e ’x’ % f) disp(’Predukce x pro dane y’) yb=16; xb=(yb-b0)/b1 % ’x’ pro dan´ e ’y’ % POZN´ AMKY: ****************** p082 ****************** % P082 - Odlehl´ a a vlivn´ a pozorov´ an´ ı % ---------------------------------intro % % x y % %
ZAD´ AN´ I Pˇ R´ IKLADU Provedte regresn´ ı anal´ yzu pro data = [1 2 3]; = [1 2.3 2.8]; Pˇ ridejte dalˇ s´ ı mˇ eˇ ren´ ı ’[x4,y4]’ a v´ ysledek porovnejte s p˚ uvodn´ ım v´ ysledkem.
% REˇ SEN´ I ii=3;
% kl´ ıˇ c, kter´ a data se pˇ ridaj´ ı % ii=1 (nic), ii=2 (odlehl´ e), ii=3 (vlivn´ e)
% ii=2 - odlehl´ a pozorov´ an´ ı if ii==2, x(4)=4; y(4)=15; end % pˇ ridan´ e mˇ eˇ ren´ ı % ii=3 - vlivn´ a pozorov´ an´ ı if ii==3, x(4)=9; y(4)=10; end % pˇ ridan´ y outlier (spr´ avnˇ e y=8.33) disp(’Koeficienty regrese’) [b1 b0 r]=reg_desc(x,y) yp=b0+b1*x; fig plot(x,y,’o’,x,yp); axis([min(x)-.1,max(x)+.1,min(y)-.1,max(y)+.1])
32
% POZN´ AMKY: ****************** p083 ****************** % P083 - M´ ıra vhodnosti dat pro regresi % ------------------------------------intro % % % %
ZAD´ AN´ I Pˇ R´ IKLADU Generujte data na elipse a provedte jejich regresn´ ı anal´ yzu. Pˇ ri experimentech mˇ ente velikosti poloos a otoˇ cen´ ı elipsy. Sledujte hodnotu korelaˇ cn´ ıho koeficientu ’r’.
% ˇ REˇ SEN´ I s=0; a=10; b=1;
% amplituda ˇ sumu pˇ ridan´ eho k elipse % poloosy elipsy % stejn´ e - nevhodn´ a, r˚ uzn´ e - vhodn´ a
d=[]; for f=0:.1:2*pi % pol´ arn´ ı souˇ radnice xf=a*cos(f)+s*randn; yf=b*sin(f)+s*randn; d=[d; [xf yf]]; end d=d*[cos(.1) sin(.1); -sin(.1) cos(.1)]; % otoˇ cen´ ı souˇ radnic x=d(:,1); y=d(:,2); [b1 b0 r]=reg_desc(x,y); % regrese yp=b0+b1*x; % predikce fig plot(x,y,’o’,x,yp) title(’Data vhodna / nevhodna pro regresi’) axis([min(x)-1,max(x)+1,min(y)-1,max(y)+1]); fprintf(’ Korelacni koeficient je: %5.3g\n’,r); % POZN´ AMKY:
15.7
Korelaˇ cn´ı anal´ yza
****************** p091 ****************** % P091 - Korelaˇ cn´ ı anal´ yza (linearni a exponencialni regrese) % ----------------------------------------------------------intro % % % %
´N´ ZADA I Pˇ R´ IKLADU Sledujeme v´ yvoj poˇ ctu automobil˚ u na jednoho obyvatele v nasi republice ve vybran´ ych roc´ ıch. Zjiˇ stˇ en´ e ´ udaje jsou (x - rok, y - poˇ cet/hlavu) x = [ 55 58 63 71 75 82 89 93 95 98 99]; y = [.42 .23 .35 .21 .35 .42 .66 .83 .95 1.2 .98]; 33
% % % % %
a) Provedte line´ arn´ ı regresn´ ı anan´ yzu tˇ echto dat a testujte vhodnost line´ arn´ ı regrese. b) Provedte exponenci´ aln´ ı regresn´ ı anal´ yzu a jej´ ı v´ ysledky porovnejte s s v´ ysledky line´ arn´ ı regrese. Pro porovn´ an´ ı pouˇ zijte koeficient determinace.
% ˇ REˇ SEN´ I % a) xp = 99; al = .05; [b1,b0,r]=reg_desc(x,y); [ie,ip,pa,pr]=reg_infe(x,y,xp,al); pv=f_test_reg(x,y); disp(’ ’) disp(’ a) Linearni regrese’); disp(’ ****************’); fprintf(’Koeficienty regresni primky: %g a %g\n’,b1,b0); fprintf(’Korelacni koeficient: %g\n’,r); fprintf(’Interval str. h. pro xp=%d je: (%g, %g)\n’,xp,ie); fprintf(’Interval predikce pro xp=%d je: (%g, %g)\n’,xp,ip); fprintf(’p-hodnota pro t-test b1: %g\n’,pa); fprintf(’p-hodnota pro t-test r: %g\n’,pr); fprintf(’p-hodnota pro F-test: %g\n’,pv); % b) p=exp_reg(x,y); yp=exp_pred(x,p); p_exp=f_test_pred(y,yp); disp(’ ’) disp(’ b) Exponencialni regrese’); disp(’ *********************’); fprintf(’Koeficienty exp. regrese: fprintf(’p-hodnota pro F-test:
%g a %g\n’,p); %g\n’,p_exp);
% graf linearni a exponencialni regrese yy=lin_pred(x,[b1 b0]); fig plot(x,y,’kx’,x,yy,’r’,x,yp,’b’); title(’Linearni a exponencialni regrese’) % POZN´ AMKY: % Pro toto zad´ an´ ı je exponenci´ aln´ ı regrese lepˇ s´ ı % (m´ a menˇ sı ´ p-hodnotu - ide´ aln´ ı je ph=0) ****************** p092 ****************** 34
% P092 - Z´ akladn´ ı inference v line´ arn´ ı regresi % -------------------------------------------intro % % % % % a b s x
´N´ ZADA I Pˇ R´ IKLADU Generujte datov´ e dvojice ’[x,y]’, pro ’x=1,2,...,10’ tak, aby splnovali rovnici ’y = ax + b + e’, kde ’e’ je norm´ aln´ ı n´ ahodn´ a veliˇ cina s nulovou stˇ redn´ ı hodnotou a smˇ erodatnou odchylkou ’s’. Volte: = 2; = 1; = 3; = 1:10;
% % % % % %
a) b) c) d) e) f)
Provedte line´ arn´ ı regresn´ ı anal´ yzu. Body ’[x,y]’ i regresn´ ı pˇ r´ ımku zobrazte. Stanovte 0.05-interval pro smˇ ernici reg. pˇ r´ ımky. Stanovte 0.05-interval pro regresn´ ı pˇ r´ ımku v ’xp=6’. Testem korelaˇ cn´ ıho koeficientu ovˇ eˇ rte vhodnost regrese. Provedte F-test regrese a porovnejte t-testem.
% ˇ REˇ SEN´ I n=length(x); e=s*randn(1,n); y=a*x+b+e;
% poˇ cet dat % ˇ sum % hodnoty ’y’
[b1 b0 r]=reg_desc(x,y); yp=lin_pred(x,[b1 b0]); if 1, axis([0,11,0,30]); % osy else axis([min(x)-1,max(x)+1,min(y)-1,max(y)+1]); end fig plot(x,y,’o’,x,yp,’+-g’); % graf hold on xp=6; alpha=.05; alt=’<>’; disp(’ is_e interval spolehlivosti pro stredni hodnotu’) disp(’ is_r interval spolehlivosti pro predikci’) disp(’ pval_a p-hodnota pro test "smernice reg. pr."=0’) disp(’ pval_p p-hodnota pro test "korelacni koef."=0’) disp(’ ’) [is_a,is_e,pval_a,pval_r]=reg_infe(x,y,xp,alpha,alt) plot(xp,is_e(1),’xb’,xp,is_e(2),’xb’); title(’Linearni regrese’) hold off % POZN´ AMKY:
35
****************** p093 ****************** % P093 - P´ as spolehlivosti % -----------------------intro % % % % % a b s x %
ZAD´ AN´ I Pˇ R´ IKLADU Generujte datov´ e dvojice ’[x,y]’, pro ’x=1,2,...,10’ tak, aby splnovali rovnici ’y = ax + b + e’, kde ’e’ je norm´ aln´ ı n´ ahodn´ a veliˇ cina s nulovou stˇ redn´ ı hodnotou a smˇ erodatnou odchylkou ’s’. Volte: = 2; = 1; = 10; = 1:10; Urˇ cete p´ as spolehlivosti pro regresn´ e pˇ r´ ımku.
% ˇ RESEN´ I n=length(x); % poˇ cet dat e=s*randn(1,n); % ˇ sum y=a*x+b+e; % hodnoty ’y’ mx=mean(x); % pr˚ umˇ er ’x’ my=mean(y); % pr˚ umˇ er ’y’ s2x=var(x)*(n-1); % souˇ cet ˇ ctverc˚ u pro ’x’ sxy=cov(x,y)*(n-1); % souˇ cet ˇ ctverc˚ u pro ’x,y’ s2y=var(y)*(n-1); % souˇ cet ˇ ctverc˚ u pro ’y’ disp(’Koeficienty regrese’) b1=sxy/s2x % smˇ ernice b0=my-b1*mx % absolutn´ ı ˇ clen r=sxy/sqrt(s2x*s2y) % korelaˇ cn´ ı koeficient se=sqrt(s2y-b1*sxy)/(n-2); % standard error t_kr=t_inv(.95,n-2); % kritick´ a t-hodnota yp=[]; y1=[]; y2=[]; for xi=x ypi=b0+b1*xi; % predikce yp=[yp ypi]; del=t_kr*se*sqrt(1/n+(xi-mx)^2/s2x); y1=[y1 ypi-del]; % intervaly y2=[y2 ypi+del]; % spolehlivosti end if 1, axis([0,11,0,30]); % osy else axis([min(x)-1,max(x)+1,min(y)-1,max(y)+1]); end fig plot(x,y,’o’,x,yp,’+-g’); % graf hold on plot(x,y1,x,y2); title(’Pas spolehlivosti’) hold off
36
% POZN´ AMKY: ****************** p101 ***************** % P101 - Regresn´ ı anal´ yza % ----------------------intro % % % % % % % % % %
V letech 1999, 2000, 2001, 2002 a 2003 byl zaznamen´ av´ an poˇ cet obyvatel, vyuˇ z´ ıvaj´ ıc´ ı praˇ zk´ e metro. Zmˇ eˇ ren´ e ´ udaje byly pˇ repoˇ c´ ıt´ any na realtivn´ ı hodnoty (poˇ cet obyvatel pˇ repraven´ ych za 1 vteˇ rinu) a jsou y = [2.1 2.5 2.4 2.8 3.6]; Vypoˇ ctˇ ete a) line´ arn´ ı regresi, b) polynomickou regresi stupnˇ e 2 a 3 a c) exponenci´ aln´ ı regresi, d) pro vˇ sechny regrese urˇ cete predikci osob pˇ repraven´ ych za jeden den pro rok 2010.
% ˇ REˇ SEN´ I x = [-1 0 1 2 3]; % a) disp(’Koeficienty linearni regrese’) p_lin=lin_reg(x,y) % b) disp(’Koeficienty polynomialni regrese radu 2’) p_pol_2=pol_reg(x,y,2) disp(’Koeficienty polynomialni regrese radu 3’) p_pol_3=pol_reg(x,y,3) % c) disp(’Koeficienty exponencialni regrese’) p_exp=exp_reg(x,y) % d) xp = 10; k = 60*60*24; disp(’Predikce s linearni regresi’) y_lin= k*lin_pred(xp,p_lin) disp(’Predikce s polynomialni regresi radu 2’) y_pol_2=k*pol_pred(xp,p_pol_2) disp(’Predikce s polynomialni regresi radu 3’) y_pol_3=k*pol_pred(xp,p_pol_3) disp(’Predikce s exponencialni’) y_exp= k*exp_pred(xp,p_exp) % POZN´ AMKY: ****************** p102 ***************** % P102 - V´ ıcen´ asobn´ a regrerse pro re´ aln´ a data 37
% ------------------------------------------intro % % % % % % % %
ZAD´ AN´ I Pˇ R´ IKLADU V strahovsk´ em tunelu nyly mˇ eˇ reny hustoty a intenzity dopravn´ ıho proudu na 10 vybran´ ych m´ ıstech. Mˇ eˇ ren´ ı se prov´ adˇ elo kaˇ zd´ ych 5 min. Zmˇ eˇ ren´ a data jsou v souboru ’data52’. a) Provedte regresi prvn´ ı mˇ eˇ ren´ e veliˇ ciny v z´ avislosti na tˇ rech dalˇ s´ ıch veliˇ cin´ ach. b) Testujte vhodnost t´ eto regrese.
% ˇ REˇ SEN´ I load data52 ii=2000:-1:1600; id=0; y=x(1,ii)’; x=x(2:4,ii-id)’; p=lin_reg_n(x,y); yp=lin_pred_n(x,p); np=length(p);
% % % % % % % %
f_test_pred(y,yp,np); wztest(y-yp);
% F-test regrese % test na nez´ avislost rezidu´ ı
2000=souˇ casnost a my bereme 400 vzork˚ u zpoˇ zdˇ en´ ı mezi ’y’ a ’x’ nez´ avisle promˇ enn´ a z´ avisle promˇ enn´ e odhad parametr˚ u predikce poˇ cet odhadnut´ ych parametr˚ u
fig plot(1:length(y),y,’ro’,1:length(yp),yp,’cx’); title(’Vicenasobna regrese pro realna data’) ylabel(’data (r), predikce (b)’) xlabel(’cas vzorkovani’) % POZN´ AMKY: ****************** p103 ***************** % P103 - V´ ıcen´ asobn´ a regrerse pro simulovan´ a data % ----------------------------------------------intro % % % % % % % %
´N´ ZADA I Pˇ R´ IKLADU Generujte 4 nez´ avisl´ e, norm´ alnˇ e rozloˇ zen´ e n´ ahodn´ e veliˇ ciny ’x’ a na nich line´ arnˇ e z´ avislou n´ ahodnou veliˇ cinu ’y’ s koeficienty ps=[.2 3 1.4 .1 8]’; kde ps(5)=8 je absolutn´ ı ˇ clen. Tj. generuj´ ıc´ ı rovnice je y=p(1)*x1+..+p(4)*x4+p(5)+e. a) Provedte line´ arn´ ı regresi ’y’ v z´ avislosti ’x’. b) Testujte vhodnost t´ eto regrese.
38
% ˇ REˇ SEN´ I n=1000; a=10; np=length(ps); x=randn(n,np-1); d=[x ones(n,1)]; y=d*ps+a*randn(n,1); p=lin_reg_n(x,y); yp=lin_pred_n(x,p);
% % % % % % % %
f_test_pred(y,yp,np); wztest(y-yp);
% F-test regrese % test nez´ avislosti rezidu´ ı
velikost v´ ybˇ eru smˇ er. odch. poruchy poˇ cet parametr˚ u simulace nez´ avisle promˇ enn´ e data z´ avisle promˇ enn´ a odhad parametr˚ u predikce
fig plot(1:length(y),y,’ro’,1:length(yp),yp,’cx’); title(’Vicenasobna regrese pro simulovana data’) ylabel(’data (r), predikce (b)’) xlabel(’cas vzorkovani’) ****************** p104 ***************** % P104 - V´ ıcen´ asobn´ a regrerse pro simulovan´ a data % ----------------------------------------------intro ´N´ % ZADA I Pˇ R´ IKLADU % Generujte 2 nez´ avisl´ e veliˇ ciny ’x’ a na nich line´ arnˇ e % z´ avislou n´ ahodnou veliˇ cinu ’y’ s koeficienty ps=[.2 3 5]’; % kde ps(3) je absolutn´ ı ˇ clen. Tj. generuj´ ıc´ ı rovnice % je y=p(1)*x1+p(2)*x2+p(3)+e. % a) Provedte line´ arn´ ı regresi ’y’ v z´ avislosti ’x’. % b) Testujte vhodnost t´ eto regrese. % !!! nebezpec´ ı kolinearity x !!! % ˇ REˇ SEN´ I n=100; a=.1; b=.00001; np=length(ps); x1=(1:n)’+b*rand(n,1); x2=(n:-1:1)’*.1; x=[x1 x2]; d=[x ones(n,1)]; y=d*ps+a*randn(n,1); p=lin_reg_n(x,y); yp=lin_pred_n(x,p);
% % % % % % % % % % % %
poˇ cet dat smˇ er. odch. ˇ sumu nez´ avislost x1 a x2 (dan´ a ˇ sumem s amplitudou b) poˇ cet parametr˚ u prvn´ ı nez´ avisle promˇ enn´ a druh´ a nez´ avisle promˇ enn´ a nez´ avisle promˇ enn´ e data (do regrese) z´ avisle promˇ enn´ a odhad parametr˚ u predikce 39
f_test_pred(y,yp,np); wztest(y-yp);
% F-test regrese % test nez´ avislosti rezidu´ ı
fig plot(1:length(y),y,’ro’,1:length(yp),yp,’cx’); title(’Vicenasobna regrese s linearni vabou v regresnim vektoru’) ylabel(’data (r), predikce (b)’) xlabel(’cas vzorkovani’)
% POZN´ AMKY:
15.8
ANOVA
****************** p111 ***************** % P111 - Anal´ yza rozptylu ANOVA (jednoduch´ e tˇ r´ ıdˇ en´ ı) % -------------------------------------------------intro % % % %
´N´ ZADA I Pˇ R´ IKLADU Pro automobily urˇ cit´ eho typu byla zjiˇ stov´ ana minim´ aln´ ı spotˇ reba pohonn´ ych hmot. Na tˇ rech vybran´ ych automobilech bylo provedeno pˇ et mˇ eˇ ren´ ı s v´ ysledky
x1=[5.32 5.24 5.47 4.98 5.16]; % prvn´ ı soubor x2=[5.88 5.31 4.86 5.45 5.12]; % druh´ y soubor x3=[5.32 4.21 5.44 5.33 5.24]; % tˇ ret´ ı soubor % Testujte tvrzen´ ı ’vˇ sechny automobily jsou na tom se spotˇ rebou % stejnˇ e’. % ˇ REˇ SEN´ I x=[x1’ x2’ x3’]; [n,a]=size(x); m=mean(x); s=var(x); sx=var(m); sp=sum(s)/a; F=n*sx/sp; pv1=1-f_cdf(F,a-1,a*(n-1));
% % % % % % % %
datov´ a matice poˇ cet mˇ eˇ ren´ ı / automobil˚ u stˇ ren´ ı hodnoty soubor˚ u x1,x2,x3 rozptyly soubor˚ u rozptyl stˇ ren´ ıch hodnot normov´ an´ ı statistika p-hodnota
fprintf(’Testovani spotreby automobilu\n’); fprintf(’H0: spotreby jsou stejne\n\n’) fprintf(’p-hodnota testu: %g\n’,pv1); % POZN´ AMKY: % Hotov´ e funkce jsou: 40
%% ekvivalentn´ ı zad´ an´ ı % pv2=anova(x) %% zad´ an´ ı vhodn´ e pro r˚ uznˇ e dlouh´ e soubory % y=[x1 x2 x3]; % g=[ones(n,1); 2*ones(n,1); 3*ones(n,1)]; % [PVAL, F, DF_B, DF_W] = anova (y,g) ****************** p112 ***************** % P112 - Anal´ yza rozptylu ANOVA (dvojn´ e tˇ r´ ıdˇ en´ ı) % ---------------------------------------------intro % ZAD´ AN´ I Pˇ R´ IKLADU % Na tˇ rech stroj´ ıch pracuje pˇ et oper´ ator˚ u. Byl zjiˇ stov´ an hodninov´ y % v´ ykon stroj˚ u pˇ ri obsluze jednotliv´ ymi oper´ atory. Namˇ eˇ ren´ e % hodnoty jsou uvedeny v tabulce (sloupce - stroje, ˇ r´ adky % oper´ atoˇ ri) V=[ 53 61 51; 47 55 51; 46 52 49; 50 58 54; 49 54 50 ]; % Zjistˇ ete, zda na hladinˇ e 0,05 existuj´ ı rozd´ ıly mezi % stroji a mezi oper´ atory. % ˇ REˇ SEN´ I anova_2(V); % POZN´ AMKY: % Pro zobrazen´ ı v´ ysledk˚ u jsme vyuˇ zili vnitˇ rn´ ı tisk. % Stroje (sloupce): % % % Oper´ atoˇ ri (ˇ r´ adky): % %
je-li prvn´ ı p-hodnota (pro sloupce) menˇ s´ ı neˇ z 0.05, jsou podtatn´ e rozd´ ıly mezy stroji je-li druh´ a p-hodnota (pro ˇ r´ adky) menˇ s´ ı neˇ z 0.05, jsou podtatn´ e rozd´ ıly mezy oper´ atory
****************** p113 ***************** % P113 - Predikce simulovan´ ych dat z line´ arn´ ı regrese % --------------------------------------------------intro % ZAD´ AN´ I Pˇ R´ IKLADU 41
% % % %
Ovˇ eˇ rte kvalitu predikce pomoc´ ı line´ arn´ ı regrse. Data simulujte pomoc´ ı statick´ eho regresniho modelu s deseti nez´ avisle promˇ enn´ ymi, generovan´ ymi jako standardn´ ı b´ ıl´ y sum. Koeficienty regrese jsou ˇ p = [.1 2 .05 .2 -3 1 -.1 .01 -.05 .2]; % a smˇ erodatn´ a odchylka ˇ sumu s = 0.2; % Vytisknˇ ete v´ ysledky F-testu a predikci zobrazte v grafu. % ˇ REˇ SEN´ I nt=1000; % poˇ cet dat Y=[]; F=[]; for t=1:nt % simulace e=s*randn; % ˇ sum x=randn(1,10); % nez´ av. promˇ enn´ a y=p*x’+e; % z´ av. promˇ enn´ a % datov´ a matice Y=[Y; y]; % v´ ystup F=[F; x]; % datov´ a matice end Fi=F(:,[2 5 6]); % data pro odhad disp(’Odhad parametru’) ps=inv(Fi’*Fi)*Fi’*Y % odhad parametr˚ u Yp=Fi*ps; % predikce ep=Y-Yp; % chyba predikce disp(’Hodnata kriteria nejmensich ctvercu’) J=ep’*ep/nt % kriterium f_test_pred(Y,Yp,length(ps)); % f-test fig plot(1:nt,Y,’r.’,1:nt,Yp,’c.’) title(’Predikce (b) realnych dat (r)’) ****************** p114 ***************** % P114 - Predikce re´ aln´ ych dat z line´ arn´ ı regrese % ----------------------------------------------intro % % % % % %
´N´ ZADA I Pˇ R´ IKLADU Ovˇ eˇ rte kvalitu predikce pomoc´ ı line´ arn´ ı regrse. Pro regresi pouˇ zijte re´ aln´ a data ze souboru "data52" - obsazenosti a intenzity ze strahovsk´ eho tunelu. Jako z´ avisle promˇ ennos volte data z prv´ eho r´ ˇ adku a za nez´ avisle promˇ enn´ e vezmˇ ete data z 2 aˇ z 5 ˇ r´ adku. Zobrazte v´ ysledky F-testu a predikci ukaˇ zte v grafu.
% ˇ REˇ SEN´ I nt=2000; load data52
% poˇ cet dat % data z disku 42
Y=x(1,1:nt)’; Fi=x(2:5,1:nt)’ ; % datov´ a matice disp(’Odhad parametru’) ps=inv(Fi’*Fi)*Fi’*Y % parametry Yp=Fi*ps; % predikce ep=Y-Yp; % chyba predikce disp(’Hodnata kriteria nejmensich ctvercu’) J=ep’*ep/nt % kriterium f_test_pred(Y,Yp,length(ps)); % f-test regrese fig plot(1:nt,Y,’r.’,1:nt,Yp,’c.’) % grag title(’Predikce (b) realnych dat (r)’) % POZN´ AMKY:
15.9
Anal´ yza hlavn´ıch komponent
****************** p121 ***************** % P121 - Redukce datovych velicin pomoci rozkladu % kovariancni matice (vlastni cisla a vektory) %-----------------------------------------------intro % % % % %
´N´ ZADA I Pˇ R´ IKLADU V souboru ’data52’ jsou v matici ’x’ uloˇ zena data (hustoty a intenzity dopravn´ ıho proudu - po ˇ r´ adc´ ıch) ze strahovsk´ eho tunelu. Provedte redukci pomoc´ ı rozkladu kovarianˇ cn´ ı maice pro prvn´ ıch pˇ et veliˇ cin souboru.
% ˇ REˇ SEN´ I load data52 al=.95; n=5; nd=1000; D=x(1:n,1:nd); pca_eig(D,al);
% % % % %
podil zachovaneho rozptylu pocet velicin pocet dat datova matice redukce
% POZN´ AMKY: ****************** p122 ***************** % P122 - Redukce datovych velicin pomoci rozkladu % datove matice (svd rozklad) % ----------------------------------------------intro % % % %
ZAD´ AN´ I Pˇ R´ IKLADU V souboru ’data52’ jsou v matici ’x’ uloˇ zena data (hustoty a intenzity dopravn´ ıho proudu - po ˇ r´ adc´ ıch) ze strahovsk´ eho tunelu. Provedte redukci pomoc´ ı SVD 43
% rozkladu datov´ e matice pro prvn´ ıch pˇ et veliˇ cin souboru. % ˇ REˇ SEN´ I load data52 al=.9; n=100; y=x(1,1:n); d=x(2:6,1:n); D=[y’ d’]; pca_svd(D,al);
% % % % % %
podil zachovaneho rozptylu pocet dat nezavisle promenna zavisle promenne rozsirena datova matice redukce
% POZN´ AMKY: ****************** p123 ***************** % P123 - Redukce ˇ r´ adu regresn´ ıho modelu % a ovˇ eˇ ren´ ı pomoc´ ı predikce % ------------------------------------intro % % % % %
ZAD´ AN´ I Pˇ R´ IKLADU a) Generujte data regresn´ ım modelem y(t)=b(1)*u(t)+a(1)*y(t-1)+b(2)*u(t-1)+a(2)*y(t-2)+ +b(3)*u(t-2)+a(3)*y(t-3)+b(4)*u(t-3)+k+e(t) kde a=[.95 -.2 .1]; b=[.2 .1 .5 -.1]; k=.3;
% b) Provedte redukci dat pomoc´ ı rozkladu datov´ e matice % a utˇ cete transformaˇ cn´ ı matici pro tuto redukci. % % % %
c) Dalˇ s´ ıch 1000 vzork˚ u generovan´ ych ze stejn´ eho regresn´ ıho modelu redukujte pomoc´ ı nalezen´ e transformaˇ cn´ ı matice a odhadujte redukovan´ e regresn´ ı koeficienty.
% d) Kvalitu odhadu z redukovan´ ych dat ovˇ eˇ rte porovn´ an´ ım % predikce redukovn´ ym modelem a simulovan´ ych dat z % regresn´ ıho modelu.
% ˇ REˇ SEN´ I na=length(a); nb=length(b)-1; ns=max(na,nb); n2=(1:na+nb+2)+1; ps=[a b k]; al=.9; ni=100; 44
y=zeros(1,ni); u=y; e=randn(1,ni); % datov´ a matice z prvn´ ıch 100 vzork˚ u a jejich redukce D=[]; for t=ns+1:ni u(t)=rand; f=[y(t-1:-1:t-na) u(t:-1:t-nb) 1]; y(t)=ps*f’+.1*e(t); D=[D; [y(t) f]]; end%for [i,dd,sn,DD,p]=pca_svd(D,al); n=1000; y=zeros(1,n); u=y; e=randn(1,n); % odhad parametr˚ u z redukovan´ ych dat (q) Dr=[]; for t=ns+1:n u(t)=rand; f=[y(t-1:-1:t-na) u(t:-1:t-nb) 1]; y(t)=ps*f’+.1*e(t); Dr=[Dr; [y(t) f*dd]]; end%for w=DD’*DD+Dr’*Dr; n1=2:length(w); q=inv(w(n1,n1))*w(n1,1); % predikce z redukovan´ eho modelu Dry=Dr(:,1); Drf=Dr(:,n1); yp=zeros(1,n-ns); for t=1:n-ns yp(t)=q’*Drf(t,:)’; end%for % porovn´ an´ ı predikce s daty disp(’Normovany soucet ctvercu chyb predikce’) ep=Dr(:,1)’-yp; pe=ep*ep’/n fig plot(1:length(Dry),Dry,1:length(yp),yp) title(’Data a jejich predikce’) % POZN´ AMKY:
45
cv01.doc
Čebyševova nerovnost Nechť X je náhodná veličina se střední hodnotou E[ X ] = µ a rozptylem D[ X ] = σ 2 . Pak platí
P ( X − E[ X ] ≥ ε) ≤
D[ X ] . ε2
Pro výběrový průměr X :
P ( X − µ ≥ ε) ≤
σ2 . nε 2
1. Příklad Jak velký musí být rozsah náhodného výběru, z náhodné veličiny se střední hodnotou µ a směrodatnou odchylkou σ , aby se výběrový průměr x nelišil od střední hodnoty µ o více než σ / 2 s pravděpodobností 0,9. Řešení:
σ σ2 P X − µ ≤ ≥ 1 − = 0,9 2 2 σ n 2 1−
1 1 n 2
2
= 0,9 ⇒ n = 40
Poznámka: Z Čebyševovy nerovnosti nás „zajímá“ jen pravá strana.
2. Příklad Počet aut přijíždějících na křižovatku v určitém časovém okamžiku má Poissonovo rozdělení se střední hodnotou 120. Určete dolní hranici pravděpodobnosti, že skutečný počet aut, která do křižovatky v tomto intervalu přijedou, bude větší než 100 a menší než 140. Řešení: X počet aut,
µ = σ 2 = 120 .
P ( X − µ ≤ 20) ≥ 1 −
120 = 0,7 . 20 2
3. Příklad Průměrná spotřeba masa v jedné domácnosti za rok je 160 kg. Jaká je pravděpodobnost, že spotřeba přesáhne 320 kg, jestliže směrodatná odchylka spotřeby je 40 kg? Řešení: X spotřeba masa.
46
P ( X − 160 ≥ 160) ≤ Pro symetrické rozdělení je
40 2 1 = 2 16 160
cv02a1.doc
Pravděpodobnosti kvantily a kritické hodnoty Definice kvantilu ζ α : P ( X < ζ α ) = α , kritické hodnoty Zobrazení pomocí hustoty pravděpodobnosti
α -5
-4
-3
-2
z α : P( X > z α ) = α
α -1
0
1
ζα z1−α
2
3
4
zα ζ 1−α
Z obrázku je patrné, že platí (jen pro symetrické rozdělení)
ζ α = −ζ 1−α , ζ α = z1−α ,
z α = − z1−α z α = ζ 1−α
Pomocí těchto vzorečků lze v tabulkách kvantilů (kde jsou jen kladné hodnoty) hledat i záporné kvantily, nebo kritické hodnoty.
Pravděpodobnost normované náhodné veličiny Z 1.
P ( Z < 1,8) = ?
[0,96]
2.
P ( Z > 1,8) = ?
[= 1 − P ( Z < 1,8) = 1 − 0,96 = 0,04]
3.
P ( Z < −1,8) = ?
[= 1 − P ( Z < 1,8) = 0,04]
4.
P ( Z > −1,8) = ?
[= P( Z < 1,8) = 0,96]
Pravděpodobnost nenormované náhodné veličiny X , E[ X ] = µ, D[ X ] = σ 2 Normování
Z= 47
X −µ σ
Je dána náhodná veličina X se střední hodnotou µ = 5 a rozptylem σ 2 = 9 . Určete následující pravděpodobnosti
5
cv02a2.doc
Pravděpodobnosti a kritické hodnoty 1. Příklad Nechť Z je normálně rozdělená náhodná veličina se střední hodnotou 0 a rozptylem 1. Určete 1. P( Z > 1,96 ) = ?
[0,025]
2. P Z > 1,96 = ?
[2.0,025 = 0,05]
3. P ( Z ≤ 1,96 ) = ?
[ 1 − P ( Z > 1,96 ) = 0,975 ]
4. P( Z > ?) = 0,05
[1,645]
5. P Z > ? = 0,05
[1,96]
(
(
)
)
6. P( Z ≤ ?) = 0,05
[-1,645]
Řešení: Excel Výpočet pravděpodobnosti a 1.96 P(Z>a)=? 0.024998 =1-NORMSDIST(B1) 0.024998 =1-NORMDIST(B1;0;1;1) Určení kritické hodnoty p 0.05 P(Z>?)=p 1.644853 1.644853
=-NORMSINV(B5) =-NORMINV(0.05;0;1)
2. Příklad Nechť X je normálně rozdělená náhodná veličina se střední hodnotou µ = 5 a směrodatnou odchylkou σ = 3 . Určete
P( X > 10.88) = ?
[0,025]
Řešení:
Z=
X −µ = 1,96 ⇒ p = 0,025 σ
Excel m s a
5 3 10.88
3. Příklad
Z p
1.96 0.024998 0.024998
=(B3-B1)/B2 =1-NORMSDIST(E2) =1-NORMDIST(B3;B1;B2;1)
(
)
Uvažujme výběr X o velikosti 50 z náhodné veličiny X ≈ N µ = 5, σ 2 = 84 . X je výběrový průměr a 48 S je výběrový součet. Určete 1. P ( X > 10 ) = ?
(
)
[0,0001] [0,22]
cv02b.doc
Pravděpodobnost, že … 1. Příklad
(
Z lékařských šetření je známo, že výška chlapců ve věku 9-10 let má normální rozdělení N µ, σ 2 kde µ = 139,1 a σ = 40,2 . Změřili jsme výšku u 15 chlapců a vypočetli její průměr. S jakou pravděpodobností bude tento průměr větší než 140?
)
2
Řešení: X výška vybraného chlapce X výběrový průměr z výšek
140 − µ P ( X > 140) = P Z > n = P( Z > 0,55) = 0, 291 σ
Excel: a n mi sig2 E[mX] D[mX] P(mX>a)=?
140 15 139.1 40.2 139.1 2.68 0.291241
=B3 =B4/B2 =1-NORMDIST(B1;B5;ODMOCNINA(B6);1)
2. Příklad
(
)
Hmotnost výrobků v gramech má normální rozdělení N µ = 2000, σ 2 = 16 . Jaká je pravděpodobnost, že průměrná hmotnost n náhodně vybraných výrobků bude větší, než 2002g, je-li počet vybraných výrobků: a) 8, b) 12, c) 16 ? Řešení: X hmotnost výrobku X hmotnost výběrového průměru
(
)
2002 − µ n P ( X > 2002) = P Z > n = P Z > σ 2
a)
n = 8, P Z > 2 = 0,079
b)
n = 12, P Z > 3 = 0,042 n = 16, P( Z > 2) = 0,023
c)
(
)
Excel (ad 1) a n mi sig2
2002 8 2000 16
E[mX] D[mX]
2000 2
P(mX>a)=?
0.07865
=B3 =B4/B2
49
=1-NORMDIST(B1;B6;ODMOCNINA(B7);1)
cv02c.doc
Centrální limitní věta Pro velký výběr (tj. n > 30 ) platí: výběrový průměr X normální rozdělení bez ohledu na to, jaké rozdělení má náhodná veličina X . Normování pro výběrový průměr X a výběrový součet S náhodné veličiny X s momenty
E[ X ] = µ , D[ X ] = σ 2 Průměr:
Z= Součet: platí X =
X −µ n σ
S po dosazení a úpravě dostaneme n S − nµ Z= σ n
Normování pro výběrový podíl p a výběrový počet n1 náhodné veličiny X s alternativním rozdělením s podílem π . Momenty jsou E[ X ] = π , D[ X ] = π(1 − π) . Pro velká n lze použít Centrální limitní větu. Platí předchozí vzorečky s korespondencí p ↔ X a n1 ↔ S a dosazením µ = π a σ 2 = π(1 − π) : Podíl:
Z=
p−π π( 1 − π)
n
Počet:
Z=
n1 − nπ nπ( 1 − π )
Náhodná veličina Z má normované normální rozdělení a kritické hodnoty z α nebo kvantily ζ α je možno nalézt v tabulkách.
1. Příklad Pravděpodobnost, že se za určitou dobu porouchá jeden přístroj je 0,2. Jaká je pravděpodobnost, že se za stejnou dobu ze 100 přístrojů porouchá a) alespoň 20; b) méně než 28; c) 14 až 26 přístrojů? Řešení:
0,2.0,8 = 0,0016 100 20 20 / 100 − 0,2 = 0.5 a) P p > = P Z > 100 0.0016 28 28 / 100 − 0, 2 = 0.977 b) P p < = P Z < 100 0.0016 E[ p ] = π = 0,2; D[ p] =
c)
50 14 26 P p ∈ , = 1 − P( p < 0.14) − P( p > 0.26 ) = 0.8664 100 100
cv03a.doc
Exponenciální třída rozdělení f ( x, θ) = exp{ Q ( θ )U ( x ) + R( θ) + V ( x ) } ;
{ x | f ( x, θ) ≠ 0}
nezávisí na θ
Alternativní rozdělení
f ( x, π) = π x (1 − π)
1− x
π = exp x log + log( 1 − π ) 1 − π
Poissonovo rozdělení
f ( x, λ ) = e =λ
λx = exp{ x log( λ ) − λ − log( x!) } x!
Normální rozdělení
(
)
f x,[µ, σ 2 ] =
1 x − µ 2 1 exp− = 2 πσ 2 σ
1 [ x 2 , x] − 2 2σ = exp µ σ 2
1 µ2 − 2 + log σ 2 2σ
( )
1 − log( 2π) 2
Exponenciální rozdělení
−1 1 x f ( x, δ) = exp− = exp x − log( δ) δ δ δ
51
cv03b.doc
Vlastnosti bodových odhadů T statistika, θ parametr
Nestrannost
E[T ] = θ
Konzistence
lim P ( Tn − θ > ε ) = 0 → lim D[Tn ] = 0
n →∞
n→∞
Vydatnost
MSE( T2 ) D[T2 ] → ( ) MSE T1 D[T1 ]
V (T1 , T2 ) =
1. Příklad Je dán výběr o rozsahu n z náhodné veličiny X s pravděpodobnostní funkcí
f ( x ) = π x (1 − π )
1− x
, x ∈ { 0, 1} . Ukažte, že statistika T = p =
odhadem parametru π .
1 n
n
∑X
i
je nestranným a konzistentním
i=1
Řešení:
E[ X ] =
∑ xπ
x
(1 − π) 1−x
2 E[ X 2 ] = π , ⇒ D[ X ] = E[ X 2 ] − ( E[ X ]) = π(1 − π )
=π,
Nestrannost
E[T ] = E[ Konzistence
1 n
n
∑
Xi ] =
i =1
1 n
n
∑ E[ X
i
]=π
i =1
D[ X ] π( 1 − π ) = lim =0 2 n → ∞ nε n →∞ nε 2
lim P( p − π > ε) ≤ lim
n →∞
Cvičení: Totéž pro T = n1 =
n
∑X
i
. (Je nestanný, není konzistentní).
i =1
2. Příklad x
1 − Je dán výběr o rozsahu n z náhodné veličiny X s hustou pravděpodobnosti f ( x ) = e λ , x > 0 . λ
Ukažte, že statistika T = X =
1 n
n
∑X
i
je nestranným a konzistentním odhadem parametru λ .
i =1
Řešení:
E[ X ] =
∫
∞
0
x
1 − 2 2 2 x e λ dx = per partes = λ , E[ X ] = 2 krát per partrs = 2λ , ⇒ D[ X ] = λ λ
Nestrannost: Konzistence:
52 E[T ] = E[ X ] = E[ X ] = λ
cv04.doc
Konstrukce bodových odhadů Metoda maximální věrohodnosti
Náhodná veličina X má hustotu pravděpodobnosti f ( x, θ ) , závislou na parametru θ . Provedeme výběr x = ( x1 , x 2 , ..., x n ) a odhadujeme parametr θ . Předpoklad: exponenciální třída
f ( x, θ) = exp{ Q( θ)U ( x ) + R ( θ) + V ( x )}
Q ' S + nR ' = 0 ,
Podmínka
kde
S=
n
∑U ( x ) i
i =1
Q '' S + nR '' < 0
Ověření
Metoda momentů Porovnáme obecné (nebo centrální) momenty souboru a výběru
E[ X k ] =
1
n
∑( x ) n i
k
, k = 1, 2, 3, ...
i =1
1. Příklad
Proveďte konstrukci odhadu parametru π náhodné veličiny s alternativním rozdělením Alt ( π ) . Řešení:
Maximální věrohodnost Exponenciální třída
f ( x ; θ ) = π x (1 − π )
π = exp x log + log(1 − π ) 1 − π
1− x
Odtud je
S=
U =x
n
∑x
i
i =1
π Q = log 1− π
Q' =
1 π(1 − π)
R = log(1 − π )
R' =
−1 1− π
Q '' =
2π − 1
π ( 1 − π) −1 R '' = (1 − π) 2 2
2
Podmínka
1 π( 1 − π ) ⇒
∑x
πˆ =
Ověření
2π − 1
π (1 − π) 2
2
∑x
i
−n
1
(1 − π)
2
< 0 ⇒ π < 1 OK
i
1 n
−n
∑x
53
i
1 =0 1− π =x=p
(dosazeno
∑x
i
= nπ )
cv05.doc
Intervaly spolehlivosti T odhadová statistika (funkce výběru, odhadující určitý parametr rozdělení) θ odhadovaný parametr (jeho bodový odhad je hodnota T spočtená z realizace výběru) Intervalový odhad je I θ = ( θ d , θ h ) pro který platí α P ( T < θ d ) = P( T > θ h ) = ⇒ P( T ∈ I θ ) = 1 − α 2 IS pro µ (známe σ 2 )
Iµ = x ± IS pro µ (neznáme σ 2 )
Iµ = x ± IS pro podíl π
n s n
zα / 2
tα / 2
p( 1 − p ) zα / 2 n
Iπ = p ± IS pro σ 2
σ
( n − 1) s 2 ( n − 1) s 2 I σ2 = 2 , 2 χ1− α / 2 χα / 2
1. Příklad (ještě pravděpodobnost že..) Předpokládáme, že velký ročník na VŠ má výsledky testu rozděleny normálně, se střední hodnotou 72 bodů a směrodatnou odchylkou 9 bodů. a) Určete pravděpodobnost, že náhodně vybraný student bude mít výsledek testu nad 80 bodů. b) Provedeme výběr o rozsahu n = 4 a ještě jeden pro n = 9 . Jaká je pravděpodobnost, že výběrový průměr z těchto výběrů bude větší než 80? Řešení: a)
80 − 72 P ( X > 80) = P Z > = 0,187 9
b) pro n = 4
80 − 72 P ( X > 80) = P Z > 4 = 0,038 9
pro n = 9
Excel a mi sig n P(X>a)=? P(mX>a)=?
80 − 72 P ( X > 80) = P Z > 9 = 0,004 9 80 72 9 4 0.187031 0.03772
54 =1-NORMDIST(B1;B2;B3;1) =1-NORMDIST(B1;B2;B3/ODMOCNINA(B4);1)
cv06.doc
Testy hypotéz Schéma: 1) Známe: veličiny, které známe nebo vypočteme 2) Testová statistika: použitá normovaná statistika a její rozdělení 3) Hodnota statistiky: dosazen naměřený výběr 4) Kritický obor: interval hodnot statistiky pro zamítnutí 0-hypotézy 5) Závěr: patří/nepatří do kr. oboru + slovní odpověď
1. Příklad (TH pro střední hodnotu) Byl proveden výběr o rozsahu n = 16 z normálního rozdělení a byl vypočten výběrový průměr X = 4518 a výběrová směrodatná odchylka s = 115 . Testujte hypotézu H 0 : µ = µ 0 = 4500 proti alternativě H A : µ > µ 0 na hladině významnosti α = 0,05 . Řešení: 1) Známe
H 0 : µ 0 = 4500; H A : µ > µ 0 ,
2) Testová statistika
T=
X − µ0 s
T=
4518 − 4500 4 = 0,626 115
3) Hodnota statistiky
4) Kritický obor
α = 0,05
n ≈ St( n − 1)
W = ( t 0, 05 (15) ; ∞ ) = (1,75; ∞ )
5) Závěr
T ∉ W ⇒ H 0 nezamítáme. Data nepřináší dostatek evidence pro zamítnutí nulové hypotézy. Excel mi0 mX s n al
4500 4518 115 16 0.05
T 0.626087 W 1.753051 až inf Hypotézu H0 nezamítáme
2. Příklad (pro střední hodnotu) Provedli jsme náhodný výběr o rozsahu 10 z normálního rozdělení se směrodatnou odchylkou 3 a vypočetli průměr 25,25. Na hladině 0,05 testujme nulovou hypotézu µ 0 = 24 proti alternativě a) µ ≠ µ 0 , b) µ > µ 0 , c) µ < µ 0 . Řešení: 1) n = 10; X = 25,25; σ = 3; α = 0,05 2)
T=
X − µ0 σ
n ≈ N ( 0;1)
55
cv07.doc
IS a TH pro dvě náhodné veličiny 1. Příklad (TH pro dvě střední hodnoty) Zjišťoval se vliv dvou katalyzátorů na konverzi plynu (v procentech). Výsledky jsou uvedeny v tabulce 1. druh 2. druh
24.1 23.2
22.8 24.5
24.4 24.1
24.7 25.1
23.9 22.6
24.2
24
Ověřte hypotézu, že výsledky obou katalyzátorů jsou na hladině významnosti 0,05 stejné. Řešení:
H 0 : µ1 = µ 0 ; H A : µ1 ≠ µ 2 , α = 0,05 , předpoklad σ1 = σ 2
X1 = X 2 T= ≈ St ( n1 + n2 − 2) ; 1 1 sp + n1 n 2 1. Druh 2. Druh
24,1 23,2
mX1 mX2 s1 s2 n1 n2 al
23,98 23,95714 0,725948 0,826352 5 7 0,05
ný sp T W
10 0,787727 0,049555 -inf H0 se
22,8 24,5
sp =
24,4 24,1
( n1 − 1) s12 + ( n2 − 1) s 22 n1 + n 2 − 2
24,7 25,1
23,9 22,6
24,2
24
=PRŮMĚR(B1:F1) =PRŮMĚR(B2:H2) =SMODCH.VÝBĚR(B1:F1) =SMODCH.VÝBĚR(B2:H2) =POČET(B1:F1) =POČET(B2:H2)
=B8+B9-2 =ODMOCNINA(((B8-1)*B6^2+(B9-1)*B7^2)/B12) =(B4-B5)/B13/ODMOCNINA(1/B8+1/B9) až -2,228 U 2,228 až inf =TINV(B10;B12) nezamítá =KDYŽ(NEBO(B14F15);"zamítá";"nezamítá")
2. Příklad (TH pro 2 střední hodnoty – párové výběry) Dvěma laboratorními metodami se zjišťoval obsah chemické látky v roztoku (v %). Předpokládáme, že výsledky mají normální rozdělení. Bylo provedeno měření pěti vzorků (každý vzorek byl měřen oběma metodami) s výsledky 1. metoda 2. metoda
2.3 2.4
1.9 2.0
2.1 2.0
2.4 2.3
2.6 2.5
Ověřte, zda existuje na hladině významnosti 0,05 podstatný rozdíl mezi oběma metodami. Řešení:
H 0 : µ1 = µ 2 ; H A : µ1 ≠ µ 2 , α = 0,05
D T= sD
56
n ≈ St ( n − 1)
