NUMERIČNE METODE V VARNOSTI

Univerza v Ljubljani
FAKULTETA ZA KEMIJO IN
KEMIJSKO TEHNOLOGIJO
Oddelek za tehniško varnost
NUMERIČNE METODE V VARNOSTI
zapiski predavanj
Jože Šrekl
Ljubljana 2013
J. Šrekl
Numerične metode v varnosti
2
Vsebina
Uvod ..................................................................................................................................................... 3
Navadne diferencialne enačbe (ODE) .............................................................................................. 5
1.1.
Diferencialne enačbe prvega reda ........................................................................................... 5
1.2.
Eulerjeva metoda..................................................................................................................... 6
1.3.
Izboljšana Eulerjeva metoda ................................................................................................. 10
1.4.
Metoda Runge-Kutta ............................................................................................................. 14
1.5.
Sistem ODE .......................................................................................................................... 16
2
Fouriereva analiza .......................................................................................................................... 19
2.1.
Razvoj v trigonometrijsko vrsto ............................................................................................ 19
2.2.
Splošni interval ..................................................................................................................... 22
3
Transportna enačba ......................................................................................................................... 25
3.1.
Teoretične osnove ................................................................................................................. 25
3.2.
Analitična rešitev .................................................................................................................. 26
3.2.
Numerično reševanje parabolične PDE ................................................................................. 28
4
Regresija ......................................................................................................................................... 31
4.1.
Korelacije .............................................................................................................................. 31
4.2.
Empirični model .................................................................................................................... 32
4.3.
Regresijska premica .............................................................................................................. 34
4.4.
Nelinearna regresija .............................................................................................................. 42
4.5.
Test hipoteze o enostavni linearni regresiji ........................................................................... 46
4.6.
Analiza variance za regresijo ................................................................................................ 48
4.7.
Multipla regresija .................................................................................................................. 54
5
Analiza variance za več obravnav .................................................................................................. 57
5.1.
Analiza variance .................................................................................................................... 57
5.2.
Razlike med karakterističnimi skupinami ............................................................................. 62
6
Modeli strukturnih enačb ................................................................................................................ 67
6.1.
Merjene in latentne spremenljivke ........................................................................................ 67
6.2.
Konstrukcija modelov regresij [5] ........................................................................................ 68
Bivariantna regresija ....................................................................................................................... 68
Multiple regresije............................................................................................................................ 68
6.3.
Strukturni in merski model [5] .............................................................................................. 70
Strukturni model (SEM) ................................................................................................................ 70
Merski model .................................................................................................................................. 70
6.4.
Skladnost statistike ................................................................................................................ 86
7
Literatura ........................................................................................................................................ 89
Razlogi za pisanje in objavo zapiskov .................................................................................................... 90
1
J. Šrekl
Numerične metode v varnosti
3
Uvod
V naravoslovju in tehniki se pogosto srečujemo s pojavi ki jih želimo opisati z
matematičnimi izrazi, ali pa vsaj izračunati vrednosti znanih ali neznanih izrazov.
Tudi področje varnosti ni izjema zlasti, ko želimo na tem področju napraviti nekaj več
kot rutinsko delo. Raziskave na področju tehnike, pa tudi na področju družboslovja so
običajno sestavljene iz dveh delov. Na eni strani poskušamo z eksperimentalnim
delom (eksperiment in meritve v naravoslovju in tehniki in meritve z anketami,
opazovanjem obnašanja, uporabo statističnih podatkov v družboslovju) ugotoviti
zakonitosti pojavov, na drugi strani pa jih poskušamo z matematičnimi modeli
spremeniti v splošna pravila in zakonitosti. Zakonitosti so lahko že razvite in
dokazane, pogosto jih je le potrebno povezati v novo zakonitost. Padanje prašnega
delca je odvisno od Newtonovega zakona, zakonitosti zračnega upora, gravitacijske
privlačnosti in morda še česa. Vse to združimo v eno samo enačbo. Same enačbe
pogosto ne omogočajo direktnega izračuna vrednosti ob danih parametrih in danih
vrednostih neodvisne spremenljivke kot so čas, krajevne vrednosti in podobno.
Za izdelavo matematičnega modela za opis pojava je potrebno nekaj standardnih
korakov, da dosežemo rezultat:
1. Določimo neodvisne in odvisne spremenljivke v sistemu pojava, ki ga
obravnavamo (čas je vedno neodvisna spremenljivka).
2. Izberemo enote in mere za posamezne spremenljivke (enote morajo biti
usklajene).
3. Določimo osnovne principe, ki opisujejo dogajanje (fizikalni zakoni itd).
4. V teh osnovnih principih (zakonih) uporabimo naše spremenljivke - to je
običajno najtežji del. S te izdelamo matematične enačbe, ki popisujejo pojav.
5. Prepričati se moramo, da imajo vsi deli enačbe enako fizikalno enoto. Ne
moremo seštevati različnih enot. Pogosto uporabljamo izraze brez (fizikalni ali
kakih drugih) enot.
6. Iz točke 4. dobimo enačbo, ki je matematični model dogajanja.
7. V primerih, ko je dogajanje kompleksno – sestavljeno, dobimo sistem večih
enačb. Primer: Newtonova enačba v prostoru .
8. Pri kompliciranih enačbah ocenimo dejanski prispevek posameznih delov
enačbe in zanemarimo dele z zanemarljivim prispevkom. (Inženirski pristop)
Model nam da splošno povezavo med neodvisnimi in odvisnimi količinami –
spremenljivkami. V praksi pa nas zanimajo konkretne rešitve pri danih vrednostih
spremenljivk in parametrov, ki jih lahko izbiramo ali merimo. Za iskanje teh rešitev
obstaja več poti. Najpogostejše so:
 Analitično reševanje do eksplicitne odvisnosti odvisne (iskane) vrednosti od
neodvisnih (merljivih ali določljivih) vrednosti.
 Numerično reševanje, ki iz kakršne koli oblike zveze pripelje do bolj ali manj
natančnih vrednosti iskanih spremenljivk.
 Eksperimentalno merjenje neodvisnih in odvisnih vrednosti, iz teh rezultatov
določimo približne vrednosti odvisnih vrednosti pri poljubni določitvi
neodvisnih.
Seveda pa so lahko vse metode pomešane in v različnih segmentih raziskave
uporabljamo različne metode (eno samo ali več hkrati).
J. Šrekl
Numerične metode v varnosti
4
Tako kot veljajo pravila v analitični matematiki, veljajo tudi pravila v numerični
matematiki. Pri tem je pomembno kakšne metode uporabljamo.
Poznavanje metod reševanja in njihovi zmožnosti (zahtevnost računanja, natančnost
računanja) je za izbor in uporabo ključnega pomena. Zato je smiselno poznati čim več
metod in medijev na katerih jih je mogoče uporabiti (računalniški programi in
računalniki druga orodja za računanje).
Zato je tudi za področje varnosti in zdravja, požarne varnosti in ekološke varnosti
pomembno poznavanje numeričnih metod zlasti s področje reševanja enačb, reševanja
vseh vrst diferencialnih enačb, sistemov enačb in statističnih rešitev. S temi orodja
lahko učinkovito raziskujemo, iščemo rešitve, potrjujemo eksperimentalne rezultate in
razvijamo nove modele za pojasnjevanje pojavov.
Knjiga ni namenjena razvijanju numeričnih metod in potrjevanju njihove pravilnosti,
ampak je namenjena predvsem uporabnemu učenju metod, ki so koristno orodje pri
raziskovalnem delu.
J. Šrekl
Numerične metode v varnosti
5
Gibanje pomeni model z diferencialno enačbo
1
1.1.
Navadne diferencialne enačbe (ODE)
Diferencialne enačbe prvega reda
Pri problemih, kjer gre za ravnotežje (sil, tlakov, napetosti, itd.), se srečujemo z
navadnimi enačbami, ki popisujejo medsebojno povezavo dveh ali več spremenljivk.
Vendar pa pri naravnih pojavih pogosto srečujemo dogajanje, kjer se nekaj spreminja.
Če se avto premika po cesti, se s spreminjajočim časom spreminja njegova prevožena
pot, morda tudi hitrost in pospešek. V zaprti posodi segrevamo tekočino, s časom se
spreminjata tlak in temperatura tekočine. Takega spreminjajočega sistemu ni mogoče
popisati z običajnimi linearnimi ali nelinearnimi enačbami, ker se poleg samih
spremenljivk (neodvisnih in odvisne) pojavijo tudi hitrosti spreminjanja ki se izražajo
kot odvodi odvisne spremenljivke (sistem postane dinamičen) in dobimo diferencialne
enačbe.
Odvod v neki točki določa smer spreminjanja funkcijske zveze. Splošna enačba smeri
povezuje odvod z neodvisno in odvisno spremenljivko, ki sta povezani z neko
funkcijsko zvezo. Grafično to lahko prikažemo v koordinatnem sistemu, kjer ima
vsaka točka v ravnini določeno smer glede na svojo lego.
dy
 f ( x, y )
dx
Slika 1.1 Smer, ki jo določa ODE v dani točki
Analitično reševanje diferencialne enačbe  ′ = (, ) poteka s pomočjo integracije
(kvadratur). Rešitev je eno parametrična družina krivulj  = (, ), kjer je c
neznana konstanta. Če poznamo vsaj eno vrednost na krivulji, običajno jo imenujemo
začetna vrednost, lahko določimo konkretno krivuljo rešitev.
J. Šrekl
Numerične metode v varnosti
6
Primer 1.1: Če je za enačbo
m
( p x)
l
ki popisuje mešanje koncentracije neke snovi v tekočino, rešitev enaka p-x = Ce-(m/l)t,
potem je končna rešitev ob upoštevanju začetnega pogoja x(0) = a enaka
x = p + (a – p)e-(m/l)t .
x 
ODE lahko rešujemo analitično ali numerično. Pri analitičnem reševanju v naslednje
koraku, ko iščemo posamezne vrednosti, pogosto moramo uporabiti numerične
metode (približno računanje funkcijskih vrednosti – to nam običajno počne kalkulator
ali računalnik, ne da bi se tega sploh zavedali). Reševanje nekaterih enačb je kar
težavno. Zato se bomo izogibali analitičnemu reševanju, če pa že dobimo analitično
rešitev, si na koncu pogosto moramo pomagati z numeričnimi približki vrednosti
funkcije. Zato bomo že v začetku izbrali numerično pot, ki bo ob dovolj veliki
natančnosti računanja zadoščala »inženirskim rezultatom«. Reševali bomo torej
začetni problem 1. reda. Enačba skupaj z začetnim pogojem ima obliko:
 ′ = (, ) in (0 ) = 0
S h bomo označili korak računanja (kar pomeni gostoto izračunanih vrednosti ali
pomik po abscisni osi za naslednjo rešitev). Rešitev bomo tabelirali v točkah
1 , 2 , 3 , . . .. pri čemer je  +1 =  + ℎ, za n = 0,1,2,.... S pomočjo ustrezne
formule bomo računali zaporedoma približke za vrednosti funkcije f (xi), označili jih
bomo z 1 , 2 , . . . .. Pravimo, da je metoda reda p, če je lokalna napaka oblike ℎ+1.
p pomeni velikostni red napake pri računanju yn+1 iz yn. Oceno za lokalno napako
metode reda p v točki x lahko izračunamo s pomočjo dveh korakov računanja.
Vzemimo za korak računanja najprej h potem pa še 2h.
ℎ =  () + 2ℎ+1
2ℎ =  () + (2ℎ)+1
Pri tem upoštevamo, da smo v prvi enačbi imeli dvakrat daljši korak, zato je tudi
napaka dvojna. Če označimo  = 2ℎ+1 , smemo zgornji dve enačbi zapisati:
ℎ =  () + 
2ℎ =  () + 2 
Iz sistema dveh enačb lahko izračunamo vrednost
 −
 = 22ℎ −1ℎ ,
ki jo uporabimo kot približek za lokalno napako.
1.2.
Eulerjeva metoda
Najpreprostejša metoda za reševanje ODE je Eulerjeva tangentna metoda. Že samo
ime nam pove, da je to linearna metoda, kjer je p=1. Algoritem metode dobimo tako,
da odvod zamenjamo s približkom diference. (To pomeni, da smerni koeficient
J. Šrekl
Numerične metode v varnosti
7
tangente v vsaki točki delitve, zamenjamo s smernim koeficientom sekante na tem
odseku.)
′ =
+1 −
ℎ
, ℎ = +1 −  , n=0,1,2,…
dobimo algoritem, ki nam po korakih računa vrednosti  . Kot smo že zapisali je
metoda reda 1, torej  = 2ℎ − ℎ .
Eulerjevo metode zapišimo še v splošnem računalniškem jeziku (ki ga lahko
prevedemo v poljuben programski jezik):
1.
2.
3.
4.
5.
6.
korak:
korak:
korak:
korak:
korak:
korak:
define f(x,y)
input začetne vrednosti x0 in y0
input velikost koraka h in število korakov n
output x0 in y0
for j od 1 do n do
k1 = f(x,y)
y = y + h*k1
x=x+h
output x in y
end
7. korak:
8. korak:
Primer 1.2: Rešimo začetni problem, ki je popisan z enačbo in pogojem:
 ′ = 1 −  + 4
(0) = 1
Analitična rešitev enačbe ima obliko:
1
3
19
 = 4 − 16 + 16 4
Rešujmo problem najprej grafično:
20
18
16
14
12
10
8
6
4
2
0
0
0,2
0,4
0,6
0,8
1
Iz začetne vrednosti (0,1) potegnemo tangento s smernim koeficientom
 = 1 −  + 4
J. Šrekl
Numerične metode v varnosti
8
za x=0 in y=1. Dobimo k=5. V točki x=0,1 dobimo novo vrednost za y na tangenti:
 =  + 1 = 5 ∙ 0,1 + 1 = 1,5
Izračunamo novi k v tej točki (0,1 ; 1,5) sestavimo novo linearno funkcijo, ki določi
novo tangento in določimo novo vrednost y v x=0,2. Nadaljujemo in dobimo lomljeno
črto, ki je približek rešitve.
Reševanje s kalkulatorjem:
1. korak: definiramo si funkcijo za izračunavanje (, ) = 1 −  + 4
2. korak: 0 = 0 in 0 = 1
3. korak: določimo korak ℎ = 0,01
4. korak: računanje vrednosti:
1 = 1 − 0 + 4 ∗ 1 = 5
 = 1 + 0,01 ∗ 5 = 1,05
 = 0 + 0,01 = 0,01
1 = 1 − 0,01 + 4 ∗ 1,05 = 5,19
 = 1,05 + 0,01 ∗ 5,19 = 1,1019
 = 0,01 + 0,01 = 0,02
Nadaljujemo do vrednosti  = 2,0
Tabela 1.1. Rezultati računanja naloge pri različnih vrednostih h in primerjava z
vrednostmi izračunanimi iz analitične rešitve. (Vir [1])
Primerjave vrednosti za različne h kažejo, da napaka računanja proti končni
vrednosti x narašča seveda odvisno od koraka. Vzemimo rezultate iz tabele in
izračunajmo oceno napake  za ℎ = 0,025 pri vrednosti  = 0,1
 = 1,5761188 − 1,5475000 = 0,0286188
Prava napaka pa je približno dvakrat večja:
 = 1,6090418 − 1,5475000 = 0,0515418
Če pogledamo še relativno napako na tem koraku, dobimo 3,2%. Pri koraku ℎ =
0,001 dobimo relativno napako 0,09%.
Primer 1.3: Določi čas iztekanja nevarne tekočine iz prevrnjene cisterne, odprtina je
pod cisterno, pretok je prost, cisterna je na začetku polna.
Podatki: Cisterna ima eliptično obliko z osema  = 40  (višina 80 cm) ,  =
65  (širina 130 cm) in dolžino  = 520 . Zaobljenost na koncih cisterne
zanemarimo.
J. Šrekl
Numerične metode v varnosti
9
Slika 1.2. Prevrnjena cisterna
Označimo z y višino tekočine v prevrnjeni cisterni in s spremenljivko t čas
iztekanja. Torricelli-jev zakon določa hitrost iztekanja v odvisnosti od višine tekočine
nad iztekom. () = √2 (μ je eksperimentalno določena vrednost, ki je odvisna
od oblike odprtine, vzeli bomo  = 1, g je gravitacijski pospešek). Sprememba
volumna tekočine v cisterni je količina iztečene tekočine in je odvisna od ploščine
odprtine () in hitrosti iztekanja (v). (Če je premer okrogle odprtine 2cm, je  =
0,00126 2).
∆ =  = −()
Če združimo obe enačbi in upoštevamo geometrijo posode, dobimo diferencialno
enačbo:
′ =

1
=−

