Centralne twierdzenie graniczne i MATLAB

Nie tak dawno odbyłem bardzo ciekawą rozmowę dotyczącą polityki. Rozmowy w tym temacie nie są zazwyczaj pasjonujące, ale nie tym razem. Mój rozmówca zastanawiał się, jak to jest możliwe, że tak ważną rzecz, jak wyniki wyborów parlamentarnych sonduje się na podstawie badania grupy jedynie 1000 osób. Czy to w ogóle ma to sens? Czy czasem nie ma tu jakiegoś przekrętu? Jak to w ogóle jest możliwe, że sondowana jest tak niewielka grupa osób, a potem podawane są przybliżone wyniki z pewnym, względnie niewielkim, błędem statystycznym? Jak ten błąd jest liczony i dlaczego jest taki, a nie inny?
Choćby nie wiem ile wysiłku włożyć w energiczną gestykulację, naprawdę ciężko jest w rozmowie odpowiedzieć w jasny sposób na postawione pytania i rozwiać wątpliwości. O wiele prościej przedstawić wszystko wykorzystując sprawdzone narzędzie, jakim jest MATLAB.

Zacząć należy od podstaw, czyli od centralnego twierdzenia granicznego. Pomijam stronę matematyczną przechodząc do praktycznej. O co chodzi w tym twierdzeniu? Otóż, jeśli jakieś zdarzenie losowe jest wynikiem kumulacji wielu zdarzeń losowych, to rozkład prawdopodobieństwa dla tego zdarzenia będzie rozkładem normalnym.  Można to sprawdzić empirycznie wykorzystując choćby zwykłą kostkę do gry. O ile ma się kilka dni wolnego czasu na ciągłe powtarzanie rzutu. Jeśli jednak mamy pod ręka MATLABa, to przeprowadzenie symulacji kilku tysięcy rzutów nie powinno stanowić problemu.

x = randi(6) % symulacja rzutu sześciościenną kością do gry
x = randi(6, 1, 10000); % symulacja 10k rzutów sześcienną kością do gry
histogram(x);
xlabel('wynik'); ylabel('liczba zliczeń')
title('Kość sześcienna, 10k powtórzeń')

rzut 1x10k

To, co widać na histogramie nie zaskakuje. Prawdopodobieństwo uzyskania każdego wyniku jest równe - mamy do czynienia z rozkładem jednostajnym.
Zobaczmy jednak co się stanie, jeśli będziemy rzucać dwoma kośćmi do gry (potraktujmy to jako dwa zdarzenia), a interesować będzie nas suma uzyskanych oczek.

x = randi(6, 2, 10000); % symulacja 10k rzutów dwoma sześciennymi kośćmi do gry
x = sum(x); % sumowanie wyników pojedynczego rzutu
histogram(x);
xlabel('wynik'); ylabel('liczba zliczeń')
title('Suma rzutu dwoma kośćmi sześciennymi, 10k powtórzeń')

rzut 2x10

Tutaj wynik też nie jest niespodzianką, nawet dla kilkuletniego dziecka ("Atakuje Cię mroczny troll,  rzuć dwoma kośćmi: 2-10 tracisz turę, 11-12 uciekasz 3 pola do przodu". Strata tury mocno prawdopodobna). Jednak jeśli wykonamy rzut dziesięcioma kośćmi, to oszacowanie prawdopodobieństwa uzyskania określonego wyniku nie będzie już tak intuicyjne. Zresztą zobaczmy. Funkcja histogram może zostać wywołana z argumentem normalization, dzięki czemu zamiast liczby zliczeń możemy obserwować prawdopodobieństwo uzyskania określonego wyniku (liczone empirycznie oczywiście).

x = randi(6, 10, 10000); % symulacja 10k rzutów dziesięcioma sześciennymi kośćmi do gry
x = sum(x); % sumowanie wyników pojedynczego rzutu
histogram(x, 'Normalization', 'probability');
xlabel('wynik pojedynczego rzutu'); ylabel('prawdopodobieństwo uzyskania wyniku');
title('Rzut 10-cioma kośćmi sześciennymi, 10k powtórzeń')

Do histogramu można dodać wykres funkcji gęstości prawdopodobieństwa wyznaczonej na podstawie statystyk dla zmiennej x.

