Analiza wrażliwości oparta na wariancjach

Analiza wrażliwości oparta na wariancjach (często nazywana metodą Sobola lub wskaźnikami Sobola , od nazwiska Ilyi M. Sobola ) jest formą globalnej analizy wrażliwości . Działając w probabilistycznych , rozkłada wariancję wyniku modelu lub systemu na ułamki, które można przypisać do danych wejściowych lub zbiorów danych wejściowych. Na przykład, biorąc pod uwagę model z dwoma wejściami i jednym wyjściem, można by stwierdzić, że 70% wariancji wyjściowej jest spowodowane wariancją pierwszego wejścia, 20% wariancją drugiego, a 10% wynika z interakcji między dwa. Te wartości procentowe są bezpośrednio interpretowane jako miary wrażliwości. Miary wrażliwości oparte na wariancji są atrakcyjne, ponieważ mierzą wrażliwość w całej przestrzeni wejściowej (tzn. jest to metoda globalna), radzą sobie z nieliniowymi odpowiedzi i mogą mierzyć efekt interakcji w systemach nieaddytywnych .

Dekompozycja wariancji

Z perspektywy czarnej skrzynki każdy model może być postrzegany jako funkcja Y = f ( X ), gdzie X jest wektorem d niepewnych danych wejściowych modelu { X ₁ , X ₂ , ... X _d }, a Y jest wybranym wyjście modelu jednowymiarowego (należy zauważyć, że to podejście bada wyniki modelu skalarnego, ale wiele wyjść można analizować za pomocą wielu niezależnych analiz wrażliwości). Ponadto zakłada się, że dane wejściowe są niezależnie i równomiernie $_$ ${\ Displaystyle i = 1,2, ..., d$ jednostkowym ja . Nie powoduje to utraty ogólności, ponieważ dowolną przestrzeń wejściową można przekształcić w ten jednostkowy hipersześcian. f ( X ) można rozłożyć w następujący sposób:

{\ Displaystyle Y = f_ {0} + \ suma _ {i = 1} ^ {d} f_ {i} (X_ {i}) + \ suma _ {i <j} ^ {d} f_ {ij} (X_{i},X_{j})+\cdots +f_{1,2,\kropki,d}(X_{1},X_{2},\kropki,X_{d})}

₀ gdzie f jest stałą, a _fi_jest funkcją X i _, fij funkcją X i _i X j _, itd. Warunkiem tego rozkładu jest to, że

i_ {s}} (X_ {i_ {1}}, X_ {i_ {2}}, \dots ,X_{i_{s}})dX_{k}=0,{\text{dla }}k=i_{1},...,i_{s}}

tj. wszystkie terminy w rozkładzie funkcjonalnym są ortogonalne . Prowadzi to do definicji terminów rozkładu funkcjonalnego w kategoriach warunkowych wartości oczekiwanych,

{\ Displaystyle f_ {0} = E (Y)}

{\ Displaystyle f_ {i} (X_ {i}) = E (Y|X_{i})-f_{0}}

{\ Displaystyle f_ {ij} (X_ {i}, X_ {j}) = E (Y | X_ {i}, X_ {j}) -f_ {0} -f_ {i} -f_ {J}}