√(2 − )
V konstanto K smo združili vse konstantne parametre.
2
=
√2
Izračunali bomo s pomočjo EXCEL-a:
Tabela 1.2. Vstavljeni podatki in uporabljene formule
Prevrnjena
cisterna
a=
b=
L=
g=
ω=
K=
0,4
0,65
5,2
9,81
0,00126
=(2*B2*B4)/(B3*B6*SQRT(2*B5))
t
0
=A10+$E$2
y
1,29
=B10+C10*$E$2
h= 2
k
=-1/($B$7*SQRT(2*$B$3-B10))
=-1/($B$7*SQRT(2*$B$3-B11))
J. Šrekl
Numerične metode v varnosti
10
Tabela 1.3. Začetna, nekaj vmesni in končna vrednost (čas je v sek, korak je 5 sek)
t
0
5
60
600
1120
1125
y
1,29
1,246398
1,104274
0,441545
0,003495
-0,00033
k
-0,00872
-0,00377
-0,00197
-0,00094
-0,00077
-0,00076
18,75
min
4247 litrov tekočine izteče v manj kot 19 minutah.
Tabela 1.4. Začetna, nekaj vmesni in končna vrednost (čas je v sek, korak je 2 sek)
t
y
0
2
60
600
1120
1130
k
1,29
1,272559
1,111668
0,446086
0,007424
-0,00024
-0,008720474
-0,005264299
-0,002009453
-0,000943698
-0,000767029
-0,000764766 18,8333 min
1,4
1,2
1
0,8
0,6
0,4
0,2
0
0
200
400
600
800
1000
1200
-0,2
Slika 1.3. Grafični prikaz rešitev
1.3.
Izboljšana Eulerjeva metoda
Eulerjevo metodo za reševanje ODE lahko nekoliko popravimo in s tem izboljšamo
red. Za n=0,1,2,… dobimo algoritem:
J. Šrekl
Numerične metode v varnosti
11
1 = ℎ ( ,  )
ℎ

2 = ℎ ( + 2 ,  + 21 )
 = −1 + 2
Rekli smo, da je metoda drugega reda, kar pomeni je p=2. Napako računamo s
 −
formulo  = 2ℎ3 ℎ
Izboljšano Eulerjevo metode zapišimo še v splošnem računalniškem jeziku (ki ga
lahko prevedemo v poljuben programski jezik):
1.
2.
3.
4.
5.
6.
korak:
korak:
korak:
korak:
korak:
korak:
7. korak:
8. korak:
define f(x,y)
input začetne vrednosti x0 in y0
input velikost koraka h in število korakov n
output x0 in y0
for j od 1 do n do
k1 = h* f(x,y)
k2 = h* f(x+h/2,y+k1/2)
y = y + k2
x=x+h
output x in y
end
Primer 1.4: Reši začetni problem iz primera 1.2 z izboljšano Eulerjevo metodo s
pomočjo kalkulatorja. Enačba in začetni pogoj:
 ′ = 1 −  + 4
(0) = 1
Algoritem za reševanje s kalkulatorjem:
1. korak: definiramo si funkcijo za izračunavanje (, ) = 1 −  + 4
2. korak: 0 = 0 in 0 = 1
3. korak: določimo korak ℎ = 0,01
4. korak: računanje vrednosti:
1 = (1 − 0 + 4 ∗ 1) ∗ 0,01 = 0,05
2 = (1 − (0 + 0,005) + 4 ∗ (1 + 0,025)) ∗ 0,01 = 0,05095
 = 1 + 0,05095 = 1,05095
 = 0 + 0,01 = 0,01
1 = (1 − 0,01 + 4 ∗ 1,05095) ∗ 0,01 = 0,051938
2 = (1 − (0 + 0,005) + 4 ∗ (1 + 0,025969)) ∗ 0,01 = 0,05098876
 = 1,05095 + 0,05098876 = 1,10193876
 = 00,01 + 0,01 = 0,02
Nadaljujemo do vrednosti  = 2,0 Vidimo, da dobimo boljše vrednosti kot pri
običajni Eulerjevi metodi. Primerjamo nekaj rezultatov obeh metod z
analitičnimi rešitvami.
J. Šrekl
Numerične metode v varnosti
12
Tabela 1.5. Rezultati računanja naloge pri različnih vrednostih h za obe metodi in
primerjava z vrednostmi analitične rešitve. (Vir [1])
Primer 1.5: Enačba o mešanju koncentracij. Bazen,
mešalo, pritok in odtok. Prostornina bazena je l =
20 m3. Doteka p =2 % raztopina česarkoli
(spreminjanje
prostornine
pri
raztapljanju
zanemarimo), mešalec meša, meša; odteka pa prav
toliko raztopine, kolikor je priteka, denimo m =3,5
litrov v minuti. V začetku (t = 0) je v posodi a =
0,01% raztopina. Kolikšna je koncentracija po t = 8
urah? Kdaj moramo ustaviti iztekanje, če je
dovoljena koncentracija pri izteku 1,2%
V začetku reševanja se ozrimo po neznanki. Kar
koncentracija tistega česarkoli naj bo. Še krstimo jo,
neznanko, y ji recimo. V dovolj majhnem času dt se
koncentracija spremeni za dy. Pomislimo, kaj vse je Slika 1.4. Mešanje koncentracij
povzročilo to spremembo, pa bo enačba pred nami. V
času dt je v posodo priteklo in iz nje odteklo mdt litrov tekočine. Vsak liter, ki je
pritekel, je vseboval pmdt/100 litrov snovi, vsak liter, ki je odtekel pa ymdt/l00 litrov
raztopljene snovi. Zato:

 ′ = ( − )

(0) = 
Analitična rešitev problema je:
 =  + ( – ) −(/)
Tabela 1.6. Vstavimo podatke in EXCEL-ove formule
enačba o mešanju
l=
20000
l
m=
3,5
l/min
p=
0,02
a=
0,0001
J. Šrekl
Numerične metode v varnosti
13
h=
10
min
t
y
k1
k2
0
0,0001
=$B$6*$B$3/$B$2*($B$4-B9)
=$B$6*$B$3/$B$2*($B$4-B9-C9/2)
=A9+$B$6
=B9+D9
=$B$6*$B$3/$B$2*($B$4-B10)
=$B$6*$B$3/$B$2*($B$4-B10-C10/2)
Tabela 1.7. Rezultat pri koraku 10 min. za 8 ur
enačba o
mešanju
l=
m=
p=
a=
h=
t
20000 l
3,5 l/min
0,02
0,0001
10 min
y
0
10
60
120
180
240
300
360
420
480
0,0001
0,000134795
0,000307857
0,000513542
0,00071708
0,000918491
0,001117799
0,001315024
0,00151019
0,001703317
k1
k2
0,000034825
3,47641E-05
3,44613E-05
3,41013E-05
3,37451E-05
3,33926E-05
3,30439E-05
3,26987E-05
3,23572E-05
3,20192E-05
3,47945E-05
3,47337E-05
3,44311E-05
3,40715E-05
3,37156E-05
3,33634E-05
3,30149E-05
3,26701E-05
3,23289E-05
3,19912E-05
Tabela 1.8. Dovoljeno koncentracijo bomo presegli po 86 urah, kar pomeni 3,58
triizmenskih dni, 5,38 dvoizmenskih dni ali 10,75 enoizmenskih dni.
t
y
0
0,0001
5160 0,011933
5220 0,012018
k1
k2
0,00020895 0,000207853
8,46997E-05 8,4255E-05
8,3815E-05 8,3375E-05
Iz grafa naraščanja koncentracije v bazenu vidimo, da koncentracija narašča skoraj
linearno.
J. Šrekl
Numerične metode v varnosti
14
Naraščanje koncentracije po NM
0,2
0,15
0,1
0,05
0
0
1
2
3
4
5
6
7
8
9
Slika 1. 5. Naraščanje koncentracije v % v odvisnosti od časa v urah (linearni graf po
Newtonovi metodi)
1.4.
Metoda Runge-Kutta
Izboljšano Eulerjevo metodo za reševanje NDE popravimo tako da dvignemo red na
p=4. Metodo poznamo pod imenom Runge-Kutta. Za n=0,1,2,… dobimo algoritem:
1 = ℎ ( ,  )
ℎ

2 = ℎ ( + 2 ,  + 21 )
ℎ

3 = ℎ ( + 2 ,  + 22 )
4 = ℎ ( + ℎ ,  + 3 )
1
 = −1 + 6(1 + 22 + 23 + 4 )
Rekli smo, da je metoda četrtega reda, kar pomeni je p=4. Napako računamo s
 −
formulo  = 2ℎ15 ℎ
Izboljšano Eulerjevo metode zapišimo še v splošnem računalniškem jeziku (ki ga
lahko prevedemo v poljuben programski jezik):
1.
2.
3.
4.
5.
6.
korak:
korak:
korak:
korak:
korak:
korak:
7. korak:
8. korak:
define f(x,y)
input začetne vrednosti x0 in y0
input velikost koraka h in število korakov n
output x0 in y0
for j od 1 do n do
k1 = h* f(x,y)
k2 = h* f(x+h/2,y+k1/2)
k3= h* f(x+h/2,y+k2/2)
k4 = h* f(x+h,y+k4)
y = y +(1/6)*( k1+2*k2+2*k3+k4)
x=x+h
output x in y
end
J. Šrekl
Numerične metode v varnosti
15
V primeru, ko funkcija f ni odvisna od y lahko zapišemo enostavnejše formule za
računanje:
1 = ℎ ( )
ℎ
2 = ℎ ( + 2 )
ℎ
3 = ℎ ( + 2 )
4 = ℎ ( + ℎ )
ℎ
ℎ
 = −1 + 6( ( ) + 4 ( + 2 ) +  ( + ℎ ))
Vidimo, da je za računanje dovolj kar zadnja formula.
Primer 1.6: Izračunaj primer 1.5 z metodo Runge-Kutta.
Tabela 1.9. Del EXCEL-ove tabele za izračun
t
y
0
10
20
30
k1
k2
k3
k4
0,0001 0,000034825 3,47945E-05 3,47946E-05 3,47641E-05
0,000135 3,47641E-05 3,47337E-05 3,47337E-05 3,47033E-05
0,00017 3,47033E-05 3,4673E-05 3,4673E-05 3,46426E-05
0,000204 3,46426E-05 3,46123E-05 3,46124E-05 3,45821E-05
480 0,001703 3,20192E-05 3,19912E-05 3,19912E-05 3,19632E-05
5200 0,01199 1,40179E-05 1,40056E-05 1,40057E-05 1,39934E-05
5210 0,012004 1,39934E-05 1,39812E-05 1,39812E-05 1,39689E-05
Čas dosežene kritične raztopine se razlikuje pri izračunu s popravljeno Newtonovo
metodo in metodo Runge-Kutta le za približno 10 min, kar v odstotkih pomeni
10
100% = 0,2%
5210
Ker smo vzeli majhen korak računanja so seveda rezultati podobni ne glede na
metodo. Vendar pa se pri podrobnem pogledu vidi, da v tem primeru koncentracija ne
narašča linearno.
J. Šrekl
Numerične metode v varnosti
16
Naraščanje koncentracije po R-K
0,014
0,012
0,01
0,008
0,006
0,004
0,002
0
0
1000
2000
3000
4000
5000
6000
Slika 1. 6. Naraščanje koncentracije v % v odvisnosti od časa v urah (nelinearni graf
po Runge-Kutta)
Naraščanje koncentracije analitično
0,014
0,012
0,01
0,008
0,006
0,004
0,002
0
0
1000
2000
3000
4000
5000
6000
Slika 1. 7. Naraščanje koncentracije v % v odvisnosti od časa v urah (nelinearni graf
po analitičnem računanju)
1.5.
Sistem ODE
Pogosto se srečujemo tudi z več diferencialnimi enačbami, ki so povezane v sistem
(spremenljivke se v enačbah prepletajo). Z vektorskimi oznakami zapišemo sistem
enačb:
⃗⃗⃗
 ′ = (, ),
(0 ) = 0
J. Šrekl
Numerične metode v varnosti
17
Za numerično reševanje uporabimo Eulerjevo formulo v vektorski obliki:
+1 =  + ℎ ,
ali pisano po komponentah za dve enačbi:
+1

( ,  ,  )
( ) = (  ) +h(    ).
( ,  ,  )
+1

Lahko pa uporabimo tudi metodo Runge-Kutta:
1 = ℎ ( ,  )
ℎ

2 = ℎ ( + 2 ,  + 21 )
ℎ

3 = ℎ ( + 2 ,  + 22 )
4 = ℎ ( + ℎ ,  + 3 )
1
+ =  + 6(1 + 22 + 23 + 4 )
S poudarjenimi (bold) oznakami so označeni vektorji.
Primer 1.7: Reši sistem enačb
 ′ =  −+ − 
 ′ = sin( − 3)
