Výpočet citlivosti parametru modelu
Uvažujme model tvaru:
kde X1, X2, X3 jsou všechny náhodné veličiny s rovnoměrným rozdělením pravděpodobnosti
na intervalu (−π, π). Proveďme globální analýzu citlivosti.
Řešení: Nejprve určíme střední hodnotu a rozptyl řešení Y. To provedeme výpočtem pomocí
systému Maple, viz zdrojový kód Maple na konci. Obdržíme tak:
Pravděpodobnostní rozdělení veličiny Y můžeme odhadnout vygenerováním realizací pro
X1, X2, X3, následným vyhodnocením a zobrazením histogramu (v Maple příkazem
Histogram). Histogram získaný v systému Maple pro M = 105
vyhodnocení modelové rovnice
je uveden na obrázku 5.1.
Obr. 5.1 Histogram vygenerovaných řešení pro model příkladu 5.5.
Podle vztahu (5.8)
vypočítáme hodnoty vlivů prvního řádu jednotlivých parametrů Xi na řešení Y (pro i = 1, 2, 3):
Následně určíme podle vztahu (5.9)
indexy globální citlivosti pro parametry X1, X2, X3:
V tomto jednoduchém příkladu bylo možné získat indexy globální citlivosti analyticky
zavedenými vztahy. Pro názornost vypočítejme nyní indexy globální citlivosti pomocí Soboľovy
metody.
Nejprve musíme vygenerovat dvě skupiny realizací parametrů X1, X2, X3. Nechť každá
skupina obsahuje M = 105
realizací pro každý z parametrů. Zvolíme jednu ze skupin, spočteme
M vyhodnocení modelu Y (pro jednotlivé realizace parametrů X1, X2, X3) a určíme jejich
průměr (tj. střední hodnotu) a rozptyl. Obdržíme1
:
Nyní přistoupíme k výpočtu odhadů vlivů jednotlivých parametrů definovaných vztahem
(5.16) až (5.18).
1
V závorce je uvedená pro porovnání zaokrouhlená hodnota přesné hodnoty V(Y).
Pro výpočet vlivu V1 parametru X1 vypočteme dalších M vyhodnocení modelu Y, a to tak,
že budeme brát realizace parametru X1 z první skupiny (té dříve zvolené) a realizace
parametrů X2, X3 z druhé skupiny na začátku vygenerovaných realizací. Zcela analogicky
vypočteme jiných M vyhodnocení modelu Y jak pro určení vlivu V2 parametru X2, tak pro
určení vlivu V3 parametru X3. Získáme níže uvedené hodnoty:
Po výpočtu vlivů V1, V2, V3 můžeme přistoupit k výpočtu indexů globální citlivosti
prvního řádu podle vztahu (5.9):
Jestliže aplikujeme korekční člen (5.21),
jehož hodnota je přibližně −0,0571. Homma a Saltelli (viz [24]) doporučují korekční člen
(5.21) aplikovat bez absolutní hodnoty2
na jednotlivé rozptyly (pro i1,…,is ∈ {1,2,3},
i1 < … < is) tak, že od rozptylů prvního řádu (Vi) se korekční člen odečte. Pak obdržíme:
Základní zdrojový kód pro výpočet výše uvedených hodnot v Maple:
Výpočet indexů globální citlivosti pomocí systému Maple
M := 10^5:
with(Statistics):
X[1] := RandomVariable(Uniform(-Pi, Pi)):
X[2] := RandomVariable(Uniform(-Pi, Pi)):
X[3] := RandomVariable(Uniform(-Pi, Pi)):
X1s := Sample(X[1], M): # M realizaci veliciny X[1]
X2s := Sample(X[2], M): # M realizaci veliciny X[2]
X3s := Sample(X[3], M): # M realizaci veliciny X[3]
Ys := sin~(X1s) + 7*(sin~(X2s))^~2 + ((1/10)*X3s^~4)*(sin~(X1s)):
f[0] := Mean(Ys): # odhad stredni hodnoty
V[0] := Variance(Ys): # odhad rozptylu
X1s2 := Sample(X[1], M): # "druha skupina" M realizaci veliciny X[1]
X2s2 := Sample(X[2], M): # "druha skupina" M realizaci veliciny X[2]
X3s2 := Sample(X[3], M): # "druha skupina" M realizaci veliciny X[3]
2
Absolutní hodnota byla přidána, aby příslušná věta měla správný matematický význam. Znaménko
členu v absolutní hodnotě je však podstatné!
Ys2 := Vector(M):
# citlivost parametru X[1]
Ys2 := sin~(X1s) + 7*(sin~(X2s2))^~2 + ((1/10)*X3s2^~4)*~(sin~(X1s)):
V[1] := Mean(Ys*~Ys2)-f[0]^2;
S[1] = V[1]/V[0];
# citlivost parametru X[2]
Ys3 := sin~(X1s2) + 7*(sin~(X2s))^~2 + ((1/10)*X3s2^~4)*~(sin~(X1s2)):
V[2] := Mean(Ys*~Ys3)-f[0]^2;
S[2] = V[2]/V[0];
# citlivost parametru X[3]
Ys4 := sin~(X1s2) + 7*(sin~(X2s2))^~2 + ((1/10)*X3s^~4)*~(sin~(X1s2)):
V[3] := Mean(Ys*~Ys4)-f[0]^2;
S[3] = V[3]/V[0];
# korekcni clen
Ysk := sin~(X1s2) + 7*(sin~(X2s2))^~2 + ((1/10)*X3s2^~4)*~(sin~(X1s2)):
V[k] := Mean(Ys*~Ysk)-f[0]^2;