_Z czego widać, że fi jest efektem zmiany _samego Xi ( _znanym jako efekt główny Xi ₎ , a fij jest _efektem jednoczesnej zmiany Xi i Xj , dodatkowo do _efektu ich indywidualnych wariacje . Jest to znane jako interakcja drugiego rzędu . Terminy wyższego rzędu mają analogiczne definicje.

Teraz, dalej zakładając, że f ( X ) jest całkowalne do kwadratu , rozkład funkcjonalny można podnieść do kwadratu i scałkować, dając:

{\ Displaystyle \ int f ^ {2} (\ mathbf {X}) d \ mathbf {X} -f_ {0} ^ {2} = \ suma _ {s = 1} ^ {d} \ suma _ {i_ {1}<\dots <i_{s}}^{d}\int f_{i_{1}\dots i_{s}}^{2}dX_{i_{1}}\dots dX_{i_{s} }}

Zauważ, że lewa strona jest równa wariancji Y , a wyrazy prawej strony są wyrazami wariancyjnymi, teraz rozłożonymi względem zbiorów X _i . To ostatecznie prowadzi do rozkładu wyrażenia wariancyjnego,

{\ Displaystyle \ nazwa operatora {Var} (Y) = \ suma _ {i = 1} ^ {d}V_{i}+\suma_{i<j}^{d}V_{ij}+\cdots +V_{12\kropki d}}

Gdzie

{\ Displaystyle V_ {i} = \ operatorname {Var} _ {X_ {i}} \ lewo (E_ {{\ textbf {X} } _ {\ sim i}} (Y\mid X_ {i})\right)}

,

{\ Displaystyle V_ {ij} = \ nazwa operatora {Var} _ {X_ {ij}} \ lewo (E _ {{\ textbf {X}} _ {\ sim ij}} \ lewo (Y \mid X_{i},X_{j}\right)\right)-V_{i}-V_{j}}

i tak dalej. Notacja X _{~ i} wskazuje zbiór wszystkich zmiennych z wyjątkiem X _i . Powyższa dekompozycja wariancji pokazuje, w jaki sposób można rozłożyć wariancję wyniku modelu na warunki przypisywane każdemu wejściu, a także efekty interakcji między nimi. Razem wszystkie warunki sumują się do całkowitej wariancji danych wyjściowych modelu.

Indeksy pierwszego rzędu

Bezpośrednia miara wrażliwości Si oparta na _wariancji , zwana „wskaźnikiem wrażliwości pierwszego rzędu” lub „wskaźnikiem efektu głównego”, jest określona następująco:

{\ Displaystyle S_ {i} = {\ Frac {V_ {i}} {\ nazwa operatora {Var} (Y)}}}

Jest to wkład w wariancję wyjściową efektu głównego X _i , a zatem mierzy efekt zmiany samego X _i , ale uśredniony względem zmian innych parametrów wejściowych. Jest standaryzowany przez całkowitą wariancję, aby zapewnić wkład ułamkowy. Indeksy interakcji wyższego rzędu S _ij , S _ijk itd. można utworzyć dzieląc inne wyrazy w dekompozycji wariancji przez Var( Y ). Zauważ, że ma to implikację, że

{\ Displaystyle \ suma _ {i = 1} ^ {d} S_ {i} + \ suma _ {i <j}^{d}S_{ij}+\cdots +S_{12\dots d}=1}

Indeks efektu całkowitego

Korzystając z podanych powyżej wskaźników Si _, Sij _i wyższych rzędów, można zbudować obraz znaczenia każdej zmiennej w określaniu wariancji wyjściowej . Jednakże, gdy liczba zmiennych jest duża, wymaga to oszacowania indeksów 2 ^d -1, co może być zbyt wymagające obliczeniowo. Z tego powodu stosowana jest miara znana jako „wskaźnik efektu całkowitego” lub „wskaźnik całkowitego rzędu”, S _{Ti .} Mierzy udział w wariancji wyjściowej X _i , w tym wszelkie wariancje spowodowane jego interakcjami, dowolnego rzędu, z innymi zmiennymi wejściowymi. Podaje się jako,

{\ Displaystyle S_ {Ti} = {\ Frac {E_ {{\ textbf {X}} _ {\ sim i}} \ lewo (\ nazwa operatora {Var} _ {X_ {i}} (Y \ mid \ mathbf { X} _{\sim i})\right)}{\operatorname {Var} (Y)}}=1-{\frac {\operatorname {Var} _{{\textbf {X}}_{\sim i }}\left(E_{X_{i}}(Y\mid \mathbf {X} _{\sim i})\right)}{\operatorname {Var} (Y)}}}

Zauważ, że w przeciwieństwie do S _i ,

{\ Displaystyle \ suma _ {i = 1} ^ {d} S_ {Ti} \ geq 1}

ze względu na fakt, że efekt interakcji między np. _Xi i Xj jest liczony zarówno w S _Ti_, jak i S _Tj . W rzeczywistości suma S _Ti będzie równa 1 tylko wtedy, gdy model jest czysto addytywny .

Obliczanie wskaźników

W przypadku funkcji dających się analizować analitycznie powyższe wskaźniki można obliczyć analitycznie, oceniając całek w rozkładzie. Jednak w zdecydowanej większości przypadków są one szacowane – zazwyczaj robi się to metodą Monte Carlo .

Sekwencje próbkowania

Przykład konstrukcji macierzy A _Bⁱ z d =3 i N =4.

Podejście Monte Carlo polega na wygenerowaniu sekwencji losowo rozmieszczonych punktów wewnątrz hipersześcianu jednostkowego (ściśle mówiąc, będą to pseudolosowe ). W praktyce powszechne jest zastępowanie sekwencji losowych sekwencjami o niskiej rozbieżności w celu poprawy wydajności estymatorów. Jest to znane jako metoda quasi-Monte Carlo . Niektóre sekwencje o niskiej rozbieżności powszechnie stosowane w analizie wrażliwości obejmują sekwencję Sobola i projekt łacińskiego hipersześcianu .

Procedura

Aby obliczyć wskaźniki za pomocą (quasi) metody Monte Carlo, stosuje się następujące kroki:

Wygeneruj macierz próbek N × 2 d , tj. każdy wiersz jest punktem próbkowania w hiperprzestrzeni o wymiarach 2 d . Należy tego dokonać w odniesieniu do rozkładów prawdopodobieństwa zmiennych wejściowych.
Użyj pierwszych d kolumn macierzy jako macierzy A , a pozostałych d kolumn jako macierzy B . To skutecznie daje dwie niezależne próbki N punktów w hipersześcianie jednostki d -wymiarowej.
Zbuduj d dalsze macierze N × d A _Bⁱ , dla i = 1,2,...,d, takie, że i - ta kolumna AB ⁱ jest równa i - tej kolumnie B , a _pozostałe kolumny są z A. _
Łącznie macierze A , B i d A B i określają N ( d +2 ) punktów _w^przestrzeni wejściowej ( po jednym dla każdego wiersza). Uruchom model w każdym punkcie projektowym w A , B , i A _Bⁱ , dając w sumie N ( d +2 ) ocen modelu – odpowiadające f( A ), f( B ) i f( A _Bⁱ ) wartości.
Oblicz wskaźniki wrażliwości, korzystając z poniższych estymatorów.

Dokładność estymatorów jest oczywiście zależna od N . Wartość N można wybrać, dodając kolejno punkty i obliczając wskaźniki, aż oszacowane wartości osiągną pewną akceptowalną zbieżność. Z tego powodu, gdy stosuje się sekwencje o małej rozbieżności, korzystne może być użycie tych, które umożliwiają sekwencyjne dodawanie punktów (takich jak sekwencja Sobola), w porównaniu z tymi, które tego nie robią (takie jak sekwencje łacińskich hipersześcianów).

estymatory

Istnieje wiele możliwych estymatorów Monte Carlo dostępnych dla obu wskaźników. Dwa, które są obecnie w powszechnym użyciu, to

{\ Displaystyle \ nazwa operatora {Var} _ {X_ {i}} (E _ {\ mathbf {X} _ {\ sim i}} (Y | X_ {i})) \ około {{\ Frac {1} {N }}\suma _{j=1}^{N}f\left(\mathbf {B} \right)_{j}\left(f\left(\mathbf {A} _{B}^{i} \right)_{j}-f\left(\mathbf {A} \right)_{j}\right)}}

I

{\ Displaystyle E _ {\ mathbf {X} _ {\ sim i}} \ lewo (\ nazwa operatora {Var} _ {X_ {i}} \ lewo (Y \ mid \ mathbf {X} _ {\ sim i} \ prawo)\prawo)\około {{\frac {1}{2N}}\suma _{j=1}^{N}\left(f\left(\mathbf {A} \right)_{j}- f\left(\mathbf {A} _{B}^{i}\right)_{j}\right)^{2}}}

do oszacowania odpowiednio Si _i S _Ti .

Koszt obliczeniowy

Do oszacowania Si _i S _Ti dla wszystkich zmiennych wejściowych wymagane jest N ( d +2) przebiegów modelu . Ponieważ N jest często rzędu setek lub tysięcy przebiegów, koszt obliczeniowy może szybko stać się problemem, gdy model zajmuje znaczną ilość czasu na pojedynczy przebieg. W takich przypadkach dostępnych jest wiele technik zmniejszających koszt obliczeniowy szacowania wskaźników czułości, takich jak emulatory , HDMR i FAST .

Zobacz też