(0) = 1,
(0) = 2
Reševanja se bomo lotili z metodo Runge-Kutta in sicer v EXCEL-u. Za računanje
bomo vzeli korak ℎ = 0,1 Zapišimo algoritme za računanje:
11 = ℎ ∗ ( (− + ) − cos())
12 = ℎ ∗ (sin( − 3 ∗ ))
21 = ℎ ∗ ((− − 11⁄2 +  + 12⁄2) − ( + 11⁄2))
22 = ℎ ∗ (( + 11/2 − 3 ∗ ( + 12⁄2))
31 = ℎ ∗ ((− − 21⁄2 +  + 22⁄2) − ( + 21⁄2))
32 = ℎ ∗ (( + 21/2 − 3 ∗ ( + 22⁄2))
41 = ℎ ∗ ((− − 31 +  + 32) − ( + 31)
42 = ℎ ∗ (( + 31 − 3 ∗ ( + 32))
( + 1) =  + (1⁄6) ∗ (11 + 2 ∗ 21 + 2 ∗ 31 + 41)
( + 1) =  + (1⁄6) ∗ (12 + 2 ∗ 22 + 2 ∗ 32 + 42)
J. Šrekl
Numerične metode v varnosti
18
Tabela 1.10. Izračun z EXCEL-ovo tabelo
h=
t
0,1
x
0
0,1
0,2
0,3
0,4
0,5
0,6
0,7
0,8
0,9
1
1,1
1,2
1,3
1,4
1,5
1,6
1,7
1,8
1,9
2
1
1,168623
1,321701
1,462956
1,594605
1,717978
1,8339
1,942912
2,045419
2,141756
2,232236
2,317165
2,396846
2,471579
2,541662
2,607381
2,669015
2,726827
2,781069
2,831976
2,879769
y
2
2,091759
2,181761
2,268411
2,350411
2,42687
2,497323
2,561673
2,620096
2,672945
2,72067
2,763756
2,802685
2,837914
2,869861
2,8989
2,925366
2,949549
2,971705
2,992057
3,010797
k11
k12
k21
k22
k31
k32
k41
k42
0,217798 0,095892 0,156347 0,094843 0,163902 0,093879 0,153445 0,077218
0,212575 0,092328 0,137799 0,091043 0,146225 0,089519 0,137841 0,086557
0,211678 0,087216 0,123048 0,085967 0,132319 0,083722 0,12512 0,093307
0,213008 0,080809 0,11067 0,079934 0,120686 0,07684 0,114172 0,09764
0,215313 0,07356 0,099833 0,073378 0,110451 0,069355 0,104359 0,099728
0,217839 0,06598 0,090064 0,066724 0,101116 0,061745 0,09533 0,099801
0,220151 0,058519 0,081101 0,060312 0,092412
0,0544 0,086902 0,098156
0,222021 0,051506 0,072806 0,054369 0,084209 0,047585 0,078986 0,095124
0,223356 0,045138 0,065108 0,049017 0,076454 0,041445 0,071544 0,091032
0,224139 0,039503 0,057962 0,044298 0,06913 0,036035 0,064558 0,08618
0,224402 0,03461 0,051342
0,0402 0,062235 0,031341 0,058018 0,080823
0,224195 0,030418 0,045219 0,036679 0,055769 0,027316 0,051913 0,075168
0,223582 0,026863 0,039568 0,033676 0,049727 0,023891 0,046229 0,069376
0,222625 0,02387 0,034359 0,031127 0,044101 0,020994 0,040949 0,063568
0,221384 0,021361 0,029564 0,02897 0,038876 0,018551 0,036053 0,057835
0,219913 0,019265 0,025151 0,027147 0,034033 0,016496 0,031519 0,05224
0,21826 0,017519 0,021092 0,025607 0,029552 0,01477 0,027324 0,046827
0,216468 0,016067 0,017358 0,024304 0,02541 0,013319 0,023446 0,041625
0,214573 0,014859 0,013921
0,0232 0,021583
0,0121 0,01986 0,03665
0,212606 0,013854 0,010756 0,022263 0,018049 0,011075 0,016545 0,031911
0,210592 0,013019 0,007839 0,021465 0,014786 0,010211 0,01348 0,027412
J. Šrekl
Numerične metode v varnosti
19
Tudi nenehno ponavljanje je lahko lepota.
2 Fouriereva analiza
2.1. Razvoj v trigonometrijsko vrsto
V uporabni matematiki pogosto naletimo na probleme približnega računanja ali
približnega izražanja funkcij. Najbolj prikladen način izražanja funkcije je
izražanje s potenčnimi vrstami
∞
() = ∑   
=1
V vsaki točki konvergenčnega območja je mogoče poljubno natančno izračunati
vrednost funkcije, če seštejemo zadostno število členov vrste. Zapletene
transcendentne funkcije lahko na tak način izražamo s preprostimi algebrajskimi
funkcijami (to je tudi način, ki ga uporablja večina računalnikov).
Zastavimo si preprosto vprašanje: Ali je mogoče zaporedje potenčnih funkcij
zamenjati z neskončnim zaporedjem kakšnih drugačnih preprostih funkcij
1 (), 2 (), 3 (), …, ki jih bi bilo mogoče sestaviti v vrsto?
Vzemimo zaporedje trigonometričnih funkcije na osnovnem intervalu [−π, π] :
1, , , 2, 2, ⋯ , , ⋯
Zaporedje teh funkcij imenujemo trigonometrično zaporedje. Za to zaporedje
poznamo naslednje lastnosti:

2
 (cos kx) dx 






2
 (sin kx) dx 


1 cos 2 x
2
1 cos 2 x
2
dx  
dx  





 (sin nx  cos mx)dx  12
 (sin(n  m) x  sin(n  m) x)dx  0


 (sin nx  sin mx)dx   (cos(n  m) x  cos(n  m) x)dx  0
1
2




 (cos nx  cos mx)dx   (cos(n  m) x  cos(n  m) x)dx  0
1
2


J. Šrekl
Numerične metode v varnosti
20
Enačbe veljajo za poljubna cela števila k, n, m, ki so različna od 0. Zaradi teh
lastnosti pravimo, da je zaporedje trigonometričnih funkcij ortogonalno. Če vse
funkcije še delimo s  , dobimo ortonormirano zaporedje. Vzemimo funkcijo
f(x) definirano na intervalu [−π, π]. Izrazimo funkcijo kot vsoto
f ( x) 

a0
2
  (ak cos kx  bk sin kx)
k 1
Neznane koeficiente v neskončni vsoti izračunamo tako, da vrsto zaporedoma
množimo s trigonometričnimi funkcijami iz zaporedja in vsakokrat produkt
integriramo. Pri integraciji upoštevamo ortogonalnost zaporedja.

 a
 
a0
 f ( x ) dx   0 dx    ( ak cos kx  bk sin kx ) dx  2 
2
k 1 





f ( x) cos nx dx 




k 1 

bk sin kx  cos nx)dx  ak


f ( x)sin nx dx 




a0
2

cos nx dx    (ak cos kx  cos nx 
a0
2

sin nx dx    (ak cos kx  sin nx 
k 1 
bk sin kx  sin nx)dx  bk
Neskončna vrsta, ki jo lahko zapišemo s pomočjo znanih koeficientov, se imenuje
trigonometrijska Fouriereva vrsta ali kratko trigonometrijska vrsta. Koeficiente
0 ,  ,  , , ki nastopajo v tej vrsti, imenujemo Fourierevi koeficienti. Če enačba
f ( x) 

a0
2
  (ak cos kx  bk sin kx)
k 1
velja na intervalu [−π, π], potem pravimo, da je trigonometrijska vrsta
konvergentna na tem intervalu. Seveda bo relacija lahko veljala le pri zveznih
funkcijah. Če bo funkcija f(x) le odsekoma zvezna in odvedljiva, bo namesto
prejšnje enačbe veljala enačba
f ( x  0  f ( x  0)) a 
 2   (ak cos kx  bk sin kx)
2
k 1
0
Primer 2.1: Usmernik je priključen na izmenično napetost U  U 0 sin t . Usmernik
zapira, ko je napetost negativna. Označimo izhodno napetost na usmerniku z U * (t )
in razvrstimo to funkcijo v trigonometrijsko vrsto. Funkcija je definirana na
intervalu   ,   in se zapiše
  t  0
 0,
U * (t )  

U 0 sin t , 0  t   
Izračunamo Fouriereve koeficiente:
J. Šrekl
Numerične metode v varnosti
21

a0    U 0 sin tdt  2U 0
1
0

ak 
1

U
0
sin t cos ktdt   ( k2U2 01) , (k  2,4, )
0
Če je k liho število, velja  = 0. Na podoben način lahko izračunamo:

bk  1  U 0 sin t sin ktdt  U20 , k  1, bk  0, k  1
0
Rezultat je mogoče zapisati v obliki vrste
1
1
U * (t )  U0  U20 sin t  2U0 (1.3
cos 2t  3.5
cos 4t
Enakost velja samo na osnovnem intervalu [−π, π]. Vrsta na desni strani
enakosti pa je definirana za vse realne t in nam popisuje napetost za usmernikom
v poljubnem času. Funkcija je seveda periodična, ker je sestavljena iz samih
periodičnih funkcij.
Zato najpogosteje srečujemo razvoje v trigonometrijske vrste ravno v
elektrotehniki, kjer je treba popisovati vsemogoče periodične napetosti. Pogosto se
ukvarjamo s funkcijami, ki so na intervalu [−π, π] simetrične ali antisimetrične,
oziroma sode ali lihe. Za take funkcije je mogoče razvoj v trigonometrijsko vrsto
zapisati enostavneje.
Poglejmo si funkcijo (), ki je na intervalu [−π, π] soda, kar pomeni
(−) = ().
Vemo, da je produkt sodih funkcij soda funkcija in produkt sode in lihe funkcije
pa liha funkcija. Za integrale potem velja:

a  2  f ( x )dx
0 
0

ak  2  f ( x) cos nx dx
0

bk  2  f ( x)sin nx dx  0
0
Ker je funkcija soda se zapiše s samimi sodimi funkcijami, torej s samimi
kosinusi.
f ( x) 

a0
2
  ak cos kx
k 1
Za liho funkcijo (), za katero velja
(−) = − (),
se razvoj zapiše s samimi lihimi funkcijami, torej s sinusi.

f ( x)   bk sin kx
k 1
J. Šrekl
Numerične metode v varnosti
Koeficiente izračunamo z enačbami:
22

a  2  f ( x )dx  0
0 
0

ak    f ( x)cos nx dx  0
2
0

 f ( x)sin nx dx
bk  
2
0
2.2. Splošni interval
Funkcij ne bomo vedno razvijali v trigonometrijsko vrsto v osnovnem intervalu,
ampak jo bomo poskušali razviti tudi v nekem splošnem intervalu, na katerem je
funkcija definirana. Osnovna ideja razvoja v trigonometrijske vrste je spremeniti
funkcijo, definirano na nekem intervalu, v periodično funkcijo, ki ima za periodo
ta interval.
Za doseganje tega cilja si poiščimo transformacijo, ki bo splošni interval [−d, d]
preslikala na osnovni interval.
x
T :x
d
Razvoj v vrsto zapišemo po uvedbi transformacije:

f ( x)  a2   (ak cos kd x  bk sin kd x )
0
k 1
Koeficiente izračunamo z enačbami:
d
a0  1  f ( x )dx
d d
d
ak 
1
d
 f ( x)cos
n x
d
dx
 f ( x)sin
n x
d
dx
d
d
bk  d1
d
Vrnimo se k osnovnemu intervalu   ,   . Funkcije je na tem intervalu mogoče
razviti v trigonometrijsko vrsto, ker so trigonometrijske funkcije ortogonalne.
Ortogonalnost funkcij velja za vsak interval  a, a  2  , saj so vse funkcije
periodične na intervalu   ,   , ali pa na kateremkoli intervalu dolžine 2 .
J. Šrekl
Numerične metode v varnosti
23
Če združimo obe posplošitvi lahko zapišemo razvoj v vrsto na splošnem intervalu
 a, b . Funkcija mora biti na tem intervalu odsekoma zvezna in odsekoma
odvedljiva (pogoj za integrabilnost funkcije in obstoj integralov za koeficiente).

f ( x)    (ak cos 2bkax  bk sin 2bkax )
a0
2
k 1
Koeficiente izračunamo po formulah:
b
a0  2  f ( x )dx
ba a
b
ak  b 2 a  f ( x) cos 2bnax dx
a
b
ak  b 2 a  f ( x)sin 2bnax dx
a
Primer 2.2: Iz grške mitologije je znan mit o Sizifu. Bogovi so ga kaznovali na
krut način. Iz kotanje je moral do vrha hriba kotaliti skalo, ki se mu je z vrha
hriba ponovno zakotalila v kotanjo. Popiši Sizifove muke z matematičnim
izrazom.
1,2
1
višina
0,8
0,6
0,4
0,2
0
0
0,1 0,2 0,3 0,4 0,5 0,6 0,7 0,8 0,9
čas
1
1,1
Slika 2.1. Grafični prikaz »Sizifovih muk«
Funkcija, ki popisuje enkratno Sizifovo muko, se zapiše:
 10t , 0  t  0.9 
h  f (t )   9 t

1  10 , 0.9  t  1 
Naš začetni interval je enak [0,1] in ga uporabimo v enačbi za razvoj v vrsto za
zapis neskončnih Sizifovih muk:
f ( x) 

a0
2
  (ak cos 2 k1 x  bk sin 2 k1 x )
k 1
J. Šrekl
Numerične metode v varnosti
24
Koeficiente izračunamo:
0.9
1
a0  2(  10t dt   (1 t ) dt )  1.08
10
0 9
0.9
0.9
1
0.9
1
ak  2(  10t cos(2k t ) dt   (1 t ) cos(2k t ) dt )
10
0 9
0.9
višina
bk  2(  10t sin(2 k t ) dt   (1 t ) sin(2 k t ) dt )
10
0 9
0.9
1,2
1
0,8
0,6
0,4
0,2
0
0
2
4
6
8
čas
Slika 2.2. Grafični prikaz trajne Sizifove muke
10
12
14
J. Šrekl
Numerične metode v varnosti
25
Večerno nebo – izziv za širjenje toplote
Transportna enačba
3
3.1. Teoretične osnove
V diferencialnih enačbah kjer imamo dve ali več neodvisnih spremenljivk, namesto
običajnih odvodov srečujemo parcialne odvode (odvode odvisne spremenljivke po
posameznih neodvisnih spremenljivkah). V koordinatnem sistemu to pomeni hitrost
spreminjanja funkcije v smereh koordinatnih osi. Kvazilinearno parcialno
diferencialno enačbo drugega reda (PDE) imenujemo enačbo:
 + 2 +  = (, , ,  ,  )
Glede na diskriminanto koeficientov  −  2 ločimo tri tipe enačb:
eliptični tip (EDE), če je  −  2 > 0 (primer: Laplaceova enačba),
parabolični tip (PDE), če je  −  2 = 0 (primer: transportna enačba)
hiperbolični tip (HDE), če je  −  2 < 0 (primer: enačba nihanja ali enačba
strune)
V enačbi pomeni:
2 
2 
2 
 =  2 ,  = ,  =  2


 =  ,  = 
V praksi, posebej še na področju požarne varnosti, se pogosto srečujemo s problemi
prevajanja toplote. Cela vrsta problemov je povezana z zagotavljanjem ustrezne
temperature ali pa ugotavljanjem vzrokov spremenjene temperature.
Razmišljanje o prenosu toplote bomo začeli na enostavnem modelu. Vzemimo palico
z dovolj majhnim presekom ω, prevodnostjo K , gostoto ρ, specifično toploto c ,
razpršenostjo κ in površinsko prevodnostjo H. Oglejmo si pretok toplote v preseku
palice pri x in pri x+dx. Z u označimo temperaturo, pretok v dani točki je odvisen
od prevodnosti, preseka in hitrosti spreminjanja temperature. Zapišimo pretok toplote
v točki x in x + dx. Pretok v točki, ki se za diferencial neodvisne spremenljivke
pomakne naprej, je povečan za diferencial pretoka.
J. Šrekl
Numerične metode v varnosti
x
−
26
x+dx



(−

2
−  2 )


Slika 3.1. Model za pretok toplote
Sprememba pretoka na odseku z dolžino dx je:
 2
 2 

Zaradi sevanja na površini se izgubi
ℎ( − 0 )
toplote, kjer je u0 temperatura okolice in p obseg toplotnega vodnika. Celotna
sprememba toplote pa je enaka:




Združimo vse ugotovitve, sprememba toplote je enaka spremembi zaradi toplotnega
toka in izgubi toplote v okolico:
ut  cK uxx  cHp
 (u  u0 )
Hp
Vstavimo še: c 

K
in c
 :
ut   uxx  (u  u0 )
Upoštevali bomo, da je prevodnost vodnika v prostor zelo majhna v primerjavi s
pretokom toplote v vodniku. Ta model lahko uporabimo tudi za steno, kjer je pri dani
globini temperatura enaka (zanemarili smo rob stena). Zato lahko krajše zapišemo
enačbo, ki je znana kot transportna enačba:
ut   uxx
3.2. Analitična rešitev
Splošna ideja za reševanje te parcialne diferencialne enačbe je razcep na neodvisni
funkciji:
(, ) = (). ()
Vstavimo v transportno enačbo in dobimo:
J. Šrekl
Numerične metode v varnosti
27
 ′′ 1  ′
=
= −


Dobimo dve ločeni navadni diferencialni enačbi s konstantnimi koeficienti
X    X  0
T   T  0
Določimo karakteristično enačbe za prvo navadno diferencialno enačbo drugega reda.
k2    0
Korena sta: k1,2  i

Splošna rešitev prve diferencialne enačbe se zapiše:
X ( x)  A cos  x  B sin  x
Rešitev druge enačbe je:
T (t )  e  t
Če želimo sestaviti celotno rešitev, moramo poznati robna pogoja in začetni pogoj. Za
začetek vzemimo palico dolžine L, ki ima izoliran plašč, tako da se toplota širi samo
vzdolž osi x.
Slika 3.2. Model prevodnika dolžine L [1]
Na krajiščih imamo enako temperaturo (0, ) = (, ) = 0, in temperatura je po
palici v začetku porazdeljena s funkcijo f(x)
V prvo enačbo vstavimo robne pogoje in dobimo:
A  0 in B sin  L  0
To pa pomeni:
   nL 
2
Celotna rešitev se zapiše:

u ( x, t )   Bne
n 1
 ( nL )2  t
sin( nL x)
Vstavimo začetni pogoj in dobimo:
L
Bn  L2  f ( x)sin nL x dx
0
J. Šrekl
Numerične metode v varnosti
28
Slika 3.3. Porazdelitev temperature v času in prostoru [1]
3.2. Numerično reševanje parabolične PDE
Ideja za numerično reševanje vseh treh tipov enačb je, da diferenciale zamenjamo z
diferencami. Splošna oblika poenostavljene transportne enačbe je:
 =  2 
PDE enačbo lahko še nekoliko poenostavimo, če zapišemo  = 1 in  = 1.
Spremenljivka x je definirana na intervalu (dolžini prevodnika) 0 ≤  ≤  in čas
teče od nič proti neskončno:  ≥ 0. Zapišimo enačbo z začetnimi (porazdelitev
temperature v prevodniku na začetku) in robnimi (gibanje temperature na konceh
vodnika).
 =  za 0 ≤  ≤  in  ≥ 0
(, 0) = ()
(0, ) = (1, ) = 0
Enostavno končno diferenčno aproksimacijo za diferencialno enačbo zapišemo:
1
1
(,+1 − , ) = 2 (+1, − 2, + −1, )

ℎ
Iz diferenčne enačbe izračunamo rekurzivno formulo za algoritem računanja
ui,j+1 = (1 − 2)ui,j + (ui+1,j + ui−1,j )
pri čemer velja, da je metoda konvergentna pri pogoju:
r=
k
1
≤
h2 2
J. Šrekl
Numerične metode v varnosti
29
Slika 3.5: Numerično računanje PDE [2]
Primer 3.1: Vzemimo kovinsko palico dolžine 1 enota z izoliranim plaščem. Na
krajiščih vzdržujemo temperaturo u = 0 (temperaturo okolja), temperatura znotraj
palice pa je v začetku porazdeljena s funkcijo () =   . Ogledali si bomo
spreminjanje temperature v času 0 ≤  ≤ 0,63. Da bo reševanje konvergentno, bomo
1
vzeli  = 4 = 0,25 . Za korak po dolžini palice bomo vzeli ℎ = 0,2, to pa pomeni, da
ℎ2
je  = 4 = 0,01.
Prvo vrstico izračunamo s pomočjo začetnega pogoja:
0
=0,2+B1
=0,2+C1
=0,2+D1
=0,2+E1
=0,2+F1
=SIN(B1*3,14159265)=SIN(C1*3,14159265)=SIN(D1*3,14159265)=SIN(E1*3,14159265) =SIN(F1*3,14159265) 0
0
Dobimo vrednosti
0
0
0
0,2
0,4
0,6
0,8
1
0,587785 0,951057 0,951057 0,587785 0
V nadaljnjih korakih si pomagamo z rekurzivnimi formulami:
0
=A2+0,25*0,04
0
=0,2+B1
=0,2+C1
=0,2+D1
=0,2+E1
=0,2+F1
=SIN(B1*3,14159265) =SIN(C1*3,14159265)=SIN(D1*3,14159265)=SIN(E1*3,14159265) =SIN(F1*3,14159265) 0
=SIN(B2*3,14159265) =0,5*C2+0,25*(B2+D2)=0,5*D2+0,25*(C2+E2)=0,5*E2+0,25*(D2+F2)=0,5*F2+0,25*(E2+G2)0
Tabela 3.1. Rezultati računanja za nekaj korakov:
0
0,01
0,02
0,03
0,04
0
0
0
0
0
0
0,2
0,587785
0,531657
0,480888
0,434967
0,393432
0,4
0,951057
0,860239
0,778093
0,703792
0,636586
0,6
0,951057
0,860239
0,778093
0,703792
0,636586
0,8
0,587785
0,531657
0,480888
0,434967
0,393432
1
0
0
0
0
0
J. Šrekl
Numerične metode v varnosti
0,05
0,06
0
0
30
0,355862 0,575797 0,575797 0,355862 0
0,32188 0,520813 0,520813 0,32188 0
Z barvo so označene celice, ki so povezane s formulo.
Celoten rezultat spreminjanja temperature s časom pa je viden na diagramu. (Časi za
posamezne krivulje so zapisani spodaj.)
Graf ohlajanja kovinske palice
1,2
1
0,8
0,6
0,4
0,2
0
0
0,2
0
0,4
0,01
0,02
0,6
0,03
0,8
0,04
Slika 3.6: Graf ohlajanja
0,05
1
0,13
1,2
0,63
J. Šrekl
Numerične metode v varnosti
31
Pričakovana ravna pot
4
4.1.
Regresija
Korelacije
V populaciji, ki jo opazujemo sta dve vrsti izidov dogodkov (npr. količino izdelkov in
dobiček ali pa količina proizvodov in velikost delovnih površin ). V prvem primeru
pričakujemo, da sta izida dogodkov odvisna, v drugem primeru pa to ni nujno. Torej
sta dve spremenljivki X in Y, ki merita izide lahko neodvisni ali odvisni
Spremenljivki združimo v urejen par (X,Y) (legi spremenljivk sta v paru določeni), ki
ga imenujemo slučajni vektor.
Povezanost spremenljivk v slučajnem vektorju imenujemo korelacija. Mera za
stopnjo korelacije je kovarianca. Kovarianca meri razliko med matematičnim
upanjem (srednjo vrednostjo) produkta spremenljivk in produktom matematičnih
upanj posameznih spremenljivk. Kovarianca je 0 pri neodvisnih spremenljivkah in je
različna od nič pri odvisnih spremenljivkah. Če označimo z 1 in  2 matematični
upanji posameznih spremenljivk, se kovarianca zapiše
K ( X ,Y )  E ( X .Y )  E ( X ).E (Y )  E ( XY )  12
Ker pa se spremenljivke merijo običajno z različnimi merami (in enotami), je veliko
bolj smiselna mera odvisnosti standardizirana korelacija ali korelacijski koeficient.
Dobimo ga tako, da korelacijo delimo s standardnima odklonoma (kvadratni koren iz
variance) obeh spremenljivk. Če sta  1 in
korelacijski koeficient izračuna iz enačbe
2
 22 varianci obeh spremenljivk, potem se
   ( X ,Y ) 
K ( X ,Y )
 1 2
Kadar govorimo o slučajnem vzorcu in slučajnih spremenljivkah v vzorcu, bo
kovarianco in korelacijski koeficient računali iz vzorčnih vrednosti. Vzorčno
kovarianco izračunamo iz vrednosti v vzorcu velikosti n

̅
̅ )
 = ∑(  − 
=
 =

−
Vzorčni korelacijski koeficient pa dobimo iz
J. Šrekl
Numerične metode v varnosti
32

̅) 
 = ∑( − 
=

̅ )
 = ∑( − 
=
 =
4.2.

√ 
Empirični model
Veliko statističnih problemov je povezano z iskanjem povezav med dvema ali več
spremenljivkami. Nekaj takih smo že srečali v prejšnjem poglavju o hipotezah.
Vendar pogosteje uporabljamo statistično tehniko, ki jo imenujemo regresijska
analiza.
Poglejmo si primer kemičnega destilacijskega procesa za pridobivanje kisika.
Spremenljivka Y čistot pridobljenega kisika, spremenljivka X pa je delež
ogljikovodika v glavnem kondenzatorju destilacijske enote.
Dobili smo naslednje podatke:
Tabela 4.1. Meritve odvisnosti čistosti kisika od deleža ogljikovodika v mešanici [3]
Podatke lahko narišemo tudi v diagram:
J. Šrekl
Numerične metode v varnosti
33
Slika 4.1. Diagram podatkov iz tabele 4.1 [3]
Če bi želeli poiskati funkcijsko zvezo med obema količinama, bi morali najti
analitično rešitev: funkcijo katere graf bi zajel vse izmerjene točke. V splošnem to ni
mogoče, kasneje bomo videli tudi razloge za to. Torej ni mogoče potegniti krivulje z
analitično enačbo, ki bi šla skozi vse točke diagrama. Zato se lotimo iskanja najbolj
primernih približkov. Dela se lotimo s statističnimi metodami. Ker vrednost čistosti
kisika ni odvisna samo od količine ogljikovodika v mešanici, ampak še od drugih
slučajnih vplivov procesa in napak pri meritvah bomo vzeli slučajno spremenljivko Y
(čistost kisika) pri posameznih vrednostih x (delež ogljikovodika). Matematično
upanje spremenljivke Y je pri vsaki vrednosti x:
E (Y | x)  Y |x  0  1x
To pomeni, da je slučajno spremenljivko Y mogoče napisati kot funkcijo x in
slučajnenapake .
Y  0  1 x  
Varianca te spremenljivke dobimo:
V (Y | x)  V (0  1x   )  V (0  1x)  V ( )  0   2
Na sliki si lahko naše zveze predstavimo z grafi:
Slika 4.2. Zveza med deležem ogljikovodika in čistostjo kisika [3]
Enačba
Y  0  1 x  
J. Šrekl
Numerične metode v varnosti
34
je model enopstavne linearne regresije, pri čemer je x regresor ali prediktor Y pa je
odvisna spremenljivka.
4.3.
Regresijska premica
Iščemo torej premico, ki se najbolje prilega vrednostim v diagramu (xy). Oceno
regresije dobimo z metodo najmanjših kvadratov (vsota kvadratov razdalj točk od
premice je minimalna)
Slika 4.3 Lega regresijske premice glede na izmerjene vrednosti [3]
Koeficiente regresijske premice ocenimo z metodo najmanjših kvadratov (to pomeni,
da je vsota kvadratov razlik med izmerjenimi in izračunanimi vrednostimi y
minimalna). Iščemo torej minimalno vrednost izraza

∑( − ̂ )2
=1
Z metodo najmanjših kvadratov dobimo koeficiente
̂0 = ̅ − ̂1 ̅
̂1 =
∑=1  ∙  − ̅ ∙ ̅ 
=
∑=1 2 −  ∙ ̅ 2

Ocena za regresijska premico se potem zapiše:
̂ = ̂0 + ̂1 
Izmerjeno vrednost dobimo kot vrednost na regresijski premici plus residual.
yi  ˆ0  ˆ1xi  ei
J. Šrekl
Numerične metode v varnosti
35
Slika 4.4. Regresijska premica skozi »oblaček« podatkov [3]
ei  yi  yˆi imenujemo ostanek (residual) in je razlika med prilagoditvijo modela in ite vrednosti spremenljivke. Ostanek se uporablja za oceno variance. Označimo vsoto
kvadratov napak:
n
n
SS E   e  ( yi  yˆi )2
i 1
2
i
i 1
Ker to ni nepristranska cenilka za varianco, bomo uporabili:
ˆ 2 
SS E
n2
Izračunamo še:
n
n
SS E   ( yi  yˆi )   ( yi  y  ˆ1 ( xi  x )) 2 
2
i 1
i 1
n
 ( y  y )
i 1
i
2
 2ˆ1 ( yi  y )( xi  x )  ˆ1 ( xi  x ) 2   SST  ˆ1S xy
n
n
i 1
i 1
SST   ( yi  y ) 2   yi2  ny 2
Primer 4.1: Prilagodimoenostavni linearni regresijski model podatkom iz naloge o
čistosti kisika v začetku poglavja
Tabela 4.2. Meritve odvisnosti čistosti kisika od deleža ogljikovodika v mešanici
Zap.št.
Delež ogljikovodikov v %
Čistost kisika v %
1
0,99
90,01
2
1,02
89,05
3
1,15
91,43
4
1,29
93,74
5
1,46
96,73
J. Šrekl
Numerične metode v varnosti
36
6
1,36
94,45
7
0,87
87,59
8
1,23
91,77
9
1,55
99,42
10
1,4
93,65
11
1,19
93,54
12
1,15
92,52
13
0,98
90,56
14
1,01
89,54
15
1,11
89,85
16
1,2
90,39
17
1,26
93,25
18
1,32
93,41
19
1,43
94,98
20
0,95
87,33
Za računanje uporabimo naslednje formule (uporabljen MATTAB):
Za računanje v EXCEL-u si lahko izdelamo tabelo (enako bi računali tudi s
kalkulatorjem):
Delež
ogljikovodikov
Čistost
Zap.št.
v%
kisika v % x*x
y*y
x*y
y^
e
e^2
1
0,99
90,01
0,9801
8101,8
89,1099
89,08132
0,928681 0,862448
4,62465
2
1,02
89,05
1,0404 7929,903
90,831
89,52974
-0,47974 0,230154
9,67521
3
1,15
91,43
1,3225 8359,445 105,1445
91,47292
-0,04292 0,001842
0,53363
4
1,29
93,74
1,6641 8787,188 120,9246
93,56556
0,174437 0,030428
2,49482
5
1,46
96,73
2,1316 9356,693 141,2258
96,10663
0,623365 0,388584
20,88033
6
1,36
94,45
1,8496 8920,803
128,452
94,61189
-0,16189 0,026207
5,24181
7
0,87
87,59
0,7569 7672,008
76,2033
87,28762
0,302378 0,091433
20,88947
8
1,23
91,77
1,5129 8421,733 112,8771
92,66871
-0,89871 0,807687
0,15249
9
1,55
99,42
2,4025 9884,336
97,45191
1,968092 3,873387
52,70034
154,101
J. Šrekl
Numerične metode v varnosti
37
10
1,4
93,65
131,11
95,20979
-1,55979 2,432932
2,21861
11
1,19
93,54
1,4161 8749,732 111,3126
92,07082
1,469185 2,158504
1,90302
12
1,15
92,52
1,3225
8559,95
106,398
91,47292
1,047084 1,096385
0,12924
13
0,98
90,56
0,9604 8201,114
88,7488
88,93184
1,628156 2,650891
2,5616
14
1,01
89,54
1,0201 8017,412
90,4354
89,38027
0,159731 0,025514
6,86702
15
1,11
89,85
1,2321 8073,023
99,7335
90,87502
-1,02502 1,050659
5,33841
16
1,2
90,39
1,44 8170,352
108,468
92,22029
-1,83029 3,349961
3,13467
17
1,26
93,25
1,5876 8695,563
117,495
93,11714
0,132861 0,017652
1,18701
18
1,32
93,41
1,7424 8725,428 123,3012
94,01399
-0,60399 0,364801
1,56125
19
1,43
94,98
2,0449
95,65821
-0,67821 0,459969
7,94958
20
0,95
87,33
0,9025 7626,529
88,48342
-1,15342 1,330378
23,33373
23,92
1,96 8770,323
9021,2 135,8214
82,9635
1843,21 29,2892 170044,5 2214,657
SSE=
21,24982 SST=173,3769
S(XY)
10,17744
S(XX)
0,68088
beta(1) 14,94748
beta(0) 74,283314
Lahko pa določimo parametre direktno (data-data analisys-regression)
Regression
Statistics
Multiple R
R Square
Adjusted R Square
Standard Error
Observations
Intercept
X Variable 1
0,936715381
0,877435705
0,870626578
1,086529053
20
Coefficients
Standard Error
74,28331424
1,593473376
14,94747973
1,31675827
Regresijsko premico lahko določimo tudi s pomočjo načrtovanja grafa v EXCEL
(Office 10) V raztreseni diagram dodamo regresijsko premico in funkcijsko zvezo za
linearno regresijo.
čistost kisika
Čistost kisika v %
105
100
95
90
85
y = 14,947x + 74,283
R² = 0,8774
Čistost kisika v %
0,8 0,9 1 1,1 1,2 1,3 1,4 1,5 1,6
Linearna (Čistost
kisika v %)
delež ogljikovodika
Slika 4.5. Rezultat pridobljen s pomočjo EXCEL
J. Šrekl
Numerične metode v varnosti
38
Primer 4.2. V reviji Journal of Sound and Vibration (Vol. 151, 1991, pp. 383-394) je v
članku opisana povezava med izpostavljenostjo hrupu in povišanim krvnim tlakom.
Priloženi so izmerjeni podatki (y – povišanje tlaka v mm, x – hrup v decibelih):
Povečan
krvni tlak
1
a)
b)
c)
Stopnja
hrupa
60
0
63
1
65
2
70
5
70
1
70
4
80
6
90
2
80
3
80
5
85
4
89
6
90
8
90
4
90
5
90
7
94
9
100
7
100
6
100
Nariši diagram. Ali je primeren enostavni linearni model regresije?
Prilagodi linearni model po metodi najmanjših kvadratov. Poišči oceno za
varianco.
Poišči pričakovano povprečno razliko krvnega tlaka pri zgornji meji
dovoljenega hrupa 85 decibelov.
Rešitev:
(a) S pomočjo grafa v EXCELU narišemo (raztreseni – scater graf):
J. Šrekl
Numerične metode v varnosti
39
Povečan krvni tlak
Odvisnost povišanja krvnega tlaka od hrupa
10
8
6
4
2
0
Odvisnost povišanja krvnega
tlaka od hrupa
50
60
70
80
90
100
Stopnja hrupa
(b) Iz slike vidimo, da je linearna regresija primerna, zato spremenimo v diagram
z regresijo:
10
y = 0,1828x - 10,764
R² = 0,7582
Povečan krvni tlak
8
Odvisnost povišanja
krvnega tlaka od hrupa
6
regresija
4
2
Linearna (Odvisnost
povišanja krvnega tlaka od
hrupa)
0
50
60
70
80
90
100
Stopnja hrupa
Ocena za varianco:
SST=
SSE=
sigma^2=
124,2
31,2664749
1,737026383
Poiščimo še pričakovano razliko krvnega tlaka pri mejni vrednosti 85 dB.
Tabela 4.3. Rezultat računanja z EXCEL-om
Povečan Stopnja
e
e^2
krvni tlak
hrupa y^
1
60 0,326098 0,673902 0,454143
0
63
1
65
2
70
5
70
1
70
4
80
0,84898
1,197568
2,069038
2,069038
2,069038
3,811977
-0,84898
-0,19757
-0,06904
2,930962
-1,06904
0,188023
0,720767
0,039033
0,004766
8,59054
1,142842
0,035353
J. Šrekl
Numerične metode v varnosti
40
6
90 5,554916 0,445084 0,198099
2
80 3,811977 -1,81198 3,283261
3
80 3,811977 -0,81198 0,659307
5
85 4,683447 0,316553 0,100206
4
89 5,380622 -1,38062 1,906118
6
90 5,554916 0,445084 0,198099
8
90 5,554916 2,445084 5,978434
4
90 5,554916 -1,55492 2,417765
5
90 5,554916 -0,55492 0,307932
7
94 6,252092 0,747908 0,559366
9
100 7,297856 1,702144 2,897295
7
100 7,297856 -0,29786 0,088718
6
100 7,297856 -1,29786 1,684429
Vse to lahko najdemo z enim korakom ob uporabi analize podatkov (Data Analisys) –
Regresija:
Dobimo osnovno statistiko podatkov:
Regression Statistics
Multiple R
R Square
Adjusted R Square
Standard Error
Observations
0,865018523
0,748257046
0,734271326
1,317962967
20
Koeficienti regresije in napake:
Intercept
Stopnja hrupa
Coefficients
-10,13153766
0,174293933
Standard Error
1,994899887
0,023828641
Predicted
Povečan krvni
tlak
0,326098326
0,848980126
1,197567992
2,069037657
2,069037657
2,069037657
3,811976987
5,554916318
3,811976987
3,811976987
4,683446653
5,380622385
Residuals
0,673901674
-0,848980126
-0,197567992
-0,069037657
2,930962343
-1,069037657
0,188023013
0,445083682
-1,811976987
-0,811976987
0,316553347
-1,380622385
t Stat
P-value
-5,07871985 7,83221E-05
7,314472274 8,56676E-07
Napake pri oceni:
Observation
1
2
3
4
5
6
7
8
9
10
11
12
Standard
Residuals
0,525332023
-0,661812345
-0,154011775
-0,053817483
2,284796786
-0,833355573
0,1465711
0,346959683
-1,412505079
-0,63296699
0,246765392
-1,076247737
J. Šrekl
Numerične metode v varnosti
13
14
15
16
17
18
19
20
5,554916318
5,554916318
5,554916318
5,554916318
6,25209205
7,297855649
7,297855649
7,297855649
0,445083682
2,445083682
-1,554916318
-0,554916318
0,74790795
1,702144351
-0,297855649
-1,297855649
41
0,346959683
1,906035863
-1,212116496
-0,432578406
0,583022734
1,326886356
-0,232189823
-1,011727913
Residuals
Stopnja hrupa Residual Plot
4
3
2
1
0
-1 0
-2
-3
20
40
60
80
100
120
Stopnja hrupa
Slika 4.6. Graf residualov –ostankov
10
5
Povečan krvni tlak
Predicted Povečan krvni tlak
0
0
20
40
60
80
100
120
Stopnja hrupa
Slika 4.7. Graf odvisnosti in regresijska prilagoditev.
Povečan krvni tlak
Povečan krvni tlak
Stopnja hrupa Line Fit Plot
Normal Probability Plot
10
5
0
0
20
40
60
80
Sample Percentile
100
Slika 4.8. Verjetnost povečanja krvnega tlaka merjena s percentili.
120
J. Šrekl
4.4.
Numerične metode v varnosti
42
Nelinearna regresija
V praksi se srečujemo z meritvami, ki ne kažejo na linearno povezavo dveh
spremenljivk, ki jih merimo. Pogosto je ta povezava eksponentna. Recimo, da
pričakujemo eksponentno zvezo med dvema spremenljivkama.
 =  ∙  
Enačbo logaritmiramo:
ln  = ln  +  ∙ 
Če zamenjamo 0 = ln ,
regresijsko premico:
1 =  in  = ln , upoštevamo še napako, dobimo
Y  0  1 x  
Če želimo iz izmerjenih podatkov določiti krivuljo, bomo meritve za odvisno
spremenljivko logaritmirali:
Tabela 4.4. Meritve za eksponentno odvisnost
x
y
 = ln 
x1
y1
Y1
x2
y2
Y2
x3
y3
Y3
x4
y4
Y4
Ko ocenimo parametra:
̂0 = ̅ − ̂1 ̅
∑=1  ∙  − ̅ ∙ ̅ 
=
∑=1 2 −  ∙ ̅ 2

Iz ocenjenih parametrov izračunamo ocenjene koeficiente enačbe:
̂1 =
̂
̂ =  0
̂ = ̂1
Primer 4.3. Število bakterij v poskusni koloniji raste po naravnem zakonu rasti, torej z
eksponentno funkcijo:
 = 0  
Pri preizkusu smo dobili naslednje rezultate:
J. Šrekl
Numerične metode v varnosti
43
Tabela 4.4. Meritve števila bakterij v eksperimentu
Čas (dan)
število bakterij
1
1229
2
3026
3
7439
4
18300
5
45008
6
110704
7
272285
8
9
669716
1647234
10
4051542
11
9965185
12
24510401
13
60285857
Ali poskus ustreza naravnemu zakonu rasti?
Reševanje: Tabelo dopolnimo z logaritmi
Čas (dan)
število bakterij
1
1229
2
3026
3
7439
4
18300
5
45008
6
110704
7
272285
8
669716
9
1647234
10
4051542
11
9965185
12
24510401
13
60285857
logaritem števila
7,11395611
8,014996894
8,91449171
9,814656339
10,71459553
11,61461525
12,51460459
13,41460902
14,31460808
15,21460811
16,11460808
17,01460812
17,91460809
In dobimo rezultat:
20
15
10
y = 0,9x + 6,2145
5
0
0
2
4
6
8
10
12
14
J. Šrekl
Numerične metode v varnosti
44
Iz enačbe:
 = 0,9 + 6,2145
kjer upoštevamo še:
̂
̂ =  0 =  6,2145 = 499,946
dobimo:
 = 499,946 ∙  0,9∙
Če si pogledamo še graf:
število bakterij
70000000
y = 499,94e0,9x
60000000
50000000
40000000
30000000
20000000
10000000
0
0
5
10
15
Vidimo, da je rast bakterij zares eksponentna po naravnem zakonu. V EXCEL-u lahko
računamo direktno, če za izmerjene podatke vzamemo raztreseni graf in dodamo
eksponentno regresijo.
Prvi korak: Za podatke izberemo raztreseni grafikon:
J. Šrekl
Numerične metode v varnosti
45
Drugi korak: V meniju orodij za grafikone – Načrt kliknemo: dodaj elemente
grafikona
V meniju izberemo Oblikuj trendno črto – možnosti trendne črte in označimo
»Eksponentna« in »Prikaži enačbo na grafikonu«.
Opomba: Uporabljali smo OFFICE 2013 vendar je podobno v OFFICE 2010 ali
starejšem.
J. Šrekl
Numerične metode v varnosti
4.5.
Test hipoteze o enostavni linearni regresiji
46
Pomembno področje testiranja hipotez je testiranje hipotez o linearni regresiji.
Predpostavimo vrednost parametra 1 in poskušamo potrditi ali zavrniti hipotezo o
vrednosti parametra.
H 0 : 1  1,0
H1 : 1  1,0
Zapisali smo dvostransko alternativno hipotezo. Ker so napake porazdeljene
neodvisno po normalnem porazdelitvenem zakonu z matematičnim upanjem nič in
varianco
zakonih
2,
so slučajne spremenljivke
Yi
porazdeljene po neodvisnih normalnih
N (0  1xi , 2 )
To pa pomeni, da je 1 porazdeljen po normalnem zakonu:
N ( 1 , Sxx )
2
Iz tega sledi, da je statistika:
T0 
ˆ1  1,0
ˆ 2 / S xx
porazdeljena po Studentovem zakonu S(n-2).
Podobno postavimo hipotezo za parameter 0
H 0 :  0   0,0
H1 :  0   0,0
Uporabimo statistiko:
T0 
ˆ0  0,0
ˆ 2  1n  Sx 
2
xx
J. Šrekl
Numerične metode v varnosti
47
Pomembna je posebna hipoteza, ki ugotavlja signifikanco ali značilnost regresije:
H 0 : 1  0
H1 : 1  0
,
ki ob potrditvi izključuje zvezo med x in Y
Slika 4.8. Pprimera, ko ničelne hipoteze ne morem zavrniti, torej ni potrjena povezava
med x in Y.
Primer 4.4. Testiramo značilnost (signifikantnost) regresije za nalogo o destilaciji
kisika.
1. Parameter, ki ga obravnavamo, je korelacijski koeficient 1..
2. Ničelna hipoteza: H 0 : 1  0 Če je korelacijski koeficient enak nič, sta
spremenljivki neodvisni. Ni povezave med dvema spremenljivkama.
3. Alternativna hipoteza: H 0 : 1  0
4. Stopnja tveganja:
  0.01
5. Testna statistika:
t0 
6. Kriterij zavračanja:
ˆ1
2
S xx
t0.005,18  2.88
7. Izračun:
8. Sklep: Hipotezo zavrnemo, kar pomen, da obstaja povezava med obema
spremenljivkama.
Poglejmo si regresijsko analizo izračunano z MINITAB:
J. Šrekl
Numerične metode v varnosti
48
Iz analize vidimo, da je P-vrednost enaka 0,kar pomeni, da pri nobeni stopnji
zanesljivosti ne moremo testa sprejeti. Vedno ga zavračamo!
4.6.
Analiza variance za regresijo
Pri vsaki vrednosti prediktorja (neodvisne spremenljivke  ) obstajajo pri dani
meritvi tri vrednosti odvisne spremenljivke:
izmerjena vrednost  ,
̅ in
srednja vrednost 
̂.
regresijska ocena 
S pomočjo analize variance ANOVE testiramo ustreznost linearnega regresijskega
modela. Vzorčna varianca spremenljivke je sorazmerna z vsoto kvadratov razlik med
posameznimi vrednostmi in povprečno vrednostjo
n
 ( y  y)
i 1
2
i
Izraz v oklepaju lahko razstavimo na dva izraza:
yi  y  yi  yˆi  yˆi  y
Če izraz kvadriramo dobimo:
J. Šrekl
Numerične metode v varnosti
49
n
n
n

2
2
(
y

y
)

(
y

y
)

 i
 i
 ( yi  yˆi ) 2
i 1
i 1
i 1
SST  SS R  SS E
V Formulah pomenijo oznake:
yi - vzorčna vrednost spremenljivke,
y - povprečje vzornih vrednosti,
yi - pripadajoča vrednost na regresijski premici,
SST - vsota kvadratov korekcij,
SS R - regresijska vsota kvadratov,
SS E - vsota kvadratov napak
Celotna vsota kvadratov korekcij je torej razdeljena na regresijska vsota kvadratov ali
in vsota kvadratov napak
Ker je:
SS E  SST  ˆ1S xy
in
SST  ˆ1S xy  SS E
je zaradi gornje enačbe
SS R  ˆ1 S xy
Če pojasnjeni in nepojasnjeni del variance delimo s njunima prostostnima stopnjama,
dobimo povprečja kvadratov:
SS R ˆ1 S xy
MS R 

1
1
MS E 
SS E
n2
MSR - povprečje kvadratov regresije - pojasnjeni del variance,
MSE - povprečje kvadratov napak - nepojasnjeni del variance.
Sestavimo statistiko:
F0 
SS R / 1
MS R

SS E /( n  2) MS E
J. Šrekl
Numerične metode v varnosti
50
ki je porazdeljena po F1,n-2 porazdelitvi (Fisherjeva porazdelitev) in ima (1,(n-2))
prostostnih stopenj.
Tabela 4.3. Izračun F-testa
ANOVA
Regression
Residual
Total
Significance
F
df
SS
MS
F
1
SS R  1S xy
MS R
MS R / MS E
SSE  SST  1S xy
MS E
n-2
n-1
SST
V tabeli pomeni:
df – prostostne stopnje (degrees of freedom),
SS – vsota kvadratov (sum of squares),
MS – povprečje kvadratov (mean square)
Seveda pa je analiza variance povezana z t-testom, ki smo ga obravnavali v prejšnjem
razdelku.
T0 
ˆ1
T02 
 / S XX
2
ˆ12 S XX
MS E

ˆ1 S XY
MS E

MS R
MS E
V F-testu hipotezo zavračamo, če je
t02  f 0  f ,1,( n2) .
Primer 4.5. Testiramo značilnost (signifikantnost) regresije za nalogo o destilaciji
kisika s pomočjo ANOVE.
Rešitev: Računali bomo s pomočjo EXCEL-a. Seveda moramo imeti vključeno Analizo
podatkov (ki jo je potrebno ustrezno aktivirati)
Vnesemo podatke v EXCEL-ov list:
J. Šrekl
Numerične metode v varnosti
V glavnem meniju odpremo meni »PODATKI«:
Odpremo :
in dobimo podokno Data Analysis in potrdimo Regression
51
J. Šrekl
Numerične metode v varnosti
52
Odpre se nam okno Regression, v katero vpišemo niz celic za Y in niz celic za X.
Zapišemo tudi ime lista, kamor se izpišejo rezultati. Označimo še residuale
Rezultati se nam izpišejo v list, ki smo ga mi označili z ANOVA_Re
J. Šrekl
Numerične metode v varnosti
53
Nas seveda zanima samo analiza variance, ki je v tabeli:
Vidimo, da je vrednost Fisherjeve statistike 128,86 in kritična vrednost 1,23 ∙ 10−9.
Hipotezo zavračamo, ker vrednost Fisherjeve statistike presega kritično vrednost.
Če vrednost 0 = 11,35 iz primera 4.3. kvadriramo, dobimo:
02 = 11,352 = 128,82
kar je približno 0 –vrednost Fisherjeve statistike. (Razlika je nastala zaradi napak
pri računanju in zaokroževanju)
J. Šrekl
4.7.
Numerične metode v varnosti
54
Multipla regresija
Mnoge aplikacije regresijske analize nas pripeljejo do primerov, ko je odvisna
spremenljivka odvisna od več neodvisnih spremenljivk - regresorjev.
Na primer življenjska doba stružnega noža je odvisna hitrosti struženja in naklonskega
kota noža ob struženju.
Model za linearno multiplo regresijo zapišemo z enačbo
 = 0 + 1 1 + 2 2 + 
Oznake v enačbi pomenijo:
 – življenjska doba stružnega noža
1 – hitrost struženja (kotna hitrost pri vrtenju)
2 – naklonski kot noža (kot med smerjo obodne hitrosti in smerjo noža)
 - slučajna napaka merjenja
Rešitev predstavlja ploskev v tridimenzionalnem prostoru.
Model lahko posplošimo v n-dimenzionalno odvisnost:
 = 0 + 1 1 + 2 2 + ⋯ +   + 
Tako kot smo v razdelku 4.4. obravnavali nelinearno regresijo s pomočjo analize
linearno, lahko tudi kakšno drugo (polinomsko npr. kubno) obravnavamo s pomočjo
linearnega modela.
Kubna regresija ima obliko:
 = 0 + 1  + 2  2 + 3  3 + 
Z uporabo zvez:
1 = 
2 =  2
3 =  3
dobimo linearno multiplo regresijo:
 = 0 + 1 1 + 2 2 + 3 3 + 
V modelih, kjer imamo tudi medsebojno učinkovanje neodvisnih spremenljivk, bomo
model zapisali
 = 0 + 1 1 + 2 2 + 12 1 2 + 
Primer 4.6. Poglejmo si primer
 = 50 + 101 + 72 + 51 2
Tridimenzionalno rešitev lahko vidimo v projekciji grafa
J. Šrekl
Numerične metode v varnosti
Rešitev si lahko pogledamo tudi v konturnem grafu
Primer 4.7. Poglejmo si še nelinearni primer
 = 800 + 101 + 72 − 8,512 − 522 + 41 2
Tridimenzionalno rešitev lahko spet vidimo v projekciji grafa
ali pa v konturnem grafu
55
J. Šrekl
Numerične metode v varnosti
Parametre v enačbi računamo z metodo najmanjših kvadratov kjer s podatki vzorca
velikosti n rešimo linearni sistem enačb
V EXCEL-u ne obstajajo izdelani algoritmi za računanje multiple regresije, zato si
bomo več o tem pogledali v zadnjem poglavju
56
J. Šrekl
Numerične metode v varnosti
57
Primerjava različnih vzorcev nas vodi do pravila v populaciji
5
5.1.
Analiza variance za več obravnav
Analiza variance
Analíza variánce (kratica ANOVA) je v statistiki zbirka statističnih modelov in z njimi
povezanih postopkov, s katerimi skupno varianco med množicami podatkov
razcepimo na posamezne komponente, navadno zato, da preverimo, ali so razlike med
vzorci razložljive kot statistična odstopanja znotraj iste populacije. variance.
Izboljšanje uspešnosti proizvodnega procesa je odvisno od preizkusov posameznih
metod v tem procesu. S primerjavami rezultatov lahko presojamo o vplivu sprememb
parametrov v procesu. Večina procesov je odvisno od več obvladljiv spremenljivk,
kot so temperatura, tlak, in pomik. Eksperimentalne tehnike podprte s statističnimi
metodami so še posebej uporabne v inženirstvu. Imajo tudi obsežno uporabno vlogo
pri razvoju novih postopkov. S poskusi, lahko inženirji določijo, katere podmnožica
procesnih spremenljivk ima največji vpliv na uspešnost procesa. Rezultati takega
poskusa lahko omogoči:
1. Izboljšana donosnost postopka
2. Zmanjša variabilnost v procesu in tesnejše usklajenosti z nominalno ali ciljano
zahtevo
3. Zmanjša čas oblikovanja in razvoja.
4. Zmanjšati stroške poslovanja
Eksperimentalne metode načrtovanja so koristni tudi pri tehničnem projektiranju, kjer
se razvijajo novi izdelki oziroma izboljšujejo obstoječi. Nekatere tipične aplikacije za
statistično obravnavo poskusov pri konstruiranja vključujejo:
1. Vrednotenja in primerjave osnovnih konstrukcijskih konfiguracij
2. Vrednotenje različnih materialov
3. Izbira konstrukcijskih parametrov, tako da izdelek deluje tudi v različnih v
pogojih (ali da bo design dovolj robusten)
4. Določitev ključnih parametrov oblikovalskih zahtev, ki vplivajo na
učinkovitost izdelka
J. Šrekl
Numerične metode v varnosti
58
Uporaba eksperimentalnega modela v procesu konstruiranja lahko privede do
izdelkov, ki so lažji za izdelavo, ki se izdelujejo na varnejši način, ki so boljši na
trgu in zanesljivejši kot njihovi konkurenti, in proizvodi, ki jih je mogoče oblikovati,
razviti in proizvesti v krajšem času.
Poskusi so običajno izvedeni v nekem smiselnem zaporedju. To pomeni, da je prvi
preskus s kompleksnim sistemom (morda proizvodni proces), ki ima veliko
obvladljivih spremenljivk, zasnovan kot presejalni poskus, da določi, katere
spremenljivke so najbolj pomembne. Kasnejši poskusi se uporabijo za izboljšanje te
informacije in določijo, katere prilagoditve teh kritičnih spremenljivk so potrebne za
izboljšanje postopka. Končno, cilj poskusa je optimizacija, to pomeni izbor vrednosti
kritičnih spremenljivk, ki da najboljše rezultate procesa.
Vsak poskus vključuje niz aktivnosti:
1. Domneva - izvirna hipoteza, ki motivira poskus.
2. Poskus - preskus za preiskavo o ugibanjih.
3. Analiza - statistična analiza podatkov iz poskusa.
4. Zaključek - pregled spoznanj kot posledic poskusa o prvotnih ugibanjih.
Faktor je spremenljivka, ki vpliva na proučevani pojav. Predpostavimo da
preučujemo en sam faktor na različnih nivojih in želimo primerjati rezultate. Vsak
nivo zase imenujemo obravnava (treatment). Beseda izhaja iz statistične obravnave v
poljedelstvu, ko se je preučevalo enake posevke z različno obdelavo (gnojenje in
škropljenje). Metode analize variance je v 20-tih in 30-tih letih 20. stoletja prvi
uporabil statistik in genetik Ronald Fisher, zato so ponekod znani kot »Fisherjeva
ANOVA« ali »Fisherjeva analiza«.
Primer 5.1: Proizvajalca papirja, ki se uporablja za izdelavo nakupovalnih vrečk,
zanima izboljšanje natezne trdnosti izdelka. Inženir tehnolog meni, da je trdnost
funkcija deleža trdega lesa v proizvodnji celuloze in da je območje deleža trdega lesa
iz praktičnih razlogov najprimernejša med 5 in 20%.
Meritev trdnosti pri različnih deležev trdega lesa je pokazala:
Tabela 5.1. Meritve štirih tretmajev pri 5%, 10%, 15% in 20% (enako veliki
vzorci)[3]
Rezultat meritev lahko nazorno prikažemo z grafikonom kvantilov (box-and-whisker
plot)
J. Šrekl
Numerične metode v varnosti
59
Slika 5.1. Grafikon kvantilov vseh štirih obravnav
Če povzamemo iz primera vidimo, da imamo a nivojev enega merjenja faktorja
(spremenljivke, ki jo želimo meriti) in v vsakem nivoju n merjenj vzorca (kadar
govorimo o enako velikih vzorcih). Po izvedbi meritev a obravnav in n opazovanj v
posamezni obravnavi dobimo tabelo podatkov.
Tabela 5.1. Tabela meritev
Vzorčne vrednosti
Obravnava
Vsota
Povprečje
1
y11
y12
y1n
y1.
y1. ,
2
y21
y22
y2n
y2.
y2. ,
a
ya1
ya 2
yan
ya .
ya . ,
y..
y.. ,
a – število obravnav (tretmajev)
n – število opazovanj v obravnavi (velikost vzorca ene obravnave)
Izmerjene vrednosti spremenljivk v poskusu lahko zapišemo z linearnim statističnim
modelom:
 =  +  + 
V enačbi pomeni:
 - celotna sredina, vrednost, ki je neodvisna od obravnave (tretmaja)
 - učinek obravnave
J. Šrekl
Numerične metode v varnosti
60
 – napaka meritve (normalno porazdeljena vrednost s srednjo vrednostjo 0)
Povprečno vrednost v posamezni obravnavi zapišemo kot vsoto celotne sredine in
učinka obravnave:
 =  +  
V tem primeru bi se enačba modela zapisala:
 =  +  .
Izračunamo vsoto kvadratov, kjer je vir variabilnosti med vzorci
a
SS B  n ( yi.  y.. ) 2
i 1
in izračunamo še vsoto kvadratov kjer je vir variabilnosti znotraj vzorcev
a
n
SSW   ( yij  yi. ) 2 ,
i 1 j 1
ki pomeni seštevek kvadratov vseh razlik med izmerjenimi vrednostmi v vzorcih in
povprečno vrednostjo v vzorcu.
Srednje vrednosti ali variance dobimo tako, da vsote kvadratov delimo s prostostnimi
stopnjami.
MSW 
SSW
SS B
, MS B 
a(n  1)
(a  1)
Statistiko za ocenjevanje hipoteze nam da
F0 
MS B
,
MSW
ki je porazdeljena po Fisherjevem zakonu z (a-1) in a(n-1) prostostnimi stopnjami.
Izberemo število obravnav in v začetni fazi običajno preverjamo hipotezo o
povprečjih obravnav. V naslednjem koraku ocenjujemo učinke obravnav (imenujemo
ga fiksni efekt modela). Razvili bomo analizo variance za model fiksnega efekta.
Učinke obravnav označimo s
povprečja  tako, da je
i
in so običajno definirani kot deviacije celotnega
a

i 1
i
0
To pomeni enako kot testirati hipotezo
H 0 : 1   2 
a  0
H1 :  i  0 vsaj za en i
Primer 5.2: Petnajst dijakov je bilo razdeljeno naključno v tri skupine, kjer so jih po
različnih metodah poučevali matematiko. Na koncu semestra je vseh petnajst dijakov
pisalo enak test. Rezultati testa so bili naslednji:
J. Šrekl
Numerične metode v varnosti
Metoda I
48
73
51
65
87
Metoda II
55
85
70
69
90
61
Metoda III
84
68
95
74
67
Ali lahko z 1% stopnjo zanesljivosti zavrnemo ničelno hipotezo, da so aritmetična
povprečja rezultatov po vseh treh metodah poučevanja enaka.
Rešitev: Podatke skupaj z naslovi vpišemo v novi delovni list v Excelu. S klikom
odpremo Tools (Orodja) Data Analysis in izberemo Anova: Single Factor
Rezultati se izpišejo v tabeli: Anova: Single Factor
SUMMARY
Groups
Metoda I
Metoda II
Metoda III
Count
5
5
5
Sum
Average Variance
324
64,8
258,2
369
73,8
194,7
388
77,6
140,3
ANOVA
Source of
Variation
Between
Groups
Within Groups
SS
df
432,1333
2372,8
2
12
Total
2804,933
14
MS
F
P-value
F crit
216,0667 1,092717 0,366464 6,926608
197,7333
J. Šrekl
Numerične metode v varnosti
62
Prva tabela nam enostavno predstavi sumarne statistične podatke za vsako skupino
posebej (število opazovanj, vsota vseh točk, povprečje, varianca)
Druga tabela je tabela po analizi ANOVA.
Prvi stolpec je med vzorčna vsota kvadratov SS B in zunaj vzorčna vsota kvadratov
SSW .
Drugi stolpec so prostostne stopnje v prvi vrstici število skupin minus ena, v drugi
vrstici produkt števila skupin krat velikost skupine minus ena.
Tretji stolpec je MS B – varianca med skupinami in
SSW – varianca brez skupin.
MS B
).
MSW
Peti stolpec je P-vrednost, to je ploščina pod F-porazdelitvijo desno od vrednosti. To
vrednost primerjamo z izbranim α. Ker je vrednost mnogo večja od   0.01 ne
moremo zavrniti ničelne hipoteze, ki trdi, da so povprečja vseh treh populacij enaka.
Z drugimi besedami različne učne metode dajo različne rezultate.
Zadnji stolpec je kritična vrednost F-porazdelitve pri   0.01 , pomeni vrednost testne
statistike pri kateri bi zavrnili hipotezo.
Četrti stolpec nam da vrednost testne statistike F ( F0 
5.2.
Razlike med karakterističnimi skupinami
Če nas zanima, ali smo v posamezni obravnavi dosegli pričakovano vrednost, bomo s
t-testom poskušali potrditi ali zavrniti pričakovanja. Pogosteje pa bomo zgolj z
iskanjem intervala zaupanja poskušali ugotoviti, ali je povprečni rezultat obravnave
znotraj naših pričakovanj.
Izračunali bomo testno statistiko, ki meri razliko
F – test v ANOVA analizira razliko med karakteristikami v skupinah. Rezultat nam
pove ali je značilna (signifikantna) razlika med parametri, ki nas zanimajo. V prejšnji
nalogi smo gledali razlike parametrov na splošno, ne glede na to kateri dve skupini
primerjamo. Sedaj pa bomo pogledali primerjavo med posameznimi parametri v
posameznih skupinah. Ugotovili bomo kako konstruiramo interval zaupanja razlik za
posamezni par spremenljivk.
Interval zaupanja za m primerjanih razlik srednjih vrednosti i   j dobimo:
Yi  Y j  t / 2 m S
1 1
 ,
ni n j
kjer je
S  MS E ,
n
m = število mogočih intervalov zaupanja m    ,
k 
J. Šrekl
Numerične metode v varnosti
t / 2m = t-statistika z
df  n  k ,
63
n  n1   nk ,
k = število obravnav (tretmajev),
ni  je velikost i-te obravnave
Primer 5.3: Poskus se izvaja za določitev primanjkljaja vlage tal, ki izhajajo iz
različnih količin ostankov lesa, ki ostane po rezanja dreves v gozdu. Trije postopki so
naslednji:
Obravnava 1: ni ostankov lesa
Obravnava 2: 4,72 m 3 ostanka
Obravnava 3: 16,88 m 3 ostanka
Meritve primanjkljaja vlage, so podane v spodnji tabeli. Opravite test ANOVA in in
konstruiraj t-intervale zaupanja za razlike obravnav s 95% stopnjo zanesljivosti.
Obravnava 1
1,52
1,38
1,29
1,48
1,63
Obravnava 2
1,63
1,82
1,35
1,03
2,30
1,45
Obravnava 3
2,56
3,32
2,76
2,63
2,12
2,78
Rešitev:
V Excelu poženemo Orodja – Data Analysis – ANOVA – Single Factor.
J. Šrekl
Numerične metode v varnosti
64
Dobimo rezultat:
Anova: Single Factor
SUMMARY
Groups
Obravnava 1
Obravnava 2
Obravnava 3
Count
5
6
6
ANOVA
Source of
Variation
SS
Sum
Average Variance
7,3
1,46
0,01705
9,58 1,596667 0,189827
16,17
2,695
0,15103
df
MS
Between Groups 5,279128
Within Groups 1,772483
2 2,639564
14 0,126606
Total
16
7,051612
F
P-value
F crit
6,34E20,84866
05 3,738892
Ker je P-vrednost izjemno majhna pomeni, da se skupine obravnav razlikujejo.
Dobili smo naslednje rezultate (označeni v tabeli):
Y1  1.46 , Y2  1.597 , Y3  2.695
Velikost vzorcev posameznih obravnav so n1  5 ,
0.05
t-statistiko za 
in df = 14.
2m

23
 0.00833
n2  n3  6 . Izračunamo
TINV(0.00833,14)=2.718
Izračunamo še
S  MS E  0.126606  0.356
Uporabimo računanje v Excelu in določimo intervale zaupanja za posamezne razlike
srednjih vrednosti.
Average
1,46
1,596667
2,695
S
0,356
Razlike Krit.vred. min
2;1
0,136667 0,585916 -0,44925
3;2
1,098333 0,558649 0,539685
3;1
1,235 0,585916 0,649084
t
2,718
max
0,722582
1,656982
1,820916
Lahko zapišemo še končni rezultat intervale zaupanja za posamezne razlike srednjih
vrednosti po dveh obravnav.
 2  1 : (-0.45, 0.72)
 3   2 : (0.54, 1.66)
 3  1 : (0.65, 1.82)
J. Šrekl
Numerične metode v varnosti
65
Primer 5.4.: Poskus so izvedli za določanje učinka štirih različnih kemikalij na
trdnost tkanine. Te kemikalije se uporabljajo kot del stalnega potiska tkanin v
zaključnem postopku. Pet tkanina v vzorci je bilo izbranih naključno, celoten blok so
dobili s preskušanjem vsake kemikalije v naključnem vrstnem redu na vsakem vzorcu
tkanine. Omenjeni podatki so prikazani v tabeli. Mi bo za test razlike uporabili
ANOVA z α = 0,01.
Tabela 5.2. Rezultati vzorčenja tkanin
V MATLAB-u so dobili rezultate:
Izračunali smo tudi kritično vrednost F-testa: 0,01;3;12 = 5,95. Glede na rezultat v
tabeli (izračunana testna vrednost 0 = 75,13) lahko sklepamo, da hipotezo o
neznačilnem vplivu obravnave moramo zavrniti. Povprečja posameznih vzorcev so
značilno različna, kemični tipi torej vplivajo na trdnost tkanin.
Značilne ali neznačilne razlike med povprečnimi vrednostmi posamezni obravnav pa
testiramo s Fisherjeva LSD metoda.
Izračunamo povprečne vrednosti obravnav:
Izračunamo še LSD vrednost
J. Šrekl
Numerične metode v varnosti
66
Uporabili smo α = 0,05 in izračunali t-kritična vrednost:
ki smo jo vstavili v gornjo enačbo.
Preverjamo razlike med povprečji:
Razlike si lahko pogledamo tudi na grafu:
Slika 5.1. Rezultat Fisherjeva LSD testa
Primer 5.5: Članek v Environment International (Vol 18, No 4, 1992) opisuje
eksperiment, v katerem se raziskuje količina radona ki se sprošča v prhah tušev. V
poskusu je bila uporabljena z radonom obogatena voda, pri in šestih različnih
premerih odprtine tušev. Podatki iz poskusa so prikazani v spodnji tabeli.
Diameter
Observation 1
Observation 2
Observation 3
Observation 4
0,37
80
83
83
85
0,51
75
75
79
79
0,71
74
73
76
77
1,02
67
72
74
74
1,4
62
62
67
69
1,99
60
61
64
66
a) Ali je velikost odprtine vpliva na povprečni odstotek sprosti radona? Uporabi
  0, 05
b) Poišči P-vrednost za F-statistiko
c) Analiza ostankov iz tega poskusa.
d) Poiščite 95% interval zaupanja za povprečno vrednost sproščenega radona,
ko je odprtina premer je 1,40.
J. Šrekl
Numerične metode v varnosti
67
Narava nam ponuja strukture
6
6.1.
Modeli strukturnih enačb
Merjene in latentne spremenljivke
V statistiki se pri preučevanju množičnih pojavov srečujemo z dvema vrstama
spremenljivk. Na eni strani imamo spremenljivke, ki jih lahko opazujemo in merimo z
merilnimi inštrumenti ali štejemo (število požarov na časovno enoto, izpostavljenost,
zaprašenost, količino nevarne snovi v prostoru, število zaposlenih, delovni čas,
osvetljenost, itd) ali pa s pomočjo vprašalnikov (znanje, pričakovana reakcija na
določen dogodek, odnos do posameznih pojavov, izobrazba, kvalifikacija, itd.).
Spremenljivke lahko merimo s poljubnimi fizikalnimi ali drugačnimi enotami
(točkovanje). Pri tem ločimo tri tipe spremenljivk [4]:
-
-
številske (numerične)
o zvezne spremenljivke (vse realne vrednosti iz nekega intervala)
o diskretne spremenljivke (posamezne realne ali cele vrednosti)
urejenostne (ordinalne) spremenljivke (vrednosti omogočajo kvečjemu
ureditev enot po velikosti npr., ocena čistoče, ocena vzdrževanja naprav);
imenske (nominalne) spremenljivke (vrednosti omogočajo le razlikovanje
z enakostjo ali neenakostjo med seboj, npr. vrsta dejavnosti);
razmernostne spremenljivke (vrednosti omogočajo tudi primerjavo
razmerij med vrednostma dvojic).
Kadar govorimo o enem samem pojavu, ki je odvisen od ene same spremenljivke je
pojav opisan, če lahko najdemo povezavo med opazovano spremenljivko in pojavom.
Take spremenljivke imenujemo merjene spremenljivke, indikatorji ali kazalci. Stanje,
ki vpliva na te kazalce imenujemo latentna spremenljivka. Vsako stanje z neznanimi
odvisnostmi vpliva na svoje indikatorje.
Primer 6.1.:
Na požarno varnost v objektu vplivajo vgrajeni varovalni sistemi, poučenost o
nevarnostih, usposobljenost za reakcije ob nevarnostih, itd. Vsakega od teh
indikatorjev lahko na delovnem mestu izmerimo, vendar nam te količine še nič ne
povedo o stopnji varnosti v objektu. Presoja stopnje požarne varnosti prešteva zgolj
prekoračitve sprejemljivih vrednosti indikatorjev in ne upošteva moči vpliva
posameznega indikatorja. Spremenljivke, ki spremljajo analizo varnosti, sestavljajo
sistem, ki ga je potrebno obravnavati kot celoto in ne zgolj kot singularnosti
posameznih spremenljivk.
J. Šrekl
6.2.
Numerične metode v varnosti
68
Konstrukcija modelov regresij [5]
Bivariantna regresija
Imamo dve spremenljivki, indikator merljivo neodvisno spremenljivko X 1 (vrednost x)
in njeno posledico odvisno spremenljivko Y1 (vrednost y).
X1
Y1
Y
D1
E
1
Slika
6.1. Shema enostavne regresije
Puščica med obema spremenljivkama predstavlja regresijsko odvisnost. Z enačbo
lahko zapišemo preprosto zvezo:
y  0  1 x  E .
(6.1)
Bivariantno regresijo je mogoče razširiti do modela tako, da med merjeno in odvisno
spremenljivko vrinemo latentno spremenljivko F1 , ki je ne merimo. Povezave kaže
diagram povezav ali diagram poti. Varianca neodvisne spremenljivke X 1 pojasnjuje
latentno spremenljivko in napake merjenja:
1
E
0
Slika 6.2. Shema enostavnega modela
Čeprav diagram vsebuje latentno spremenljivko, je v osnovi identičen prejšnjemu.
Predpostavili smo, da je napaka merjenja enaka 0, kar ima za posledico, da sta
latentna spremenljivka in indikator enaka. To pomeni, da je varianca indikatorja v
celoti (100%) pojasnjena z latentno spremenljivko. Če velja enostavna regresija, je
mogoče sklepati, da velja tudi model povezav. (Npr. Naj bi veljalo, da je velikost
požara odvisna samo od količine goriva.) Običajno ni mogoče spremenljivke določati
z enim samim merjenim indikatorjem. (Primer: Požarno obremenitev stavbe ni
mogoče meriti zgolj s količino vgrajenih gorljivih snovi. Različne snovi in oblike
različno vplivajo na možnost požara.) Poleg tega pri merjenju nastajajo napake. Zato
enostavna regresija v praksi proučevanja požarnega tveganja ne velja in je potrebno
uporabiti multivariantno analizo podatkov.
Multiple regresije
Če več indikatorjev določa neko odvisno spremenljivko, potem govorimo o multipli
regresiji. Diagram zarišemo:
J. Šrekl
Numerične metode v varnosti
69
E
Y
Slika 6.3. Shema multiple regresije
Enačba za tako multipla (mnogoterna) regresijo se zapiše:
y  0  1 x1  2 x2  3 x3  E
(6.2)
Model vsebuje spet zgolj merjene spremenljivke, ki pa se med sabo povezujejo. Če
odvisno ali pojasnjeno spremenljivko, ki jo pojasnjujejo indikatorji od X 1 do X 3
nadomestimo s konstruktom (latentno spremenljivko) F1 , se diagram povezav
spremeni. S tem ne pojasnjujemo več odvisne spremenljivke, ampak iščemo vplive
skrite latentne spremenljivke na merjene indikatorje. Ne pričakujem, da bomo lahko
izmerili indikatorje brez napake, zato se v modelu pojavijo tudi napake merjenja
( E1 , E2 , E3 ).
E1
E2
E3
Slika 6.4. Shema merskega modela
Latentne spremenljivke, imenovane tudi faktorji, dobimo v ne merjenem konstruktu
ali latentnem konstruktu. Tak konstrukt smo že srečali pri eni sami neodvisni
spremenljivki (glej sliko 3.2). Če je neodvisnih spremenljivk več, dobimo strukturni
model (SEM), kjer konstruiramo latentne spremenljivke s pomočjo potrjevalne
faktorske analize (confirmatory factor analysis – CFA). Spremenljivke od
E1 , E2 , E3
predstavljajo napake merjenja, X 1 , X 2 in X 3 so merjeni indikatorji in F1 je latentna
spremenljivka.
V čem se model na sliki 6.4. razlikuje od modela na sliki 6.3? Prvi pojasnjuje vplive
posameznih indikatorjev na merjeno odvisno spremenljivko (v kakšnih deležih
vplivajo posamezni indikatorji na vrednost odvisne spremenljivke), drugi model pa
pojasnjuje vplive skrite (neopazovane, latentne) spremenljivke na posamezne
opazovane indikatorje. (Visoka požarna ogroženost vpliva na večje vrednosti
J. Šrekl
Numerične metode v varnosti
70
indikatorjev požarne ogroženosti.) Diagram lahko razširimo na več latentnih
spremenljivk, ki so odvisne od merjenih indikatorjev.
6.3.
Strukturni in merski model [5]
Strukturni model (SEM)
Analiza povezav je približevanje modeliranja pojasnjenih povezav med opazovanimi
spremenljivkami. Vse spremenljivke (neodvisne, odvisne) so merjene. Oglejmo si
primer sistema s tremi neodvisnimi in dvema odvisnima spremenljivkama.
Spremenljivke so standardizirane:
y1   11 x1   12 x 2   13 x3  E1
(6.4)
y 2   21 x1   22 x2   23 x3  E2
γ11
x1
γ12
γ21
E1*
y2
E2*
γ22
x2
γ13
y1
γ23
x3
Slika 6.5. Strukturni model opazovanih spremenljivk
x1
x3 so neodvisne spremenljivke, y1 in y 2 sta odvisni spremenljivki.
Sistem ima šest parametrov  11 do  23 , ki jih želimo oceniti. Ti parametri določajo
do
delne regresijske koeficiente in pojasnjujejo povezave med neodvisnimi in odvisnimi
E1 , E2 v enačbah predstavljajo
spremenljivkami. Zadnji dve spremenljivki
residualni varianci ali varianci napak. To niso zgolj napake merjenja, ampak tudi
vpliv medsebojnega delovanja posameznih spremenljivk. Sistem (5.4) lahko splošno
zapišemo v matrični obliki:
Y = ΓX + E .
(6.5)
Y in X sta vektorja odvisnih in neodvisnih spremenljivk, Γ je matrika delnih
regresijskih koeficientov. Tu gre za posplošitev multiple regresije, kjer imamo
namesto ene odvisne spremenljivke lahko več odvisnih spremenljivk, kjer so
medsebojni vplivi povezani preko indikatorjev.
Merski model
J. Šrekl
Numerične metode v varnosti
71
Merski model pomeni pretvorbo hipotetičnega konstrukta povezave latentnih
spremenljivk z opazovanimi indikatorji, ki jih merimo, v matematični jezik.
Yi so
merjene odvisne spremenljivke (indikatorji) in Fi so latentne odvisne spremenljivke,
ki vplivajo na indikatorje. Obstaja tudi medsebojni vpliv latentnih spremenljivk.
Model je nadgradnja faktorske analize, kjer poleg vplivov latentnih spremenljivk na
indikatorje pojasnjuje tudi medsebojne odvisnosti latentnih spremenljivk. Zato
tovrstni merski model pogosto imenujemo tudi pojasnjevalna faktorska analiza (ang.
confirmatory factor analisys) – CFA. Z diagramom zapišemo model odvisnosti, kjer
z  i označujemo faktorske uteži (ang. factor loading) in z
faktorje ali napake merjenja.
Ei označimo specifične
Y1
E1*
Y2
E2*
λ3*
Y3
E 3*
λ4*
Y4
E 4*
Y5
E5*
Y6
E 6*
λ1*
λ2*
F1
*
λ5*
F2
λ6*
λ7*
Y7
E 7*
Slika 6.6. Merski model z dvema latentnima spremenljivkama
Neznane količine v modelu smo označili z zvezdicami. Na podlagi diagrama
zapišemo sistem enačb:
J. Šrekl
Numerične metode v varnosti
72
Y1  1 F1  E1
Y2   2 F1  E2
Y3   3 F1  E3
Y4   4 F2  E4
(6.6)
Y5   5 F2  E5
Y6   6 F2  E6
Y7   7 F2  E7
Neznane količine v modelu imenujemo parametre modela. Za njihovo določanje
veljajo naslednja pravila [5]:
Pravilo 1. Vse variance odvisnih spremenljivk so parametri modela.
Pravilo 2. Vse kovariance med odvisnimi spremenljivkami so parametri modela.
Pravilo 3. Vse faktorske uteži, ki povezujejo latentne spremenljivke z indikatorji so
parametri modela.
Pravilo 4. Vsi regresijski koeficienti med indikatorji in latentnimi spremenljivkami so
parametri modela.
Pravilo 5. Za vsako latentno spremenljivko je treba določiti mersko lestvico.
Pravilo 5 je potrebno, ker v samem sistemu niso določene merske lestvice za
spremenljivke, le te so odvisne so od različnih indikatorjev z različnimi merskimi
lestvicami.
Lastnosti računanja z variancami in s kovariancami lahko združimo štiri pravila.
Naj bodo X, Y, Z in U slučajne spremenljivke in a, b, c in d realna števila. Potem
velja:
(1) Cov( X , X )  Var ( X ),
(2) Cov( a  X  b  Y , c  Z  d  U ) 
 ac  Cov( X , Z )  bc  Cov(Y , Z )  ad  Cov( X ,U )  bd  Cov(Y ,U ),