% wyznaczanie statystyk
m = mean(x);
s = std(x);
y = 10:0.1:60; % najmniejszy możliwy wynik to 10, największy to 60
f = exp(-(y-m).^2./(2*s^2))./(s*sqrt(2*pi)); % funkcja gęstości prawdopodobieństwa rozkładu normlanego
hold on
plot(y,f,'LineWidth',1.5)
hold off

rzut 10x10

Proste rzuty kością doprowadziły do uzyskania rozkładu normlanego, opisanego stosunkowo zawiłym wzorem:

P(x) = \frac{1}{{\sigma \sqrt {2\pi } }}e^{- ( {x - \mu } )^2 / {2\sigma}^2 }

Nie musimy wyznaczać sumy wyników poszczególnych zdarzeń losowych, by obserwować takie "czary". Wróćmy do początkowego rzutu pojedynczą kością. Wiemy, że prawdopodobieństwo uzyskania w wyniku szóstki wynosi 1/6. Spróbujmy to jednak sprawdzić empirycznie na tej zasadzie, że wykonamy 600 rzutów i sprawdzimy, ile razy wypadła szóstka.

n = 600;
x = randi(6,n,1);
x6 = (x == 6); % interesuje nas, czy wypadła 6 czy nie
s6 = sum(x6)
p6 = mean(x6)

s6 = 97
p6 = 0.1617

No cóż, w pojedynczym eksperymencie udało mi się uzyskać powyższy wynik. Jeśli wykonam symulację jeszcze raz, to uzyskam najpewniej nieco inną wartość wyjściową. Na pewno wynik będzie oscylował wokół wartości s6 = 100 i p6 = 0.1667,  pytanie tylko w jaki sposób. Można to sprawdzić.

n = 600; % liczność próby, na podstawie której szacuje prawdopodobieństwo uzyskania 6
N = 10000; % ilość przeprowadzonych eksperymentów
x = randi(6,n,N); % macierz losowa
x6 = (x == 6); % czy w danym rzucie wypadło 6? 
s6 = sum(x6);
p6 = mean(x6) ;
bins = (65:135 ); % przedziały
subplot(2,1,1);
histogram(s6, bins);
xlabel('szóstki w próbie');
ylabel('ilość prób');
subplot(2,1,2);
histogram(p6, bins/n, 'Normalization', 'probability');
xlabel('prawdopodobieństwo wypadnięcia 6') ; 
ylabel('pr. "trafienia na próbkę"');

matrix

Znów rozkład normalny! To jest ten moment, w którym bez pomocnika takiego jak MATLAB wszelkie tłumaczenia nie mają sensu. To trochę tak jak z filmem Incepcja - to po prostu trzeba zobaczyć. Tutaj też mamy "sen w śnie". Pojedynczy eksperyment, w którym szacujemy prawdopodobieństwo, daje nam jakiś wynik. Wynik ten jest jednak niepewny. Jeśli przeprowadzimy ileś kolejnych eksperymentów, to ich wyniki ułożą się zgodnie z rozkładem normalnym, a "prawdziwego" wyniku będziemy mogli doszukiwać się gdzieś w środku. Mamy tu więc pewne prawdopodobieństwa oszacowania określonego prawdopodobieństwa.
Jeśli chcemy wiedzieć jakie jest prawdopodobieństwo wyrzucenia 6 w pojedynczym rzucie kością, to najskuteczniej jest wykonać jeden porządny eksperyment i przeprowadzić symulację 6 milionów rzutów. Czasem jednak nie można sobie na to pozwolić.

Możemy zatem rzucić kością 600 razy i powiedzieć: tym razem wyszło 97/600, ale trzeba uwzględnić błąd. Zakładamy optymistycznie, że nasz wynik to najlepszy możliwy, jednak jeśli nasz kolega też rzuci 600 razy i wyjdzie mu co innego, to też może być w porządku. Nasza ocena poprawności wyniku eksperymentalnego kolegi będzie zależeć od tego, o ile jego wynik będzie oddalony od naszego. Pytanie tylko, o ile może być oddalony od naszego, żeby uznać go za w miarę prawdopodobny. Możemy w tym momencie posłużyć się rozkładem normlanym i na jego podstawie orzec. Jako że nasz wynik jest najlepszy (trzeba mieć wiarę w własne wyniki), to przyjmujemy, że wartość średnia, czyli punkt centralny "dzwonka" wynosi 97. Do wyznaczenia funkcji potrzebna będzie jeszcze wartość sigma. Tę szacujemy na podstawie wzoru:

{ {{\sigma}^2}_{se}} = {\sigma}^2 / n

 \sigma po prawej stronie to odchylenie standardowe dla naszej próbki. Sigma po lewej - błąd standardowy.

Ostatecznie wszystko to wygląda tak:

 

n = 600; % liczność próby, na podstawie której szacuje prawdopodobieństwo uzyskania 6
x = randi(6,n,1); % rzut kością (600 x)
x6 = x == 6; % czy w danym rzucie wypadło 6?
p6 = mean(x6); % wartość średnia (czyli prawdopodobieństwo wypadnięcia 6 dla próbki)  
s = std(x6); % odchylenie standardowe dla próbki

se = s / sqrt(n); % błąd standardowy
% zdefiniowanie funkcji gęstości prawdopodobieństwa (funkcja anonimowa)
f = @(x)(exp(-(x-p6).^2./(2*se^2))./(se*sqrt(2*pi))); 

x = 0.1:0.001:0.3;
x2s = (p6-2*se):0.001:(p6+2*se); 
plot(x, f(x), x2s, f(x2s));
hold on;
stem(1/6, max(f(x)));
xlabel('prawdopodobieństwo');
title('funkcja gęstości prawdopodobieństwa');

wykres f. gestosci pr.

Zaznaczony powyżej obszar określa wartości wyznaczonego prawdopodobieństwa, które pojawią się na ~95%, jeśli powtórzymy eksperyment (odległość 2\sigma od średniej) przy założeniu, że nasza wartość średnia, to ta prawidłowa. Jeśli nasz kolega przyjdzie ze swoimi wynikami i będą się one mieścić w tym zakresie, to przyjmujemy, że faktycznie przeprowadził eksperyment, a nie wyciągnął liczbę z kapelusza. To znaczy - jest szansa, że na 600 powtórzeń ani razu nie wypadła mu szóstka, ale bez przesady. Za potencjalnie prawdziwe przyjmujemy tylko te wyniki w obszarze. Do czasu. W końcu ktoś rzuca kością 6000000 razy i wychodzi na to, że nasz wynik nie jest jednak taki super (pionowa linia na wykresie to "prawdziwe" prawdopodobieństwo). Nie szkodzi, jesteśmy kryci. Jeśli prawdziwy wynik jest w odległości  2\sigma od naszego,  to znaczy,  że my jesteśmy w odległości 2\sigmaod prawdziwego wyniku, czyli w dopuszczalnym zakresie.

Przenosząc to na szacowanie wyników wyborów, przyglądnijmy się partii PSL, która w ostatnich wyborach parlamentarnych balansowała w na krawędzi progu wyborczego. Ostatecznie partia zdobyła 5,13 procent głosów i do sejmu się dostała. Czy badanie exit poll prowadzone na grupie 90000 wyborców ma szansę prawidłowo określić czy partia dostanie się do sejmu?

n = 90000;
p = 5.13 / 100;
% wyznaczenie wariancji (rozkład Bernouliego)
v = p*(1-p); 
se = sqrt(v / n); % błąd standardowy
f = @(x)(exp(-(x-p).^2./(2*se^2))./(se*sqrt(2*pi))); 
x = 0.045:0.0001:0.055;
plot(x, f(x));
hold on;
stem(0.05, max(f(x)));
stem(p - 1.96*se, f(p-1.96*se));
xlabel('szacowany wynik')
legend({'f. gęstości prawd.', 'próg wyborczy', 'odległość 2\sigma od \mu'}, 'Location', 'best')

exit poll psl

Udało się trafić na skrajny przypadek. Nawet gdyby założyć, że badanie było przeprowadzone idealnie, to jest względnie duża szansa, że nie udzieliłoby prawidłowej odpowiedzi na pytanie, czy PSL dostanie się do sejmu. Tzn. gdybyśmy przeprowadzili 100 takich badań, to w kilku wypadkach wynik byłby poniżej progu 5%. Choć z drugiej strony patrząc, te 100 badań to też próba losowa, więc jest pewna mała szansa, że wszystkie dałyby wynik poniżej progu...

W związku z powyższym, na zakończenie polecam kolejny ciekawy film: Matrix. Poza tym zagadnienia związane ze statystyką dobrze tłumaczy Sal Khan z KhanAcademy. Tłumaczy bez MATLABa i robi to dobrze.

(Visited 682 times, 1 visits today)

Dodaj komentarz

Twój adres email nie zostanie opublikowany. Pola, których wypełnienie jest wymagane, są oznaczone symbolem *