Mając dany zestaw \(N\) wartości \(x_1 \ldots x_N\) z dowolnego rozkładu, możemy zastosować metodę bootstrap w celu znalezienia przedziału ufności dla wariancji. Jako miarę pewności oszacowania, ustalamy poziom istotności \(\alpha\) jako maksymalne dopuszczane przez nas prawdopodobieństwo, że nasz przedział nie obejmie prawdziwej wartości wariancji. Zazwyczaj przyjmujemy w tym celu wartość \(\alpha = 0.05\) lub mniejszą.
W tym celu wielokrotnie powtarzamy następujące kroki:
import numpy
import scipy.stats
# poziom istotności
alpha = 0.01
# dane wejściowe
N = 5
X = numpy.array([1.2, 2.3, 3.1, 1.7, 2.5])
ile_powtorzen = 10000
# --- WERSJA 1 ---
wariancje = numpy.zeros(ile_powtorzen)
for k in xrange(ile_powtorzen) :
Y = numpy.zeros(N)
for j in xrange(N) :
Y[j] = X[numpy.random.randint(N)]
wariancje[k] = numpy.var(Y, ddof=1)
# --- WERSJA 2 ---
wariancje = numpy.zeros(ile_powtorzen)
for k in xrange(ile_powtorzen) :
Y = X[numpy.random.randint(N, size=N)]
wariancje[k] = numpy.var(Y, ddof=1)
# --- WERSJA 3 ---
wariancje = [ numpy.var(X[numpy.random.randint(N, size=N)], ddof=1) for k in xrange(ile_powtorzen) ]Wynikiem powyższego etapu powinna być lista estymowanych wariancji z poszczególnych „syntetycznych” prób (zestawów danych). Z tej listy możemy obliczyć kwantyle rzędu \(\frac{\alpha}{2}\) i \((1-\frac{\alpha}{2})\); są one już krańcami szukanego przedziału ufności.
# znajdujemy przedział ufności,
# pamiętajmy że funkcja scoreatpercentile przyjmuje wartości w procentach
xL = scipy.stats.scoreatpercentile(wariancje, 100.0 * alpha/2)
xP = scipy.stats.scoreatpercentile(wariancje, 100.0 * (1-alpha/2))
print [xL, xP]