(3) Var (a  X  b  Y )  a2Var ( X )  b2Var (Y )  2ab  Cov( X ,Y ),
(4) Var (a  X  b  Y )  a2Var ( X )  b2Var (Y ),
če sta spremenljivki X in Y nekorelirani.
Model vsebuje indikatorje, ki so med sabo povezani ali korelirani. Zato izračunamo
variančno-kovariančno matriko modela, ki nam daje mero povezanosti med
indikatorji. Ker je variančno-kovariančna matrika Σ simetrična (simetrični elementi
glede na glavno diagonalo so enaki), predstavljajo vrednosti v trikotniku nad glavno
diagonalo odvečne količine (ang. redundant informations). Če je p indikatorjev, je
matrika dimenzije p×p in vsebuje
p( p  1)
potrebnih ali neredudantnih vrednosti.
2
Matrika je spodnja trikotna in ima za elemente kovariance med spremenljivkami. Po
diagonali zaradi lastnosti (1) dobimo variance spremenljivk.
 Var ( Y1 )
Cov( Y ,Y ) Var ( Y )
1
2
2
Σ
 Cov( Y1,Y3 ) Cov( Y2 ,Y3 ) Var ( Y3 )










(6.7)
J. Šrekl
Numerične metode v varnosti
73
Po pravilu 5 je
Var ( Fi )  1 ,
iz česar sledi, da je po (2):
Cov(Yi , Yj )  Cov(i Fk  E i ,  jFk  E j ) 
  i   j  Cov(Fk , Fk )   i  Cov(Fk , E j )   j  Cov(E i , Fk )  Cov(E i , E j ) 
