Wnioskowanie Statystyczne

Jak wyznaczyć przedział ufności dla proporcji?

Przy pomocy symulacji!

Mając daną pojedynczą proporcję \(p\) wynikającą z serii danych, np.

W badaniach nad cholesterolem u ludzi stwierdzono, że w grupie 135 badanych z wysokim poziomem cholesterolu 10 osób przeszło zawał serca.
możemy zastosować metodę bootstrap w celu znalezienia przedziału ufności dla tej proporcji. Jako miarę pewności oszacowania, ustalamy poziom istotności \(\alpha\) jako maksymalne prawdopodobieństwo, że nasz przedział nie obejmie prawdziwej wartości. Zazwyczaj przyjmujemy w tym celu wartość \(\alpha = 0.05\) lub mniejszą.

Jedną z metod jest wygenerowanie serii danych \(X\) jako wektora składającego się z zer i jedynek; wprzywołanym zadaniu, należałoby za \(X\) przyjąć wektor składający się z 135 elementów: 10 jedynek i 125 zer.

# proszę zwrócić uwagę na podwójne nawiasy,
# funkcja concatenate przyjmuje pojedynczy argument będący listą tablic
X = numpy.concatenate(( numpy.ones(10), numpy.zeros(125) ))

Na podstawie otrzymanego wektora można wtedy wyznaczyć przedział ufności dla wartości oczekiwanej przy użyciu metody bootstrap. Można bowiem zauważyć, że średnia z takiego wektora ma właśnie postać poszukiwanej proporcji.

Alternatywnie, można w ogóle nie generować tablicy \(X\), zaś „syntetyczne” próby \(Y\) losować z odpowiedniego rozkładu: w tym przypadku byłby to rozkład zero-jedynkowy (Bernoulliego) z prawdopodobieństwem \(p = 10/135\).

W tym celu wielokrotnie powtarzamy następujące kroki:

import numpy
import scipy.stats

# poziom istotności
alpha = 0.01

# dane wejściowe
N = 135
p = 10.0 / 135.0
ile_powtorzen = 10000

srednie = numpy.zeros(ile_powtorzen)
for k in xrange(ile_powtorzen) :
  Y = scipy.stats.bernoulli(p).rvs(N)
  srednie[k] = numpy.mean(Y)

Wynikiem powyższego etapu powinna być lista średnich 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(srednie, 100.0 * alpha/2)
xP = scipy.stats.scoreatpercentile(srednie, 100.0 * (1-alpha/2))
print [xL, xP]
autor: Piotr Różański, ostatnia modyfikacja: 10.04.2016