(6.8)
  i   j  Cov(Fk , Fk )   i   j  Var(Fk , Fk ) 
 i   j
za vse indikatorje, ki so odvisni od iste latentne spremenljivke
sta odvisna od različnih latentnih spremenljivk pa velja:
Fk . Za indikatorja ki
Cov(Yi , Yj )  Cov(i F1  E i ,  jF2  E j ) 
 i   j  Cov(F1 , F2 )   i  Cov(F1 , E j )   j  Cov(E i , F2 )  Cov(E i , E j ) 
(6.9)
 i   j  Cov(F1 , F2 )   i   j  12
pri čemer je Cov( Fi , F j )  ij in
Var( Fi )  1 .
Izračunajmo še varianco indikatorjev, kjer upoštevamo
Var ( Ei )   i in dobimo:
Var(Yi )  Cov(i Fk  E i ,  i Fk  E i ) 
 i2  Cov(Fk , Fk )   i  Cov(Fk , E i )   i  Cov(E i , Fk )  Cov(E i , E i ) 
 i2  Var(Fk , Fk )  Var(E i , E i ) 
(6.10)
 i2  i
Matriko varianc in kovrianc  lahko izpišemo:
 12   1

 22   2
 1  2
 1  3
 2 3
 32   3


 2  412  3  412  24   4
Σ   1 4 12
 1  512  2  512  3  512  4  5

 1  612  2  612  3  612  4  6
 1  7 12  2  7 12  3  7 12  4  7


 52   5
5 6
5 7
 26   6
6 7
 27   7





 (6.11)






Indikatorje izmerimo v vzorcu velikosti n, ki vsebuje n izmerjenih vrednosti za vsak
indikator. Iz vrednosti v vzorcu izračunamo vzorčne variance in vzorčne kovariance.
Sestavimo vzorčno variančno-kovariančno matriko, ki jo označimo z S. Z
nadomestitvijo matrike
Σ z matriko S dobimo sistem
p( p  1)
2
nelinearnih
enačb za (2p+r) neznank – parametrov modela (p faktorskih uteži, p residualnih
J. Šrekl
Numerične metode v varnosti
74
varianc in r kovarianc latentnih spremenljivk). Za rešljivost sistema mora biti
izpolnjen pogoj:
p( p  1)
 2 p  r , ali
2
p( p  3)
 r.
2
(6.12)
(6.13)
Iz enačbe sledi, da morajo biti več kot trije indikatorji. Razlika med številom enačb in
številom neznank je prostostna stopnja modela. Če ima model manj enačb kot
neznank, ni določen. Minimalne zahteve sta dva indikatorja na eno latentno
spremenljivko. Če ima model enako število enačb kot je neznank (parametrov
modela), ga imenujemo saturiran ali natančno določljiv model. V splošnem računamo
s sistemi, kjer nastopa določeno število prostostnih stopenj. Tak model ni natančno
določljiv, lahko pa izračunamo prileganje (ang. fit) rešitev na model. Iščemo torej
najboljšo rešitev modela, ki jo je mogoče iz danih podatkov izračunati. Sistem enačb
rešujemo na primer z metodo najmanjših kvadratov.
Potrjevalna faktorska analiza nam izračuna najverjetnejše vrednosti iz matrike Σ.
Govorimo o verjetnih rezultatih, ker ti rezultati temeljijo na izbranem statističnem
vzorcu. Dobimo torej matriko za kovariance med latentnimi spremenljivkami
1
Φ
12
matriko faktorskih uteži
in matriko residualnih varianc
0
,
1
(6.14)
Λ  1 2  7 
(6.15)
Θ  1  2   7 .
(6.16)
CFA nam oceni relacije med konstrukti. Nobena specifična zveza ni predpostavljena
med konstrukti, iščemo le korelacije med njimi. Pogosto srečujemo modele, kjer že
vnaprej predpostavljamo specifično pojasnjevalno zvezo med konstrukti. Za računanje
takih zvez je potrebno sestaviti strukturni regresijski model.
Primer 6.2.: Računanja pojasnjevalne faktorske analize.
Želimo pokazati zvezo med latentnimi spremenljivkami: SAMOSPOŠTOVANJE,
DEPRESIVNOST in IMPULZIVNOST
 Podatki vsebuje 204 meritev 12 zveznih spremenljivk.
 12 spremenljivk so kazalniki latentnih spremenljivk samospoštovanje,
depresija in impulzivnost.
 12 spremenljivk je merjeno s 5-točkovno Likertovo lestvico.
Oblika merjenja je tipičnih Likertova 5-točkovna lestvica:
Na vprašanje ali se strinjam s ponujeno trditvijo je lestvica petih odgovorov
1. se absolutno ne strinjam
2. se delno ne strinjam
J. Šrekl
Numerične metode v varnosti
75
3. ne vem (nedoločno)
4. se strinjam
5. se močno strinjam
Podatke dobimo v datoteki depres.PSF znotraj programa LISREL v mapi TUTORIAL
SELF1 za SELF5 so kazalniki latentne spremenljivke samospoštovanje (selfest).
DEPRES1 za DEPRES4 so kazalniki depresije latentne spremenljivke (depress).
IMPULS1 za IMPULS3 so kazalniki impulzivnosti latentne spremenljivke (impulz).
Sestavimo ustrezni shematski diagram poti in povezav med spremenljivkami kot jih
želimo povezati. Diagram poti – povezav (angl. path diagram) modela je prikazan na
naslednji sliki.
Do diagrama bomo prišli z ustreznimi definicijami spremenljivk v podoknu „Label”.
Na levi imamo faktorje (indikatorje), na desni latentne spremenljivke.
Analiza:
J. Šrekl
Numerične metode v varnosti
76
1. Izberite novo možnost v meniju File (Datoteka) za nalaganje New dialog box
(Novo pogovorno okno).
2. Izberite možnost Path diagram (Diagram poti) s seznama polja za novo
pogovorno okno za nalaganje pogovornem oknu Save As ( Shrani kot).
3. Vnesite ime depressSEM.pth na področju File name (Ime datoteke) niza.
4. Kliknite na gumb Save (Shrani), da odprete prazno okno PTH.
5. Izberite Title and commments (Naslov in Komentarji) v meniju Setup za
odprtje okna.
6. Vnesite naslov A Model for Self-Esteem v področje naslova nizov.
7. Kliknite na gumb Next (Naprej) za nalaganje v pogovorno okno Group
Names( imena skupine).
8. Kliknite na gumb Next (Naprej) za nalaganje v pogovorno okno Labels
(imena spremenljivk).
9. Kliknite na Add/Read Variables (Dodaj / Preberi spremenljivke) gumb za
nalaganje v pogovorno okno.
10. Izberite File System PRELIS možnost Read from file: na spustnem seznamu
polja branje iz datoteke.
11. Kliknite na gumb Browse (Prebrskaj) za nalaganje Prebrskaj okno.
12. Izberite datoteko depress.PSF v mapi TUTORIAL.
13. Kliknite na gumb Open, da se vrnete na Add / Read spremenljivke v
pogovornem oknu.
14. Kliknite na gumb OK (V redu), da se vrnete v pogovorno okno za imena
spremenljivk.
15. Kliknite na gumb Add Latent Variables ( Dodaj Latentna Spremenljivke), da
lahko zapišemo imena latentnih spremenljivk.
16. Vpišite okolje selfest za samospoštovanja v polje niza področju.
17. Kliknite na gumb OK (V redu), da se vrnete v pogovorno okno.
18. Kliknite na gumb Add Latent Variables ( Dodaj Latentna Spremenljivke), da
lahko zapišemo imena latentnih spremenljivk..
19. Vpišite oznak impuls za impulzivnost v nizu področju.
20. Kliknite na gumb OK (V redu), da se vrnete v pogovorno okno.
21. Kliknite na gumb Add Latent Variables ( Dodaj Latentna Spremenljivke), da
lahko zapišemo imena latentnih spremenljivk..
22. Vpišite oznak depress za impulzivnost v nizu področju.
Kliknite na gumb V redu, da se vrnete v pogovorno okno nalepke.
J. Šrekl
Numerične metode v varnosti
77
Po ustrezni definiciji spremenljivk in procesa računanja v „diagram” izrišemo
povezave, tako, da prenesemo vse spremenljivke in zarišemo povezovalne puščice.
Kliknite na gumb V redu, da se vrnete v okno za PTH depress SEM.PTH.
Označite polja Y za SELF1 - SELF5.
Označite polja Y za IMPULS1 - IMPULS3.
Označite Eta kvadratek selfest in impuls.
Kliknite imena spremenljivk, povlecite in spustite v PTH okna, kot kaže slika
diagrama
Izberite ikono puščice v orodni vrstici Risanje.
Kliknite in povlecite 5 indikator poti iz latentne spremenljivke selfest v SELF1,
SELF2, SELF3, SELF4 in SELF5.
Kliknite in povlecite 4 indikator poti iz latentne spremenljivke depress v
DEPRES1, DEPRES2, DEPRES3 in DEPRES4.
Kliknite in povlecite indikator 3 poti iz latentnih spremenljivk impulz v IMPULS1,
IMPULS2 in IMPULS3.
Kliknite in povlecite 2 strukturne poti iz latentne spremenljivke depress v selfest in
impuls.
Kliknite in povlecite 1 strukturno pot iz latentnih spremenljivk impulz v selfest,
tako, da dobimo okno PTH na sliki
J. Šrekl
Numerične metode v varnosti
78
Sam program se nam bo shranil v *.SPJ. Izberi opcijo Build SIMPLIS
Syntax v Setup menu za odpiranje naslednjega SPJ okna.
Z zagonom RUN Lisrel bomo pognali iteracijski program računanja, ki nam bo v
nekaj korakih izračunal povezave med posameznimi spremenljivkami. Povezave so
lahko kovariance ali pa korelacije v standardizirani obliki. Povedo kako latentna
spremenljivka vpliva na posamezni idikator, kar je mogoče izračunati tudi z običajno
faktorko analizo. Pri računanju z LISREL pa dobimo tudi medsebojne odvisnosti
latentnih spremenljivk (glej vrednosti na diagramu. Dodane so tudi usnovne mere
zanesljivosti rezultatov o katerih bomo govorili kasneje.
Poženi Run LISREL ikono na glavnem meniju, da ustvariš naslednje PTH
okno.
J. Šrekl
Numerične metode v varnosti
79
LISREL model
Model strukturnih enačb kot ga je postavil Jöreskog [7], je sestavljen iz neodvisnih in
odvisnih indikatorjev ter neodvisnih in odvisnih latentnih spremenljivk. Metode
proučevanja SEM modelov temeljijo na pojasnjevanju povezav med latentnimi
spremenljivkami, ki so merjene z indikatorji.
Multipla regresija vsebuje odvisnost med endogenimi spremenljivkami in eksogenimi
faktorji, ne upošteva pa vplivov povezav med faktorji. Te povezave upošteva
modeliranje s strukturnimi enačbami (SEM) s pomočjo modela LISREL (linear
structural relationship), ki upošteva tudi vse notranje medsebojne vplive faktorjev.
LISREL metoda ocenjuje neznane koeficiente množice linearnih strukturnih enačb
[8]. LISREL model so konstruirali v zgodnjih sedemdesetih letih Keesling, Jöreskog
in Wiley [9] . Model je razširil in dopolnil Sörbom (1974) [10]. (Programsko opremo
za računanje LISREL modelov so razvili v Scientific Software International, Inc.
7383 N. Lincoln Avenue, Suite 100 Lincolnwood, IL 60712, U.S.A.)
Popolni LISREL model za slučajni vzorec je sestavljen iz treh osnovnih komponent:
strukturnega modela in dveh merskih modelov. Merska modela pojasnjujeta zveze
med indikatorji in latentnimi spremenljivkami za endogeni in eksogeni del modela.
Strukturni model pa pojasnjuje povezave med endogenimi in eksogenimi latentnimi
spremenljivkami. Daje nam torej zvezo med neodvisnimi in odvisnimi
spremenljivkami tako, da upošteva tudi zveze med neodvisnimi spremenljivkami.
J. Šrekl
Numerične metode v varnosti
80
X1
Y1
X2
Y2
X3
Y3
X4
Y4
X5
Merski model
Strukturni model
Merski model
Slika 6.7. Strukturni model LISREL
Model opišemo s tremi matričnimi enačbami. Strukturni model opišemo z enačbo:
(6.17)
η  Bη  Γξ  ζ ,
in oba merska modela z enačbama:
y  Λy η  ε
(6.18)
x  Λ xξ  δ .
(6.19)
V enačbah sta η in ξ vektorja latentnih spremenljivk; x in y vektorja opazovanih
(observed) spremenljivk; ε in δ sta vektorja merskih napak; in ζ je vektor
strukturnih napak. y in η sta vektorja endogenih spremenljivk ( spremenljivke, ki so
pojasnljive s spremenljivkami obravnavanega modela; odvisne spremenljivke); ξ in x
sta vektorja eksogenih
spremenljivk (spremenljivke, ki so pojasnljive s
spremenljivkami zunaj obravnavanega modela, neodvisne spremenljivke); ε , δ , in ζ
so eksogene spremenljivke.
B (beta, BE) je matrika regresijskih koeficientov za η , Γ (gama, GA) je matrika
regresijskih koeficientov med η in ξ . Λ y in Λ x sta matriki faktorskih uteži med y in η
ter x in ξ .
Iz enačbe (2.13) dobimo:
η  (1  B)1 Γξ  (1  B)1ζ
(6.20)
Primer 6.3.: Računanja LISREL modela
Želimo pokazati kako neodvisna latentna spremenljivka DEPRESIVNOST vpliva na
latentni spremenljivki SAMOSPOŠTOVANJE in IMPULZIVNOST
 Podatki vsebuje 204 meritev 12 zveznih spremenljivk.
J. Šrekl
Numerične metode v varnosti
81
 12 spremenljivk so kazalniki latentnih spremenljivk samospoštovanje,
depresija in impulzivnost.
 12 spremenljivk je merjeno s 5-točkovno Likertovo lestvico.
Oblika merjenja je tipična Likertova 5-točkovna lestvica:
Na vprašanje ali se strinjam s ponujeno trditvijo je lestvica petih odgovorov
1. se absolutno ne strinjam
2. se delno ne strinjam
3. ne vem (nedoločno)
4. se strinjam
5. se močno strinjam
Podatke dobimo v datoteki depres.PSF znotraj programa LISREL v mapi TUTORIAL
SELF1 SELF5 so indikatorji za SAMOSPOŠTOVANJE (selfest).
DEPRES1 - DEPRES4 so indikatorji za DEPRESIVNOST (depress).
IMPULS1 - IMPULS3 so indikatorji za IMPULZIVNOST (impuls).
Poskušamo ugotoviti kako DEPRESIVNOST vpliva na IMPULZIVNOST in
SAMOSPOŠTOVANJE
Diagram poti – povezav (angl. path diagram) modela je prikazan na
naslednji sliki.
J. Šrekl
Numerične metode v varnosti
82
Analiza:
1. Izberite novo možnost v meniju File (Datoteka) za nalaganje New dialog box
(Novo pogovorno okno).
2. Izberite možnost Path diagram (Diagram poti) s seznama polja za novo
pogovorno okno za nalaganje pogovornem oknu Save As ( Shrani kot).
3. Vnesite ime depressSEM.pth na področju File name (Ime datoteke) niza.
4. Kliknite na gumb Save (Shrani), da odprete prazno okno PTH.
5. Izberite Title and commments (Naslov in Komentarji) v meniju Setup za
odprtje okna.
6. Vnesite naslov A Model for Self-Esteem v področje naslova nizov.
7. Kliknite na gumb Next (Naprej) za nalaganje v pogovorno okno Group
Names( imena skupine).
8. Kliknite na gumb Next (Naprej) za nalaganje v pogovorno okno Labels
(imena spremenljivk).
9. Kliknite na Add/Read Variables (Dodaj / Preberi spremenljivke) gumb za
nalaganje v pogovorno okno.
10. Izberite File System PRELIS možnost Read from file: na spustnem seznamu
polja branje iz datoteke.
11. Kliknite na gumb Browse (Prebrskaj) za nalaganje Prebrskaj okno.
12. Izberite datoteko depress.PSF v mapi TUTORIAL.
13. Kliknite na gumb Open, da se vrnete na Add / Read spremenljivke v
pogovornem oknu.
14. Kliknite na gumb OK (V redu), da se vrnete v pogovorno okno za imena
spremenljivk.
15. Kliknite na gumb Add Latent Variables ( Dodaj Latentna Spremenljivke), da
lahko zapišemo imena latentnih spremenljivk.
16. Vpišite okolje selfest za samospoštovanja v polje niza področju.
17. Kliknite na gumb OK (V redu), da se vrnete v pogovorno okno.
18. Kliknite na gumb Add Latent Variables ( Dodaj Latentna Spremenljivke), da
lahko zapišemo imena latentnih spremenljivk..
19. Vpišite oznak impuls za impulzivnost v nizu področju.
20. Kliknite na gumb OK (V redu), da se vrnete v pogovorno okno.
21. Kliknite na gumb Add Latent Variables ( Dodaj Latentna Spremenljivke), da
lahko zapišemo imena latentnih spremenljivk..
22. Vpišite oznak depress za impulzivnost v nizu področju.
23. Kliknite na gumb V redu, da se vrnete v pogovorno okno nalepke.
J. Šrekl
Numerične metode v varnosti
83
Kliknite na gumb V redu, da se vrnete v okno za PTH depressSEM.PTH.
Označite polja Y za SELF1 - SELF5.
Označite polja Y za IMPULS1 - IMPULS3.
Označite Eta kvadratek selfest in impuls.
Kliknite imena spremenljivk, povlecite in spustite v PTH okna, kot kaže slika
diagrama
Izberite ikono puščice v orodni vrstici Risanje.
Kliknite in povlecite 5 indikator poti iz latentne spremenljivke selfest v SELF1,
SELF2, SELF3, SELF4 in SELF5.
Kliknite in povlecite 4 indikator poti iz latentne spremenljivke depress v
DEPRES1, DEPRES2, DEPRES3 in DEPRES4.
Kliknite in povlecite indikator 3 poti iz latentnih spremenljivk impulz v IMPULS1,
IMPULS2 in IMPULS3.
Kliknite in povlecite 2 strukturne poti iz latentne spremenljivke depress v selfest in
impuls.
Kliknite in povlecite 1 strukturno pot iz latentnih spremenljivk impulz v selfest,
tako, da dobimo okno PTH na sliki
J. Šrekl
Numerične metode v varnosti
84
Izberi opcijo Build SIMPLIS Syntax v Setup menu za odpiranje naslednjega
SPJ okna.
Poženi Run LISREL ikono na glavnem meniju, da ustvariš naslednje PTH
okno.
J. Šrekl
Numerične metode v varnosti
85
Iz diagrama, vidimo da depresija veliko bolj vpliva na samopodobo kot na
impulzivnost. Med samopodobo in impulzivnostjo ni značilne statistične povezave. S
povečano depresivnostjo se slabša tudi samopodoba človeka in s tem povezani
kazalniki (faktorji, indikatorji) samopodobe.
Ni pa nujno, da v modelu nastopajo vsi tipi spremenljivk. V tem primeru govorimo o
podmodelih LISREL modela (angl. submodels). Pri manjšem številu indikatorjev
iščemo povezavo direktno med eksogenimi indikatorji in endogenimi latentnimi
spremenljivkami. Indikatorji prevzamejo vlogo latentnih spremenljivk in dobimo
model, ki je sestavljen iz enega strukturnega in enega merskega modela.
Y1
X1
Y2
Y3
X2
Y4
Strukturni model
Merski model
J. Šrekl
Numerične metode v varnosti
86
Slika 6.8. Strukturni model LISREL brez eksogenega merskega modela
Tak model opišemo samo z dvema matričnima enačbama. Namesto enačbe (6.20)
zapišemo enačbo:
(6.21)
η  Bη  Γx  ζ
in enačbo merskega modela (2.15). Seveda se tudi enačba (2.16) ustrezno spremeni:
η = (1- B)-1 Γx + (1- B)-1 ζ
(6.22)
V primeru, ko imamo eno samo endogeno latentno spremenljivko, lahko model še
nekoliko poenostavimo.
Y1
X1
Y2
X2
Strukturni model
Merski model
Slika 6.9. Strukturni model LISREL z eno samo latentno spremenljivko
Iz enačbe (3.18) pa lahko dobimo navadno enačbo:
η = Γx + ζ
6.4.
(6.23)
Skladnost statistike
V LISREL modelu računanje parametrov modela spremlja tudi preverjanje skladnosti
modela oziroma vzorčna potrditev modela. Skladnost merimo z različnimi indeksi,
LISREL program preverja 15 indeksov, nekateri drugi programi za računanje SEM
modela celo več. Oglejmo si nekaj najpomembnejših indeksov za preverjanje
skladnosti modela. Računanje mer skladnosti temelji na prilagajanju parametrov
modela momentom vzorca, kar pomeni primerjava opazovane oziroma merjene
kovariančne matrike z ocenjeno [11]. Pri tem predpostavljamo, da testirani model
J. Šrekl
Numerične metode v varnosti
87
zares velja. Mere skladnosti pogosto imenujemo tudi funkcije diskrepance (funkcija
neskladnosti).
Hi-kvadrat mero modela (model chi-square) pogosto imenujemo tudi diskrepanca
modela in je eden najenostavnejših testov skladnosti modela. Hi-kvadrat ima
neznačilno (not significant) vrednost, če je skladnost dobra. Značilnost te številke
opozarja na slabo skladnost modela z vzorčnimi vrednostmi, kar pomeni značilne
razlike med vzorčnimi kovariancami in kovariancami v modelu. Zato običajno ne
govorimo o zavračanju ali sprejemanju modela na podlagi vrednosti  2 , ampak samo
o boljšem ali slabšem prilagajanju testa glede na vrednost  2 [7]. LISREL hi-kvadrat
modela preprosto imenuje hi-kvadrat (sinonimi: chi-square fit index, chi-squre
goodness of fit). Okrajšano ga označujemo s CFI.
Relativni hi-kvadrat ali normiran hi-kvadrat (normed chi-square fit index – NCFI) je
CFI deljen s prostostno stopnjo in s tem postane manj občutljiv na velikost vzorca.
Različni avtorji zahtevajo različno mejno vrednost NCFI za sprejemljivost modela.
Nekateri trdijo, da je vrednost večja od 5 povsem ustrezna vrednost za sprejemanje
modela, drugi sprejemajo model pri 2 ali manj.
Indeks skladnosti (goodness of fit index – GFI) je razmerje med fit funkcijo modela in
fit funkcijo modela, kjer so vsi parametri enaki 0. Model sprejmemo, če je GFI > 0.90.
Urejen indeks skladnosti (adjust goodness of fit index – AGFI) je izpeljan indeks iz
GFI. Pri majhnih vzorcih (n < 200) ga lahko direktno izračunamo:
AGFI  1 
(1  GFI )  p  ( p  1)
2  df
(3.20)
Pri čemer je p število parametrov v modelu in df število prostostnih stopenj. Model
bo sprejemljiv, če bo vrednost AGFI večja od 0.9.
Primerjalni fit indeks (comparative fit index – CFI) imenovan tudi Bentlerjev CFI
primerja obstoječi model z ničelnim modelom (kjer so vse spremenljivke
nekorelirane). Vrednosti CFI so med 0 in 1, vrednosti bliže 1 pomenijo boljšo
prilagoditev.
Normiran fit indeks (normed fit index – NFI) se je razvil kot alternativa CFI.
Zavzema vrednosti med 0 in 1, pri čemer pomeni 1 popolno prilagoditev modela. Za
vrednosti med 0.90 in 0.95 je model sprejemljiv, za vrednosti nad 0.95 zelo dober, pri
vrednostih pod 0.90 pa je treba model popraviti.
Ne normirani fit indeks (not normed fit index – NNFI) se izračuna s formulo:
NNFI 

chisqn
chisq
dfn
df
chisqn
dfn
1
,
(3.21)
kjer sta chisqn in dfn hi-kvadrat statistika in število prostostnih stopenj za model in
chisq in df hi-kvadrat statistika in število prostostnih stopenj za ničelni model (vse
spremenljivke nekorelirane).
J. Šrekl
Numerične metode v varnosti
88
Koren srednje vrednosti kvadratov ostankov (root mean square residual – RMSR ali
RMR) je koren iz srednje vrednosti kvadratov razlik med vzorčnimi in ocenjenimi
kovariancami).
Stasndardiziran RMSR (SRMS) je izračunan RMSR nad standardiziranimi razlikami.
Čim manjši je SRMS, tem večja je skladnost modela z vzorcem.
J. Šrekl
Numerične metode v varnosti
7
[1]
[2]
[3]
[4]
[5]
[6]
[7]
[8]
[9]
[10]
[11]
[12]
[13]
[14]
89
Literatura
Boyce, E.W., DiPrima, R.C., Elementary Differential Equations and
Boundary Value Problems, (7th edition), Wiley, New York, 2001
Kreyszig, E., Advance Engineering Mathematics (9th edition), Wiley, New
York, 2006
Montgomery, D. C., Runger, G.C.,Applied Statistics and Probability for
Engineers, (third edition), Wiley, New York, 2003
Dowdy, S., Wearden, S., Chilko, D.,Statistics for Research, Third edition,
Wiley, New Jersey, 2004
Raykov T., Marcoulides G. A., A First Course in Structural Equation
Modeling, LEA, London, 2000,
Pugesek B. H., Concept of structural eqution modeling in biological reserch,
Structural Eqution Modeling, edited by Purgesek B.H. at al., Cambridge
University Press, 2003, 42-59
Jöreskog K. G. A general methods for estimating a linear structural equation
system. Structural Equation Model in the Social Science. New York :
Seminar Press, (1973)
Hershberger S.L., Marcoulides G.A., Parramore M.M., Structural equation
modeling: an introduction; Structural Equation modeling, edited by Purgesek
B.H. at al., Cambridge UP, Cambridge, 2003,. 3-41.
Bentler, P. M., Multivariante analysis with latent variables causual
modeling, Annual Rewiew of Psychology, 31 (1980), 419-456.
Sörbom, D., A general methods for studying differences in factor means and
factor structures between groups, British Journal of Mathematical and
Statistical Psychology, 27 (1974), 229-239.
http://www2.chass.ncsu.edu/garson/pa765/structur.htm
ŠREKL, J., Izbrana poglavja iz matematike. Ljubljana: Fakulteta za kemijo
in kemijsko tehnologijo, Oddelek za tehniško varnost, 1997.
ŠREKL, J., Zapiski predavanj IPMS leto 2005-06. http://www.fkkt.unilj.si/si/msg.php?6116
McKibben, M. A., Kirchner West, J., Excel Manual for Statistics, J. Wiley,
New York, 2005
J. Šrekl
Numerične metode v varnosti
90
Razlogi za pisanje in objavo zapiskov
Zapiski predavanj so namenjeni študentom druge stopnje bolonjskega študija
Tehniške varnosti. To ni učbenik, ki bi razjasnjeval vse teoretične vidike snovi, ki jo
obravnava. Je zgolj priročnik za boljše razumevanje predavanj, poglobljeno znanje je
potrebno iskati v našteti literaturi iz katere povzemam predavanja.
Namen predavanj je študentom pokazati vso pestro uporabo numeričnih in statističnih
metod posebej na področju njihovega raziskovalnega in strokovnega dela. Z dovolj
primeri je prikazan način delovanja orodij. Uporabljam predvsem priročna
računalniška orodja, ki so večini na voljo, le tam kjer je nujno potrebno se sklicujem
na specialne programe za reševanje problemov.
Doc. dr. Jože Šrekl, univ. dipl. mat.