Rachunek prawdopodobieństwa i statystyka Home
I rok geoinformatyki uruchomione 2020 03 12, ostatnia modyfikacja 2020 03 28

Na czas zamkniętej Uczelni. Pytania kierujemy na a@adamwalanus.pl; temat maila: zajęcia.

Rachunek prawdopodobieństwa

Chaos deterministyczny     Rozkład dwumianowy     Inne rozkłady prawdopodobieństwa     Gęstość prawdopodobieństwa     Histogram     Rozkład normalny     Skala i przesunięcie     Skumulowany rozkład     Rozkład wykładniczy     Rozkład Poissona     Parametry rozkładów     Zależność    

Statystyka

Dlaczego statystyka, to nie matematyka?     Średnia     Odchylenie standardowe     Prawda statystyczna     Test t     Test statystyczny     Test χ2     Regresja     Szeregi czasowe i ciąg dalszy, w R

Fake news a statystyka     Zaraza


Chaos deterministyczny

1. Można (nie zagłębiając się) googlować hasła "chaos deterministyczny" i "efekt motyla".

2. Googlujemy hasło: 4x(1-x). Analizujemy (=przyglądamy się) otrzymaną parabolę.

3. Ewentualnie przeglądamy poniższy obrazek

który można ściągnąć w postaci arkusza Excela.

4. Bierzemy jakiś arkusz kalkulacyjny, np. Arkusz Google, MS Excel, OpenOffice Calc. Zapoznajemy się z nim.

5a. W komórce A1 (lub którejś innej) wpisujemy liczbę z otwartego przedziału: (0, 1), ale nie ćwierci: 0.25, 0.5, 0.75. Uwaga, liczby na tej stronie będą pisane z kropką jako separatorem dziesiętnym, w swoim arkuszu robimy jak chcemy.

5b. W komórce niżej (A2) wpisujemy formułę
=4*A1*(1-A1)
i przeciągamy ją w dół, tak żeby mieć w sumie z 50 komórek.

5c. Tworzymy (Wstaw, Insert) wykres z tych pięćdziesięciu komórek; liniowy, żeby wyglądał, jak ten na powyższym rysunku.

5d. Zmieniamy dowolnie pierwszą liczbę (A1), obserwujemy wykres.

5e. Ewentualnie podczytujemy Odwzorowanie logistyczne lub Logistic map.

6a. Do A1 wpisujemy jakąś prostą liczbę, np. 0.2.

6b. Kopiujemy nasze 50 komórek do sąsiedniej kolumny (B). Zmieniamy liczbę w B1, np. na 0.21.

6c. Wyrzucamy wykres i tworzymy nowy (można nie wyrzucać, tylko modyfikować, ale wyrzucić łatwiej). Ma to być wykres liniowy, z obu kolumn. Powinniśmy otrzymać na nim nałożone na siebie dwie podobne linie wykresu.

7a. Obserwujemy od którego punktu wykresy "rozjeżdżają się", tzn. wartości w kolumnach A i B zaczynają się poważnie różnić, np. o 0.3 (może tu być trochę nie matematycznej subiektywności, w końcu mamy do czynienia z chaosem). W przypadku A1=0.2, B1=0.21 będzie to w szóstym wierszu.

7b. Zapisujemy w arkuszu, np. w kolumnach C i D różnicę B1-A1 i numer wiersza. Tak by to mogło wyglądać: w C1: 0.01 a w D1: 6. Uwaga! Jak rozumieć słowo "zapisujemy"? Możemy zapisac formułę albo wartość. Ponieważ za chwilę będziemy zmieniać B1, musimy zapisywać (kolejno w C) wartości. Możemy nawet ręcznie je przepisywać.

7c. Zmieniamy liczbę w B1 na 0.201. Wykonujemy 7b. Dodajemy dalsze zera (0.2001), do oporu (jakiego rodzaju opór napotkamy?).

8a. Tworzymy wykres rozrzutu (x,y) z kolumn C i D.

8b. W kolumnie E liczymy logarytm dziesiętny z C, do F przepisujemy D (żeby łatwiej się robiło wykres).

8c. Tworzymy wykres rozrzutu (x,y) z kolumn E i F. Dopasowujemy do wykresu prostą y=ax+b. Zapamiętujemy a - to nasz Wynik.

9. Uruchamiamy ten programik. Można podejrzeć jego kod, ale poza technicznymi sprawami jest tam tylko jedna istotna linijka: y=R*y*(1-y). Otóż tu czynnik R można zmieniać. Dotąd mieliśmy R=4. Jasne, że w arkuszu również możemy sobie R zmienić. Wartości R>4 i R<0 są bez sensu, ciekawie jest od 3, najciekawiej w okolicy 3.8.

9a. Zadanie: znaleźć najciekawsze R.

9b. Tu jest podsumowanie; R = 2.9 do 4. (Warto by dorobić do tego programiku skalę poziomą, czyli oś R.)


Rozkład dwumianowy ↑↑↑

1. Rzucamy monetą. Spadnie na podłogę orłem lub reszką do góry. Oba wyniki są jednakowo prawdopodobne, bo ta drobna różnica kształtu powierzchni nie może przecież mieć wpływu (oczywiście warto by to SPRAWDZIĆ, ale nie wiemy ile razy trzeba by rzucać). Tak więc przyjmujemy, że prawdopodobieństwo otrzymania orła (reszki) p=0.5.

2. Rzucamy tą monetą N razy. Lub rzucamy raz N monetami, wszystko jedno. To N mogłoby wynosić 2, 10, 100. Otrzymujemy k orłów i N-k reszek. Ile wynosi k? Może być k=0, 1, 2, … N. A ile wyjdzie? Jakie jest prawdopodobieństwo, że przy N=10 otrzymamy k=5? Pewnie największe. Pewne jest jedno, że suma wszystkich prawdopodobieństw będzie równa 1; p(0)+p(1)+ … +p(N)=1, bo jakieś k otrzymamy na pewno.

3. Zróbmy model w arkuszu kalkulacyjnym, nie będziemy przecież N razy rzucać monetą. Moglibyśmy użyć naszych liczb 4x(1-x). Na przykład, można by zrobić tak: orzeł jeśli kolejna liczba jest mniejsza od 0.5, reszka w przeciwnym wypadku. Jednak lepiej będzie użyć profesjonalnego, sprawdzonego generatora liczb losowych. W Excelu nazywa on się LOS(), w JavaScript jest to Math.random(), a w arkuszu Googla rand(). Jak to działa? Podobnie jak 4x(1-x). To są deterministyczne generatory losowości, czyli generują pseudolosowość. Nam to wystarczy. W bankach potrzebują prawdziwej losowości i używają drogich, zewnętrznych, fizycznych urządzeń.

Ciekawa rzecz, te funkcje (jak to funkcje) mają nawiasy, żeby umieścić w nich argument, a przecież żadnego argumentu nie potrzeba.

3a. Wpisujemy więc do którejś komórki =rand() i przeciągamy w dół N=10 razy. Patrzymy na te liczby. Mając wybraną którąś pustą komórkę naciskamy klawisz Delete, żeby otrzymać nowe liczby (w Excelu F9).

0.5375650976 0.3254018231 0.5787474213 0.1045142063 0.1116005375 0.1271713004 0.8296205742 0.6400735714 0.6014123909 0.2384220984

3b. Modyfikujemy formułę: =rand()>0.5 i przeciągamy. Wyniki są: Prawda lub Fałsz. Nie musimy mieć: "Orzeł", "Reszka".

3c. Policzmy Prawdy. Wpisujemy w którejś komórce (może w A11) =COUNTIF(A1:A10,true) i otrzymujemy liczbę. Może nam wyjść od 0 do 10, ale wyjdzie raczej 4, 5, 6, no może 3 lub 7.

3d. A może jednak poczekamy na wynik 0 lub 10. Wystarczy trzymać naciśnięty klawisz Delete. Gorzej z zauważeniem tego 0. Jakie jest prawdopodobieństwo otrzymania 10 prawd?

Prawdopodobieństwo otrzymania w pierwszej komórce Prawda wynosi 0.5 (co łatwo zobaczyć, mniej więcej). Żeby również w drugiej komórce była prawda, musi drugi raz zajść zdarzenie, którego prawdopodobieństwo jest 0.5. Czyli w każdej z dwóch niezależnych prób musimy otrzymać konkretny wynik.

Prawdopodobieństwo takiego zdarzenia jest iloczynem prawdopodobieńśtw tych dwóch zdarzeń
p = 0.5 * 0.5 = 0.25
to ważne. Jeżeli prawdopodobieństwo uszkodzenia dysku z naszymi bezcennymi danymi wynosi p=1e-6, to prawdopodobieństwo uszkodzenia obu dysków, również tego z BackUpem, jest kwadratem tego prawdopodobieństwa: p=1e-12. Pamiętać musimy tylko, żeby te uszkodzenia były niezależne od siebie, czyli dyski nie mogą leżeć w tej samej szufladzie, którą zaleje nam woda z kwiatków.

Prawda w trzech kolejnych komórkach? p = 0.5 * 0.5 * 0.5 = 0.5^3 = 0.125.
We wszystkich dziesięciu? p = 0.5^10 = 1/1024. To jak długo trzeba trzymać wciśnięty Delete? Tak, żeby z 1000 razy sie przeliczyło. Jak liczyć ile razy nam się przeliczyło? I czy wyszło 10 prawd? Googlujemy hasło "Iterative calculation".

3e. Wracamy do pytania o prawdopodobieństwo otrzymania k=5 prawd. Zacznijmy jednak od k=1. Zauważmy, że przypadek k=0 (same Fałsze) realizuje sie na jeden sposób
FFFFFFFFFF. Natomiast w przypadku k=1 mamy 10 różnych układów

PFFFFFFFFF
FPFFFFFFFF
FFPFFFFFFF
...
FFFFFFFFFP.

Każdy z powyższych przypadków jest tak samo prawdopodobny i p=1/1024. Po prostu, żeby otrzymać FPFFFFFFFF 10 razy musi wypaść co trzeba.
Każdy z tych dziesięciu przypadków jest unkalny, jak ON wypadł, to nie wypadł inny.
Tak więc mamy tu sumę prawdopodobieństw. Prwadopodobieńtwo otrzymania któregokolwiek z powyższych przypadków (zwierających jeden fałsz) jest sumą:

p = 1/1024 + 1/1024 + 1/1024 + 1/1024 + 1/1024 + 1/1024 + 1/1024 + 1/1024 + 1/1024 + 1/1024 = 10/1024 = 0.00977

p(k=1) = N * 0.5^N

Tak samo dla k=9. Ale jak z k=2, 3, itd. To już jest KOMBINATORYKA.

4. Kombinatoryka. Na ile sposobów można zapisać 10 liter alfabetu?
ABCDEFGHIJ - 1
BCDEFGHIJA - 2 itd. Zasada jest taka: A umieszczamy na jednej z 10 dostępnych pozycji. B na jednej z pozostałych 9, itd. Tak więc sposobów tych jest 10*9*8*7*6*5*4*3*2*1 = 10! = 3628800 ( silnia = factorial = FACT() ).

4a. Tylko, że my mamy tylko dwie litery: P i F. Mamy mieć k Prawd. Każde P jest takie samo jak inne P (a przynajmniej my ich nie rozróżniamy, mamy k tych P). Liczbę k RÓŻNYCH P można rozmieścić na k! sposobów. Natomiast N-k RÓŻNYCH F można rozmieścić na (N-k)! sposobów. Więc to 10! = N! trzeba dwa razy podzielić.

N! / ( k! * (N - k)! )      - to się nazywa symbol Newtona (patrz też: trójkąt Pascala)

5. Liczbę RÓŻNYCH sposobów mnożymy przez prawdopodobieństwo p=0.5^N i mamy odpowiedź na pytanie o prawdopodobieństwo otrzymania k orłów w N rzutach monetą:

p(k)=N!/(k!*(N-k)!)*0.5^N     

Dla N=10 i k (np. k=5) umieszczonym w komórce, powiedzmy C1, formuła w arkuszu Googla będzie taka
=0.5^10*FACT(10)/(FACT(C1)*FACT(10-C1))
czy wyszło 0.246? Dla k=5.

5a. Wpiszmy do C1 wartość 0, a dalej (niżej) 1, 2, ... 10 - wszystkie możliwe k. Przeciągamy formułę. Mamy ROZKŁAD PRAWDOPODOBIEŃSTWA. Warto zrobić wykres z zakresu C1:D11. Arkusz Googla sam się domyśla, jaki to ma być wykres.

5b. SPRAWDŹMY czy =SUM(D1:D11) = 1.

5c. Tak to wygląda, najbardziej prawdopodobne jest 5, potem spada, oczywiście symetrycznie. Ciekawy kształt - dzwon - krzywa dzwonowa.

6. A gdyby moneta był krzywa i p(orzeł)=0.51. Albo, gdyby chodziło o mecz piłkarski i naszej drużynie przydzielamy p=0.9?
Sprawa jest prosta, zamiast 0.5^N będzie
p^k * (1-p)^(N-k) bo ma być k "sukcesów", a więc musi być też N-k "porażek", których prawdopodobieństwo wynosi 1-p (jeżeli p(sukcesu)=p). Wszystkie N prawdopodobieństw mnożymy, np. 0.9 * 0.9 * 0.9 * 0.9 * 0.1 * 0.1 * 0.1 * 0.1 * 0.1 * 0.1.

6a. Sprawa symbolu Newtona pozostaje bez zmian, tak więc tak wygląda
DWUMIANOWY rozkład prawdopodobieństwa w ogólnym przypadku

p(k) = N! / ( k! * (N-k)! ) * pk * (1-p)N-k

p(k) to prawdopodobieństwo k "sukcesów" w N próbach (meczach), przy prawdopodobieństwie sukcesu równym p.

6b. Zróbmy taki arkusz:

W kolumnie A mamy N i p, w kolumnie B mamy kolejne k, w kolumnie C mamy p(k). Zamiast męczyć się z silniami użyliśmy gotowej formuły na rozkład dwumianowy (binomial distribution). Jaki sens ma ostatni jej argument, który tu wynosi: false? Zostawmy to na później. Teraz bawimy się zmieniając N i p. Warto zacząć od dużego N, zrobić N+1 wierszy i narysować wykres. Potem jak zmniejszymy N, wykres będzie działał. Przy większych N (20, 30) i p niezbyt odległym od pół, otrzymujemy piękny dzwon.

7. Dla zaawansowanych. Formuła =rand() zwraca losową liczbę z przedziału (0, 1). Formuła =rand()>0.5 zwraca wynik rzutu monetą. A jak użyć rand do otrzymywania od razu losowej liczby k (przy danych N i p)?
kostka
Modyfikujemy formułę BINOM z false na true. Mamy kumulację. W C1 wpisujemy zero. w D2:D6 wpisujemy formułę
=IF(AND(A$5>C1,A$5<=C2),B2,"")
W A5 oczywiście siedzi RAND. Naciskając klawisz Delete obserwujemy generowanie losowej liczby sukcesów (jak w D9 wpiszemy =SUM(D2:D8), to wynik będzie w jednej komórce). Uwaga! to nie jest, po prostu losowa liczba z przedziału od 0 do 6, ta liczba ma rozkład dwumianowy; 0 i 6 nam się szybko nie pojawią, w przeciwieństwie do 2 i 3.


Inne rozkłady prawdopodobieństwa ↑↑↑

1. W poprzednim temacie było o rozkładzie dwumianowym prawdopodobieństwa. Dlaczego dwumianowy? Bo ma coś wspólnego z dwumianem (a+b)N. Ale nazwy można sobie zawsze wygooglować. Oczywiście, była tam już mowa o najprostszym rozkładzie, o rozkładzie jedynki na dwie połówki. Że prawdopodobieństwa orła i reszki wynoszą 0.5. Była też mowa, że może być tak, że 1=0.51+0.49, jeśli moneta jest krzywa. Albo, że p=0.1 a wtedy (1-p)=0.9 – bo grają słabsi z silniejszymi.

2. Nowością będzie kostka do gry, bo ma nie dwie, ale sześć ścianek. Jak jest równa, ro rozkład jest taki, że każda z ponumerowanych ścianek (1, 2, … 6) ma takie samo prawdopodobieństwo, i musi to być więc 1/6 = 0.166667. Przyjrzyjmy się liczbie 0.166667, która nie jest dokładna, dokładnie byłoby z kropeczkami: 0.1666…, które co oznaczają? Oznaczają nieskończoność. Z tym trudnym tematem (z rachunkiem różniczkowym i całkowym) będziemy musieli się zmierzyć w następnym temacie, kiedy trzeba będzie coś powiedzieć o gęstości prawdopodobieństwa. A trzeba będzie, bo zdarzają się doświadczenia losowe, które mają (w zasadzie, nie w realu) nieskończenie (nieprzeliczalnie) wiele wyników. Mieliśmy już z tym do czynienie; używaliśmy rand().
tak wygląda rozkład prawdopodobieństwa dla sześciennej kostki.

2a. Jak realizujemy rzut kostką w arkuszu? Tak: =TRUNC(rand()*6+1,0)

2b. Jakie mamy jeszcze kostki do gry? Jest ich 5, nazywają się "bryły Platońskie". Jakie tam są prawdopodobieństwa?

2c. Jakie mamy jeszcze inne rozkłady prawdopodoobeństwa? Np. taki: w gonitwie startuje 10 ponumerowanych koni. W ocenie pewnego fachowca mają takie prawdopodobieństwa zwycięstwa: 1-5%, 2-10%, 3-5%, 4-60%, 5-20%, 6-2%, … i jakoś tak dalej, byleby procenty sumowały się do 100%. Oczywiście prawdopodobieństwo możemy wyrażać w %, bo to tylko inny zapis ułamka.

2d. Jak generujemy w arkuszu liczbę losową z dowolnego rozkładu? Jak ten, zapisany w A i B.

W C robimy kumulację, dodając zero na początku. W D dajemy rand. W E2:E6 skomplikowaną formułę, jak widać. W E1 dajemy =SUM(E2:E6), żeby mieć wynik w konkretnej komórce.

3. A co z rzutem dwiema kostkami. Podobnie jak z rzutem dwiema monetami, ale nie całkiem tak samo. Sumaryczna liczba oczek może być 2, 3, ... 12. Wszystkie wyniki są tak samo prawdopdobne? Nie, bo 2 i 12 otrzymamy, TYLKO gdy na obu kostkach będzie co trzeba. Obie kostki muszą dać 1 (6). Tak więc prawdopodobieństwo otrzymania 2 (12) będzie 1/6 * 1/6 = 1/36 = 0.27778.
Za to sumaryczną liczbę oczek 7 możemy otrzymać aż na 6 sposobów: 1+6, 2+5, 3+4, 4+3, 5+2, 6+1. Przecież 1+6 i 6+1 to to samo. Jednak to są dwie RÓŻNE kostki. Albo jedna kostka rzucona raz i POTEM drugi raz. Tylko cząstki elementarne są identyczne.
Tak więc prawdopodobieństwo p(7) = 6 * 1/36 = 1/6. Podobnie liczymy inne. Rozkład będzie trójkątny.

Jak będzie przy 3 i więcej kostkach? Sumowanie to ważna rzecz. Większość zdarzeń naszego życia jest sumą wielu składników (DNA, wychowanie, spóźniony tramwaj ...). Z kostkami wygląda to tak:

Nazywa się to "splot" ("convolution"). Zostawmy sumowanie na później.


Gęstość prawdopodobieństwa ↑↑↑

1. Pytanie: Ile waży woda? Zrozumiane dosłownie, nie ma odpowiedzi, bo nie wiemy ile wody pytający ma na myśli.
Jaka jest gęstość wody? Odpowiedź: 1 kg/dm3.
Jakie jest prawdopodobieństwo 3 oczek na kostce? Odp.: 1/6.
Jakie jest prawdopodobieństwo, że rand()=0.2? Odp.: zero. Ścisła odpowiedź: ok. 1e-15. Zostańmy przy zerze.
Dlaczego zero? Bo o ile kostka ma 6 ścianek, to różnych liczb w przedziale (0, 1) jest nieskończenie wiele.
Filozofia: W arkuszu, ani nigdzie na ziemi nie ma nieskończoności. Liczby podwójnej precyzji, jakich używa nasz arkusz mają 15-16 cyfr znaczących. Różnych ciągów 15 cyfr (tych za zerem i kropką) jest 1e15. Gdyby to były liczby poczwórnej precyzji ... nie to nas interesuje. Dlatego zadajemy takie pytanie:
Jakie jest prawdopodobieństwo trafienia w przedział (0.1, 0.3). Odp.: 0.2.
To ważne, przetłumaczmy to zdanie na obliczenia: A1 =rand(); A2 =AND(A1>0.1,A1<0.3). Pytanie teraz brzmi: jakie jest prawdopodobieństwo otrzymania prawdy w A2? Dlaczego wynik jest p=0.2, czyli tyle ile długość przedziału (0.1, 0.3)? Bo GĘSTOŚĆ PRAWDOPODOBIEŃSTWA wynosi 1. Czy gęstość prawdopodobieństwa zawsze wynosi 1? Nie. Weźmy taki przykład: =20*rand(), teraz gęstość prawdopodobieństwa wynosi 0.05. Oczywiście tylko w przedziale (0, 20), poza nim wynosi 0. Zauważmy, że 0.05*20=1.

2. Weźmy inny przykład =rand()^2 . Pogenerujmy trochę. Jakoś częste są małe liczby, a przekroczenie 0.5 na pewno jest rzadsze niż połowa przypadków. No tak, prawdopodobieństwo trafienia rand'em w (0, 0.1) albo (0.9, 1) jest 0.1. Ale jak rand podniesiemy do kwadratu, to te przedziały przekształcą się w (0, 0.01) i (0.81, 1). Ten pierwszy jest krótszy od oryginału 10 razy, a ten drugi dłuższy ok. 2 razy. To oznacza, że gęstość prawdopodobieństwa blisko zera jest 10 = 0.1/0.01, a blisko 1 jest 0.5 = 0.1/0.2. Czujemy, że trzeba przejść na krótsze przedziały.
Poniżej, w A są granice przedziałów, w B jest PRAWDOPODOBIEŃSTWO dla danego przedziału (gęstość byłaby 0.01/0.01=1).

W kolumnie C są granice przedziałów podniesione do kwadratu. W D jest prawdopodobieństwo z B podzielone przez aktualną długość przedziału (C2-C1 itd.) - to jest GĘSTOŚĆ PRAWDOPODOBIEŃTWA.
Podzieliliśmy odcinek (0, 1) na sto mniejszych. Na ile powinniśmy podzielić? Na nieskończenie wiele. Tak się to robi w matematyce, konkretnie w rachunku różniczkowym. W D powinna być funkcja. Jaka? Najwyraźniej f(x → 0) → ∞. f(1)=0.5. Może jest to po prostu odwrotność pochodnej z x^2, czyli f(x)=1/(2x)? Dodajmy do powyższego arkusza kolejną kolumnę i wpiszmy do E2 =0.5/((A1+A2)/2). Zgadza się. Dlaczego średnia z sąsiednich granic przedziału? To drobiazg. Dlaczego odwrotność pochodnej?

3. Podsumujmy. Co to jest funkcja gęstości prawdopodobieństwa f(x)? Otóż f(x) to taka liczba (dla danego x), że prawdopodobieństwo trafienia w przedział (x, x+dx) wynosi dp=f(x)*dx. Uwaga, dx jest maaaałe. Jeżeli jest nieskończenie małe, to nie piszemy gwiazdki i mamy różniczki dp=f(x)dx. Wtedy dp też jest nieskończenie małe, a f(x)=dp/dx (mówi nam to coś?).

4. Tak wygląda najważniejsza (najważniejsza to ona nie jest, ale jest najprostsza) funkcja gęstości prawdopodobieństwa (nasze rand):

Czy występuje gdzieś w praktyce, a nie tylko w obliczeniach? Tak, w mechanice, jak są luzy, albo gdy ktoś chce się precyzyjnie wyrazić o czasie, a mógł tylko prze moment rzucić okiem na zegar, na którym zobaczył 12:05. Co może powiedzieć, wobec ewidentnej niepewności? Że t, wyrażone w minutach po godzinie 12, ma rozkład równomierny w przedziale (5, 6). Czyli jednakowo prawdopodobne jest 5.1 minuty, jak i 5.9. Tu uwaga, zmienna na osi x ma wtedy wymiar, jednostką jest minuta. Jaki wymiar będzie miała f(t) na osi pionowej. To dziwne, ale będzie to 1/min. Tak, żeby prawdopodobieństwo, że czas był w przedziale (5.0, 5.3) wynosiło 0.3, a nie 0.3 minuty.

5. Trochę wyżej jest wykres z daszkiem albo trójkątem. Rzucaliśmy dwiema kostakami. Teraz zróbmy taka sumę: =rand()+rand(). Też będzie trójkąt? Popatrzmy w arkuszu. Może wyjść od 0 do 2, najbardziej prawdopodobne będzie 1? Jak to porządniej sprawdzić? Da się w arkuszu, ale może by użyć czegoś mądrzejszego?
Rozwińmy Tools i wybierzmy Script editor. Zdefiniujmy własną funkcję:
function troj(x,dx,n) {var j=0
 for (var i=0;i<n;i++){var y=Math.random()+Math.random()
  if (y>x && y<x+dx){j++} }
 return j/n }
Trzeba ją zapisać, a potem można wpisać do arkusza =troj(A1,0.1,100000). W A1 podajemy x, np. 0, 1, 1.9. Funkcja zwraca, jaka część generowanych liczb sum mieści się w przedziale (x, x+0.1). Dla x=1 wynik będzie ok. 0.1, dla x=0 ok. 0.005 - to jest pole małego, prostokątnego trójkącika równobocznego o bokach 0.l.


Histogram ↑↑↑

1. W poprzednim temacie użyliśmy zaawansowanej techniki (skryptu), żeby zobaczyć jak wygląda funkcja gęstości prawdopodobieństwa. Konkretnie mogliśmy sprawdzić, jaką, mniej więcej ma wartość dla dowolnego x. Uwaga, wyżej liczyliśmy prawdopodobieństwo, żeby to była gęstość, trzeba by jeszcze podzielić wynik przez dx.

2. Histogram - dziwna nazwa, ale unikalna, Googluje się jak trzeba. To te słupki. Zróbmy tak. Włączmy ten programik. Mamy jakiś histogram. W arkuszu wygenerujmy z 10000 liczb =rand()+rand() i wklejmy je do pola tekstowego, na miejsce tamtych danych (Ctrl+A, Ctrl+V). Jest trójkąt? Można wrócić do pojedynczego =rand(), a potem wziąć =rand()+rand()+rand()+rand(), żeby powoli przechodzić do następnego tematu.

3. Histogram zlicza ile mamy wyników w różnych przedziałach. Te przedziały zwykle są tej samej długości. Nie może ich być za dużo, bo będą zawierały mało liczb i słupki będą nierówne. Wypróbujmy to.

4. Zastosowanie histogramu w fotografii.
Dlaczego poważniejsi fotografowie patrzą na histogram? Bo patrząc na zdjęcie na ekranie można się pomylić, co do jego jasności czy ciemności. Co prawda możemy mieć na zdjęciu sygnalizację obszarów, które wyszły poza zakres (były ciemniejsze niż 0, a z konieczności zapisały się jako 0, albo były jaśniejsze niż 255, a zapisane zostały z tą maksymalną wartością). Ale to mało. Chcemy (my, poważni fotografowie) obiektywnie widzieć, czy mamy więcej jasnych, czy ciemnych obszarów.
Czy na pewno potrzebujemy osobno trzech histogramów RGB? Taki przykład: czerwone maki. Ich kolor dokładnie trafia w R (matrycy i naszego oka). Jeżeli płatki zajmują tę część obrazu, gdzie akurat mamy ustawiony pomiar światła, to ekspozycja będzie poprawna, ale całe światło pójdzie do R. Potem będziemy się dziwili płaskim, nieciekawym obszarom z R=255. A na histogramie widzimy to od razu. Dajemy korektę ekspozycji (na minus) i powtarzamy zdjęcie.
W punkcie 3 powyżej, jest mowa, żeby przedziały, w których grupujemy wartości do osobnego słupka nie były za wąskie. Tak z 10 pomiarów, to dobrze by było, by było. Tu, w fotografii, każda wartość (0, 1, 2, ..., 255) ma osobny słupek (kreskę), a wykresy i tak są dość gładkie. Jak to możliwe? Pomiarów mamy 20 milionów. Tak, aparat fotograficzny potrafi w 0.1 sekundy odczytać ze światłoczułej matrycy i zapisać w pamięci 30MB (cały Sienkiewicz i Żeromski).


Rozkład normalny ↑↑↑

1. Normalne, że najważniejszy rozkład. Opisuje prawie wszystko. Ma też więcej nazw: rozkład Gaussa, krzywa dzwonowa. Dawniej był na banknocie 10DM (Google: "10DM Gauss").

2. Opisuje wszystko bo wszystko jest sumą drobnych czynników. Wpiszmy taką, długą formułę:
=rand()+rand()+rand()+rand()+rand()+rand()+rand()+rand()+rand()+rand()+rand()+rand()-6
Rand'ów jest 12, to nie przypadek, ale na tym etapie nie da się tego wytłumaczyć, pięć by wystarczyło. Natomiast, dlaczego na końcu odejmujemy 6, jest dość proste. Chcemy wynik mieć wokół zera. Popatrzmy (we własnym arkuszu) na serię takich liczb losowych normalnych unormowanych. Zróbmy ich, jak poprzednio 10tys. i wklejmy do histogramu, który, być może mamy już otwarty.

3. Kształt otrzymany w histogramie opisuje funkcja:
f(x)=(2*pi)-1/2*exp(-x2/2) czyli e do -x2.
Ta funkcja f(x) jest gęstością prawdopodobieńtwa, bo pole pod nią jest równe 1. Na poniższym fragmencie banknotu funkcja jest ładniej napisana, ale ponadto ma dwa parametry, zwyczajowo oznaczane greckimi literami μ i σ, o tym później.
To z tym polem powierzchni pod f(x) jest ważne, zróbmy wykres. Użyjemy gotowej funkcji: =NORMDIST(A1,0,1,false). Teraz policzmy małe prostokąciki, których pole wynosi 0.2*0.01=0.002, powinno ich więc być (tych niebieskich) 1/0.002=500 sztuk. Do dzieła.
Oczywiście trzeba to zrobić inaczej. Mając kilkadziesiąt wartości funkcji (dla x = -4, -3.9, -3.0, ... 0, 0.1, ... 4) posumujmy je, a sumę pomnóżmy przez 0.1. Wyjdzie coś takiego: 0.999949. OK! Co zrobiliśmy? Zsumowaliśmy pola pionowych prostokącików. Gdybyśmy wzięli trapeziki, byłoby dokładniej, ALE, przede wszystkim ograniczyliśmy się do przedziału od -4 do 4. W tym przedziale pole, na prawdę jest 0.999936. Nasz błąd numerycznego całkowania wynosi tylko 0.001%.

4. W Wikipedii, w haśle IQ, jest mowa, że wielkość ta podlega rozkładowi normalnemu ze średnią 100 i odchyleniem standardowym 15. To są te parametry na dziesięciomarkówce (i nie tylko): μ=100, σ=15. Następująca formuła:
=NORMDIST(A1,100,15,true)*100 będzie nam odpowiadała na pytanie, jaki procent ludzi ma IQ powyżej A1.
Jeżeli na Ziemi jest 1e10 ludzi, to ten z nich, który ma najwyższe IQ stanowi 1e-8% ludności. Jakie ma IQ? Takie, żeby
=(1-NORMDIST(A1,100,15,true))*100 wyniosło ok. 1e-8. Szukamy A1.


Skala i przesunięcie ↑↑↑

1. Funkcja rand daje nam rozkład równomierny na przedziale (0, 1). Każdy wynik z tego przedziału jest jednakowo prawdopodobny. Przez "wynik" możemy rozumieć liczbę zapisaną z dokładnością do np. czterech cyfr po kropce, np. 0.1234.
Pomijając małą praktyczną przydatność rozkładu równomiernego, powiedzmy, że potrzebny nam jest rozkład na przedziale (25000, 30000). Pewien ankieter musi wylosować 100 respondentów z takiego zakresu pewnej listy. Chodzi o oczywiście o liczby całkowite, o numery, ale to nie ma większego znaczenia, jeżeli jest ich 5000. Po prostu, na pewnym etapie zaokrąglimy wynik.
Rzecz nie jest trudna, jak pomnożymy rand przez 5000, otrzymamy liczby od 0 do 5000. Teraz wystarczy dodać 25000.
=rand()*5000+25000
Formuła =NORMDIST(A1,100,15,false) jest równoważna formule =NORMDIST(A1,0,1,false)*15+100. Sprawdźmy.
Poniższa, używana już tu formuła:
=rand()+rand()+rand()+rand()+rand()+rand()+rand()+rand()+rand()+rand()+rand()+rand()-6
generuje liczbę losową o rozkładzie normalnym właśnie z μ=0 i σ=1, czyli tzw. standardowym. Gdyby składników było 13 sigma byłaby większa itp. Oczywiście, w dalszym ciągu nie wiemy DLACZEGO. Jeżeli tak wygenerowaną liczbę pomnożymy przez dowolne σ i dodamy dowolne μ, to otrzymamy liczbę losową z rozkładu normalnego o dowolnych parametrach μ i σ. Zapisuje się to czasem: N(μ, σ), a czasem N(μ, σ2), to kwestia gustu.

2. Kierunek odwrotny jest taki: z=(x-μ)/σ. Zauważmy, że x, μ i σ mogą mieć wymiar. Na przykład, x może być jakąś długością w mm, a μ może wynosić 15.5mm, a σ=1.2mm. Jednak zmienna z nie ma już wymiaru i "kręci się" w okolicy zera nie oddalając się na dużo bardziej niż 1. Przerabianie x na z (oczywiście nie chodzi o te, czy inne litery) nazywa się czasem standaryzacją.

To, że jak krzywa wyższa, to węższa ilustruje ten programik, działający na ruch kursora (po co ma czerwone buciki, patrz niżej).


Skumulowany rozkład ↑↑↑

1. Skumulowana funkcja rozkładu prawdopodobieństwa, albo krótko, ale mniej zrozumiale: dystrybuanta.
To ta niebieska linia. Zawsze zaczyna się od 0, a kończy na 1. Na osi pionowej ma, po prostu prawdopodobieństwo, żadną tam gęstość. W przeciwieństwie do f(x), zwykle pisze się ją F(x).
Jaki jest związek pomiędzy f(x) i F(x)? Druga jest całką pierwszej, pierwsza pochodną drugiej.
F(x) odpowiada na pytanie jakie jest prawdopodobieństwo, że otrzymamy wartość mniejszą od x. Inaczej: F(x) = P(X < x), gdzie X to zmienna losowa, a x to liczba ;-). Na przykład, jak widać z wykresu (a dotyczy on rozkładu normalnego unormowanego): P(X<-1)=0.16, a P(X<2)=0.977.

2. Zapamiętać by należało trzy liczby dotyczące najważniejszego rozkładu prawdopodobieństwa; rozkładu normalnego. Pierwsza to: F(1)-F(-1) = 0.68, czyta się to tak: prawdopodobieństwo, że zmienna losowa normalna wylosuje się w przedziale między -1 i 1 wynosi 2/3 (dokładniej 0.68). Dla szerszego przedziału, od -2 do 2 jest to już 0.95 (taka ładna liczba), a dla (-3, 3) aż 0.997. To ostatnie, to tzw. reguła 3 sigma, streszcza się ona do tego, że zmienna normalna (czyli niemal każda) nie oddala się od μ dalej niż na 3σ, tzn. oddala się, ale tylko raz na 1/(1-0.997)=370 przypadków (średnio). Całość nazywa się 68–95–99.7 rule - patrz takież hasło Wikipedii.

3. Krzywej kumulacyjnej można użyć do generowania liczb losowych o dowolnym rozkładzie. Dla rozkładu normalnego mamy, na razie tylko sposób z dwunastoma rand'ami. Niezbyt elegancki.
Na powyższym wykresie jest blado niebieska linia. Idzie ona od lewej, od wartości F(x)=0.9 i "szuka" odpowiedniego x, znajdując go w x=1.2. Tak się to robi. Bierze się rand, umieszcza się na osi pionowej i odczytuje się x. To x jest już zmienną losową o rozkładzie normalnym. Jak chwilę pomyśleć, to to działa. Np. żeby otrzymać ekstremalne x=3 lub więcej, trzeba trafić rand'em powyżej 0.998 - to nie jest łatwe. To jest dokładnie tak mało prawdopodobne, jak mało prawdopodobne jest otrzymanie x>3 wg N(0,1). Potrzebna nam tylko funkcja odwrotna do =NORMDIST(A1,0,1,true), tak się nazywa: =NORM.S.INV(A1). Robimy więc tak:
=NORM.S.INV(rand()), a jak potrzebujemy IQ, czy coś innego, to tak:
=NORM.S.INV(rand())*15+100


Rozkład wykładniczy ↑↑↑

1. Normalny jest najważniejszy, ale nie jest jedyny. Tytułowy rozkład tego tematu jest zupełnie inny. Jak ktoś pamięta, co to ex, to ten rozkład, to po prostu e-x, od 0 do ∞. Nasi znajomi rozjechali się po całym świecie, do różnych stref czasowych. Od czasu do czasu, jak im przyjdzie ochota dzwonią do nas. Kiedy dzwoni nasz telefon? O każdej porze. Czy to będzie tak, że jak mamy tych telefonów w ciągu dnia, średnio 24 (powiedzmy), to dzwonek będzie mniej więcej co godzinę? Zupełnie nie. Oni kompletnie się nie umawiają, może i wiedzą o sobie, a nawet do siebie też dzwonią, ale załóżmy, że do nas telefonują niezależnie, jak coś ciekawego zobaczyli, jedni w Hongkongu, inni w Lizbonie. Jak to będzie wyglądało?

2. Zróbmy model w arkuszu. W A1 piszemy =rand(). Ciąąąąągniemy w dół, kilkaset wierszy. Zaznaczamy kolumnę i wybieramy Format i dalej Conditional formatting. Rozwijamy Format rules i wybieramy Greater than, wpisujemy wartość 0.95. Zielone komórki to telefony od znajomych. One się grupują. Chętnie przychodzą jeden po drugim. Może to i dobrze, bo siłą rzeczy muszą też wystąpić dłuższe przerwy, może akurat w nocy? Dawne centrale telefoniczne miały z tym kłopot, żeby obsłużyć niemal jednoczesne żądania. Zamiast telefonów mogą być (inne) nieszczęścia. Przysłowie "Nieszczęścia chodzą parami", jak większość przysłów ma pewne podstawy. Podstawą tego akurat jest rozkład wykładniczy (nie wytłumaczymy tego babci, która mówi "popatrz, sąsiad z góry złamał rękę, a dopiero co wczoraj ciocia Bronia złamała nogę").
Ta wartość przy formatowaniu warunkowym, 0.95 jest dowolna, pierwsza z brzegu. Jeżeli jeden wiersz, to miałaby być minuta, żeby było dokładnie, to przy średnim odstępie między telefonami = 1 godzina, powinno tam być 1-1/60=0.98333. Taka sobie uwaga.

3. Kompletny wzór, z jednym parametrem, T którego sens jest taki, że jest to średni odstęp, jest taki:
f(t) = 1/T * e-t/T dla t>0
zmienna nazywa się tu t, nie x, bo mowa była o czasie, ale to nie ma żadnego znaczenia.
Funkcja wykładnicza bardzo łatwo się całkuje. Całka z f(t), czyli F(t)=1-e-t/T. A to oznacza, że funkcja odwrotna F-1(p)=-T*Ln(1-p), gdzie p to, oczywiście prawdopodobieństwo. A to oznacza, że formuła =-LN(rand()) będzie nam generowała losowy odstęp miedzy niezależnymi zdarzeniami. Dlaczego LN(rand()), a nie LN(1-rand())? Bo wychodzi na to samo. Gdzie T? W tej formule T=1 (może jedna sekunda, może godzina), jak chcemy inne T, to musimy pomnożyć.

Ile będziemy mieli telefonów dziennie? Średnio 24. A ile może być? Tyle: 0, 1, 2, ... 24, ... 100, ... Wynik 24 jest zapewne najbardziej prawdopodobny, 0 jest mało prawdopodobne, podobnie jak 100, a tym bardziej 1000. Jednak nie ma górnej granicy. Tym różni się ta sytuacja od rozkładu dwumianowego. Ta nowa sytuacja, dzienna liczba połączeń, nazywa się ...


Rozkład Poissona ↑↑↑

1. Zauważmy, że ze zmiennej ciągłej, czasu t (albo x) wracamy tu do liczby sztuk k. Litery jak litery, ale sprawa jest istotna. Liczby rzeczywiste, a całkowite. Informatykowi nie jest wszystko jedno.

2. Weźmy w rozkładzie dwumianowym małe prawdopodobieństwo "sukcesu"; p=0.1. Kibicujemy słabej drużynie. Jeżeli jednak prób będzie dużo; N=100, to z 10 sukcesów będzie. (Są słabi, ale bardzo lubią grać.)
Powiedzmy, że interesuje nas liczba sukcesów, ma być 10. Czyli iloczyn N*p musi być równy 10. Można to osiągnąć na wiele sposobów: 100*0.1, 1000*0.01, 10000*0.001 itd. Rozkład dwumianowy będzie niemal identyczny:
BINOM.DIST(10,  100,  0.1, false)=0.13187
BINOM.DIST(10, 1000, 0.01, false)=0.12574
BINOM.DIST(10,10000,0.001, false)=0.12517


Czytelnik, na tym etapie znajomości z Arkuszem i jego Formułami jest w stanie wybaczyć powyższy zapis.
Różnica jest jedna, BINOM.DIST(101, 100, 0.1, false) zwróci błąd, a BINOM.DIST(101, 1000, 0.01, false) zwróci zero (7e-66). Jednak to tylko formalność.
Konsekwentnie napiszmy to tak, do A1 wpiszmy 1000, a do A2: =BINOM.DIST(10,A1,10/A1,false).
Otrzymamy prawdopodobieństwo=0.1257. Teraz, do A3 wpiszmy: =POISSON(10,10,false). Dostajemy 0.1251, drobna różnica.
Funkcja BINOM.DIST ma, matematycznie dwa parametry. Pierwszy parametr "informatyczny" to zmienna k, drugi to właściwy parametr N, trzeci to p, a czwarty decyduje, czy chcemy prawdopodobieństwo skumulowane. Funkcja POISSON (dlaczego nie POISSON.DIST?) ma, matematycznie jeden parametr: N*p. Zwykle oznacza się go grecką lambdą: λ = Np.

3. Wzór jest taki:
p(k) = λk * e / k!

Coś wspólnego z dwumianowym ma ... k! Googlujmy wzór Stirlinga, pamiętając, że N → ∞, a p → 0, tak, żeby N*p było stałe (np. 10).
p(k) = N! / ( k! * (N-k)! ) * pk * (1-p)N-k

4. W przemyśle, nawet z lekką przesadą rozróżniają rozkład dwumianowy od Poissonowskiego. To się tam nazywa karty kontrolne (dział SPC - Statistical Process Control). Chodzi o to, że jak operator bada 100 sztuk produktu i wykrywa tam 4 wadliwe, to to jest rozkład dwumianowy. Bardzo interesuje ich p - prawdopodobieństwo "porażki", czyli wadliwej sztuki. Żeby nie wzrosło nadmiernie. Natomiast jak inny operator ogląda stół i widzi na blacie 4 (drobne) rysy, to ma do czynienia z rozkładem Poissona. Tu nie ma żadnego N, ani p, jest tylko λ - średnia liczba rys (która, oby nie rosła).

5. Znamy trzy ciekawe rozkłady, dwumianowy, Poissona i normalny. Dwa dyskretne, jeden ciągły - ważna różnica, ale pomijamy ją. Weźmy w dwumianowym N=100, p=0.3, w poissonowskim λ=N*p=30, a w normalnym weźmy μ=30 (mniej więcej rozumiemy), natomiast
σ = (N*p*(1-p))1/2 = 4.58 (DLACZEGO TAK?). Zróbmy, ponadto jeszcze jeden normalny, z σ = (N*p)1/2 = 4.48 (też nie wiemy dlaczego?). Tak to wygląda:
Gdybyśmy wzięli N=1000, p=0.03 wszystkie by się zlały, bo 1-p byłoby już bardzo bliskie 1. Wszędzie NORMALNY! Wszędzie, wszystko się sumuje, oddziałuje ze sobą.

6. Czy są jeszcze jakieś inne rozkłady prawdopodobieństwa? Patrz "List of probability distributions". Czy da się je jakoś opisać, jednym, dwoma, góra trzema parametrami? Na przykład, czymś takim jak μ i σ - przesuwanie i przemnażanie (skalowanie)?


Parametry rozkładów ↑↑↑

1. Jeden jest najprostszy, oznacza się go zwykle μ, a nazywa WARTOŚĆ OCZEKIWANA. Ta nazwa dużo mówi. To parametr przesunięcia, pokazuje, czy rozkład jest w środku, po lewej, czy po prawej. Czy losujemy wartości z okolicy 0, czy raczej 100 (jak IQ), albo 165cm (wzrost).

1a. Weźmy dziwną kostkę do gry, która na trzech ściankach ma 4 oczka.
Możliwe wyniki: 1, 2, 3, 4.
Rozkład prawdopodobieństwa: 1 - 1/6, 2 - 1/6, 3 - 1/6, 4 - 1/2.
Średnia z 1, 2, 3, 4 wynosi (1+2+3+4)/4=2.5, inaczej zapisane: 1*1/4 + 2*1/4 + 3*1/4 + 4*1/4. Ale przecież 4 wypada częściej niż pozostałe liczby. Tak ma być:
μ = 1*1/6 + 2*1/6 + 3*1/6 + 4*1/2 = 3.
Wyszło 3, ok!, ale poprzednio 2.5, a nie ma ścianki z dwoma i połową oczka, nieważne. Dla zwykłej kostki μ=3.5, dla monety (orzeł=0, reszka=1) μ=0.5. Oczekiwana wartość może nigdy nie wypadać. Mały paradoks.

1b. Ogólny wzór jest taki:
μ = k1 * p1 + k2 * p2 + k3 * p3 + ile trzeba.

1c. Dwumianowy: μ = N*p. Sprawdźmy w arkuszu. Ambitniejsze zadanie? Policzmy sumę dla k od 0 do N iloczynów k*p(k), gdzie
p(k) = N! / ( k! * (N-k)! ) * pk * (1-p)N-k

Poissona: μ = λ
Normalny: μ = μ ale chwila, jak liczymy sumę "po" x, które jest ciągłe? Liczymy całkę od -∞ do ∞ z iloczynu x*f(x). Czyli nieskończona sumę z nieskończenie małych xf(x)dx. Pamiętamy, że takaż całka z f(x)dx wynosi 1 (jeden).
Wykładniczy: μ = λ     - to akurat względnie łatwo policzyć.
Uwaga o liczeniu całek. Całka (nieoznaczona) z f(x) dla rozkładu normalnego, po prostu nie istnieje. Całkę oznaczoną da się dokładnie policzyć, ale trzeba się napracować, żeby otrzymać to 1. Natomiast całka z x*f(x) jest w zasięgu analizy matematycznej uczelni technicznej.

2. Drugi parametr jest znacznie bardziej skomplikowany niż pierwszy. Nawet nazwę ma pokręconą. Oznaczenia jako, tako stabilne: σ. A nazwy takie: wariancja, dyspersja, odchylenie standardowe. Pierwsza oznacza właściwie σ2. Ostatnia jest nieprawidłowa, ale, jak to bywa, najczęściej używana i sama się trochę objaśnia.

2a. Ogólny wzór dla zmiennej dyskretnej (k) jest taki:
σ2 = (k1 - μ)2 * p1 + (k2 - μ)2 * p2 + (k3 - μ)2 * p3 + ile trzeba.
Mamy tu odchylenia, odchylenia od wartości oczekiwanej. Dlaczego podniesione do kwadratu? Bez tego wyszło by 0, bo część odchyleń jest na +, a część na -. Dlaczego nie wartość bezwzględna? Kwadrat JEST lepszy, chodzi o efektywność estymatora (!).

2b. Tym razem na próbę weźmy już zwykłą kostkę (bo i tak pojęcia nie mamy co ma wyjść). Wiemy natomiast (?), że μ=3.5. Liczmy:
=((1-3.5)^2+(2-3.5)^2+(3-3.5)^2+(4-3.5)^2+(5-3.5)^2+(6-3.5)^2)/6 wychodzi 2.9, koniecznie bierzemy z tego pierwiastek, i mamy σ=1.71. Narysujmy to sobie.

Trzy pionowe linie, środkowa k=μ, boczne k=μ-σ i k=μ+σ. Oś pozioma nazywa się k, są to liczby całkowite, liczby oczek, natomiast μ i σ już takie nie są, czyli pisać k=μ jest niezręcznie. Ten sam kłopot.

2c. Dwumianowy: σ2 = N*p*(1-p). Dla małych p będzie więc w przybliżeniu σ2 = μ. To ciekawy i ważny wynik, przydatny w statystyce. Domyślamy się, wobec tego, że dla rozkładu Poissona dokładnie: σ2 = μ = λ.
Wyprzedzając nieco materiał (wchodząc w statystykę); Ankietowano 1200 respondentów, i 16 odpowiedziało, że lubi X. Czyli p, że ludzie lubią X p=16/1200 (mało, ale nie o to chodzi). Jaka jest dokładność naszego badania? Ma podstawie naszego wyniku sądzimy, że μ = 16, tak więc σ = 4. Czyli dokładność tego 16 jest ±4. Wyszło nam 16, ale mogło wyjść (z inną, podobną grupą) 12 albo 20, albo nawet 8, czy 22, czy jakoś tak. (Statystyka nie daje pewnych odpowiedzi, po prostu nie mamy lepszej metody. Lepiej znać ograniczenia swojej wiedzy.)

2d. Wzór dla zmiennej ciągłej (x) to całka z (x-μ)2f(x).
Dla normalnego wychodzi (domyślamy się) σ2 = σ2.
Dla wykładniczego: σ2 = T2.

2e. Dla równomiernego (rand) policzymy. Oczywiście μ=0.5 - środek symetrycznego rozkładu. Gęstość prawdopodobieństwa jest 1 na odcinku (0, 1), a poza nim 0. Czyli liczymy całkę od 0 do 1 z (x-0.5)^2*1 = (x-0.5)^2 = x^2 - x + 0.25
funkcja pierwotna: 1/3*x^3 - 1/2*x^2 + x/4, podstawiamy x=1: 1/3-1/2+1/4 = (4-6+3)/12 = 1/12, x=0 nie musimy podstawiać.
σ2 = 1/12, σ = 0.29. To jest ta dwunastka z sumy dwunastu rand'ów. Trzeba ich aż tyle, żeby osiągnąć σ2. Bo wariancje (σ2) się sumują. Dlaczego?

3. Czy są jeszcze jakieś ogólne parametry rozkładów? Dlaczego nie?
σ2 = Σ(ki - μ)2 * pi może jakiś sens miałyby jeszcze potęgi 3 i 4? Σ(ki-μ)3*pi. Miałyby. Nazywałoby się to skośność i kurtoza i opisywałoby kształt rozkładu. Skośność wiadomo, np. wykładniczy jest zdecydowanie niesymetryczny, skośny. Gorzej z kurtozą. Zacznijmy od negacji. "Naukowiec.pl" pisze: Kurtoza jest miarą koncentracji wyników. Kurtoza informuje nas o tym, na ile nasze obserwacje, wyniki są skoncentrowane wokół średniej. Ten tekst oczywiście dotyczy parametru σ. Co to więc kurtoza? Kurtoza jest tak zdefiniowana, że dla rozkładu normalnego wynosi 0. Normalny, jak zawsze jest punktem odniesienia. Nie zależy to od σ, od szerokości rozkładu. Kurtoza opisuje kształt. Będzie ona dodatnia, jeżeli rozkład, przy nie zmienionej sigmie będzie miał dłuższe ogony. Oczywiście wtedy musi być smukły w centrum. Kurtoza szybko wykrywa tzw. wartości oddalone. Gdy poza grupa porządnych pomiarów mamy jeden lub kilka wartości "odlecianych", nie wiadomo skąd. Rozkład równomierny, który w ogóle nie ma ogonów ma ujemną kurtozę (-6/5).
Normalny ma tu σ2 = 1/12, tyle co widoczny równomierny. Różnią się tylko kurtozą.


Zależność ↑↑↑

1. Coś musi zależeć od czegoś innego. Potrzebujemy dwóch zmiennych losowych. Może rzut jednocześnie kostka i monetą? Niezbyt ciekawe, ale na początek...

2. Dwuwymiarowy rozkład prawdopodobieństwa. Jak będzie wyglądał w tym przypadku. Wykres musi mieć trzy osie k-kostka, m-moneta, p-prawdopodobieństwo; k=(1, 2, ..., 6), m=(0, 1), p=1/12 - mamy rozkład równomierny. Dwanaście możliwych wyników. Tylko tytułowej zależności tu nie ma. Kostka nie zależy od monety, moneta nie zależy od kostki. Może i się trochę obijają o siebie, ale co z tego wynika? Jak zrobić zależność? Pomijając jakieś skomplikowane mechanizmy w rodzaju magnesów na ściankach i monecie, zróbmy tak.
Rzucamy monetą i kostką, ale jak moneta 1, to kostką rzucamy drugi raz. Rozkład? Piszmy tak: p(m,k)=p, będzie więc:
p(0,1)=1/12, p(0,2)=1/12, ..., p(0,6)=1/12,
p(1,2)=1/(6*6*2)=1/72, p(1,3)=2/72, p(1,4)=3/72, ..., p(1,12)=1/72. Jest tu zależność? Co to jest zależność?
Zależność mamy wtedy, kiedy znajdziemy choć jedną taką parę m,k, że p(m,k1)≠p(m,k2≠k1) lub p(m1,k)≠p(m2≠m1,k). Ta definicja wygląda dość abstrakcyjnie, ale przecież każdy wie, co to jest zależność.
No to aktualnie jak jest? Sprawa prosta p(1,12)=1/72, a p(0,12)=0.
Ale, co zależy od czego? Czy to k zależy od m, czy odwrotnie? Jasne, że k od m, taki jest mechanizm. A jeżeli nie znalibyśmy mechanizmu, a tylko wynikowy rozkład prawdopodobieństwa? Przecież p(0,2)=1/12, a p(1,2)=1/72, też różne.
Zwykle na wykresach rozrzutu X-Y na osi poziomej (x) mamy zmienną niezależną, a na pionowej (y) zależną. Okazuje się, że nasz rozkład możemy narysować tak i tak, m-k albo k-m.
Czytelnik zechce wybaczyć, że oś pionowa na lewym wykresie ma niewłaściwy kierunek, rośnie w dół.

2. Tak to jest z zależnością. Duża część badań (i życiowych research'ów) polega na poszukiwaniu kierunku zależności. Kto od Kogo? Czasem jest tak, że ani x od y, ani y od x, tylko oba od z! A tego z nie znamy, więc szukamy. To sobie zaraz pomodelujemy. Ale najpierw, obowiązkowo ...

3. Dwuwymiarowy normalny. Czyli funkcja gęstości dwuwymiarowego rozkładu normalnego.
f(x)=(2*pi)-1/2*exp(-x2/2) to jeden wymiar (patrz),
f(y)=(2*pi)-1/2*exp(-y2/2) to drugi.
Jeżeli mamy otrzymać x i y (takie, a takie x i takie, a takie y), to musimy pomnożyć gęstości prawdopodobieństwa.
f(x,y)=f(x)*f(y)=1/(2*pi)*exp(-(x2+y2)/2) Poniżej dołączamy 4 należyte parametry:
f(x,y)=1/(2*pi*σxy)*exp(-(((x-μx)/σx)2+((y-μy)/σy)2)/2) Ten programik ilustruje powyższe (4 parametry sobie jeżdżą).

3a. BRAKUJE tu piątego parametru piątej esencji (kwintesencja), piątej kolumny (E), piątego elementu (B (bor)), planety (Jowisz), wymiaru (teoria strun), symfonii (c-moll op. 67), ulicy (NYC), Nostromo (eee, ósmy był), smak (umami), krąg piekła (za Styksem).
Brakuje ZALEŻNOŚCI. Ten programik ją ma (natomiast owe 4 są tu: 0, 0, 1, 1). Piąty parametr nazywa się ρ (ro).
f(x,y)=1/(2*pi*σxy*(1-ρ2)1/2)*exp(-1/(2*(1-ρ2))*(((x-μx)/σx)2+((y-μy)/σy)2)-2*ρ*(x-μx)/σx*(y-μy)/σy))
Jeżeli wyrzucić cztery parametry, które już znamy, będzie krócej.
f(x,y)=exp(-(x2+y2-2*ρ*x*y)/(2*(1-ρ2)))/(2*pi*(1-ρ2)1/2) tu już coś widać. Przede wszystkim dla ρ=0 dostajemy iloczyn f(x)*f(y) - brak zależności. Właściwie, trochę lepiej byłoby mówić o KORELACJI zamiast zależności. Bo tu nie ma żadnego kierunku x→y, ani y→x. One, te dwie zmienne losowe są ze sobą skorelowane (albo nie, ρ=0). Parametr korelacji ρ mieści się w przedziale (-1, 1). Dla ρ=-1 mamy wykres \ dla ρ=1 / dla ρ=0 0 (kółko, elipsa, ale pozioma albo pionowa, nie skośna).

4. Czy może być więcej (skorelowanych lub nie) zmiennych losowych niż dwie? Tak, może.
Tu mamy ich 10. Ciekawe, ale między PM2.5 i PM10 nie widać korelacji. Antykorelację (ρ<0) widać między wilgotnością i wiatrem. Temperatura też się z tymi dwiema koreluje, ale zmienne czystości powietrza, jakieś takie chaotyczne. W tych danych jest autokorelacja, ale to temat szeregów czasowych, później.
Ile mamy współczynników korelacji, jak mamy n zmiennych? Mamy ich n*(n-1)/2, każda zmienna z każdą, ale tylko raz. Jaka jest korelacja zmiennej z samą sobą? Oczywiście ρ=1. Mamy więc macierz korelacji. Jest ona kwadratowa n*n, symetryczna (ρijji), a na głównej przekątnej ma jedynki. Ponadto, jeżeli na przykład ρ2,3=0.9, a ρ3,4=0.8, to ρ2,4 nie może być 0, ujemne, ani bardzo małe.

5. Model korelacji: A1: =NORMSINV(rand()), B1: =A1+0.1*NORMSINV(rand()), i mamy między A1 i B1 silną, dodatnią korelację. Warto to zobaczyć na wykresie rozrzutu, przeciągnąwszy wpierw dwie komórki trochę w dół. B to, po prostu, A plus trochę szumu.


Statystyka


Dlaczego statystyka, to nie matematyka? ↑↑↑

1. Rachunek prawdopodobieństwa jest działem matematyki. Kołmogorow podaje formalną definicję prawdopodobieństwa, ale nie zająknie się, jak to prawdopodobieństwo obliczyć, skąd je wziąć. Moneta jest symetryczna, ok, ale bywają trudniejsze przypadki.
Zadanie z rachunku prawdopodobieństwa: wzrost dzieci w klasach pierwszych podlega rozkładowi normalnemu z wartością oczekiwaną μ = 133 cm i dyspersją (odchyleniem standardowym) σ = 17 cm. Oblicz jakie jest prawdopodobieństwo, że losowe dziecko będzie wyższe niż 150 cm (jaka część dzieci >150 cm). (Zgodnie z wcześniejszą umową, zadanie to rozwiązujemy w pamięci. Bo znamy zasadę 68–95–99.7.)
ALE SKĄD mamy te parametry rozkładu μ i σ i skąd wiemy, że normalny? Dzieci nie są symetryczne. Nic z tych rzeczy. Trzeba je zmierzyć. Zmierzyć wszystkie?
Jeżeli zmierzymy wszystkie obiekty jakimi się interesujemy, nie potrzebujemy statystyki. Wtedy "krzywa kumulacyjna" mówi nam wszystko.
Opis poniższego rysunku. Widoczna formuła wygenerowała nam 30 losowych wzrostów, z odpowiednimi parametrami. Kolumnę A kopiujemy do B jako wartości. Bo musimy ją posortować, a sortowanie liczb, które ciągle się losują, to syzyfowa praca. Sortujemy B (najwyższe dziecko ma 171 cm!?). W C1 wpisujemy =1/60 (obiektów jest 30), w C2 =C1+1/30 i przeciągamy. Zaznaczamy dane w B i C, wstawiamy wykres, Arkusz Googla powinien się domyślić, że chodzi nam o wykres rozrzutu. Czy punkty układają się w kształt podobny do niebieskiej krzywej? Jeżeli tak, to dzieci są "normalne".
Na powyższym wykresie widzimy wszystkie dzieci (w interesującym nas aspekcie), uporządkowane, dla jasności. TYLKO ŻE nam nie chodzi zwykle o konkretne dzieci, a o pierwszoklasistów w ogóle. Może nie o wszystkich na świecie, ale w naszym województwie? Czy jakoś tak.
Statystyka jest ważna tam, gdzie mamy poważne zadania. Np. poważny fotograf potrzebuje histogramu. W produkcji przemysłowej mierzymy produkty nie tylko dlatego, żeby odrzucić te, które nie spełniają kryteriów, ale żeby obliczyć μ i σ, żeby wiedzieć jak ma się nasz proces produkcji (SPC - Statistical Process Control). Czy nie zmierza w niebezpieczną bliskość dopuszczalnych granic. Chcemy to wiedzieć wcześniej, zanim będziemy musieli wyrzucać gotowe produkty.

2. Mamy 30 pomiarów. Mamy konkretne liczby (w kolumnie B). To ile wynosi μ? (Ile σ, czy rozkład JEST normalny?)
To wcale nie jest prosta sprawa. Potrzebujemy jakiejś "oceny" μ, jakiegoś "estymatora".
Posłużmy się pewnym "oszustwem". Osoba, która mierzyła obiekty NIE WIE jakie jest μ, a my wiemy. Mianowicie μ = 133 cm.
Wartość μ powinna być gdzieś w środku kolumny B. No tak, ale środek jest tu, między 15 i 16.
To może wziąć (140.1 + 140.9)/2 = 140.5 ? Taka wielkość nazywa się mediana.
A może tak: 119.4+120.5+122.7+123.3+123.5+124.6+128.7+129.6+129.9+130.9+132.1+135.9+138.4+138.7+140.1 +140.9+142.9+143.8+144.4+146.4+146.5+147.1+148.1+148.8+150.3+153.7+154.1+156.9+161.0+171.0=139.8. To by była średnia.
Czy μ = mediana, albo μ = średnia? Widzimy że nie, 133 ≠ 140.5 ≠ 139.8. Trzy różne liczby.
Na tym właśnie polega kłopot ze statystyką.
Ale my MUSIMY coś wiedzieć o μ. Nigdy nie będziemy znali μ ale COŚ chcemy wiedzieć. Być jak najbliżej μ.
Może "moda" byłaby rozwiązaniem? Ktoś mógł słyszeć takie statystyczne pojęcie. NIE, nie ma idealnego rozwiązania, jest tylko najlepsze.
Średnia jest lepsza, czy mediana (modzie dajmy spokój)? W tym konkretnym przypadku bliżej jest średnia (pamiętajmy o naszym "oszustwie", nielegalnie znamy μ). Ale oba "estymatory" sa bardzo bliskie sobie, a względnie odległe od μ.
To może być przypadek, że średnia jest bliżej. Możemy próbować poanalizować nasze estymatory. Rozpatrzmy taki przypadek, mamy jedną wartość dużo większą od reszty. Biedne dziecko musiałoby mieć 2 metry, ale najczęściej takie wyniki sa skutkiem "grubego" błędu pomiarowego. Otóź mediana jest nieczuła, czy to są 2, 3 czy 5 metrów, a średnia jak najbardziej. Mediana jest dobra, jak mamy "wartości oddalone". No dobrze, skupmy się na czymś konkretnym. Załóżmy, że mamy rozkład normalny. O rozkładzie normalnym można pisać w ∞. Wiemy, że ma on ogony ciągnące się do ± ∞, KAŻDA wartość jest możliwa. Jednak duże odstępstwa są bardzo mało prawdopodobne (reguła 3σ).
W kolumnie A mamy 30 liczb losowych o μ = 133. Dodajmy w A31 =AVERAGE(A1:A30), a w A32 =MEDIAN(A1:A30). Natomiast w A33 wpiszmy =abs(A$31-133)<abs(A$32-133), czy średnia jest bliżej (TRUE) czy mediana (FALSE). I liczmy. Na palcach. (Chyba, że ktoś poradzi sobie z "iterative calculation".) Bliżej częściej będzie średnia. Średnia jest lepsza. Jest estymatorem najefektywniejszy. Efektywność estymatora jest trochę inaczej zdefiniowana. Otóż średnia jest lepsza bo σśrednia < σmediana Odchylenie standardowe średniej jest nieco mniejsze niż odchylenie standardowe mediany. Przy tym, co bardzo ważne, wartości oczekiwane obu estymatorów są równe 133 (w tym wypadku), czyli tyle co trzeba. Jak średnia jest lepsza, to po co nam mediana. Była już o tym mowa. Jest lepsza na okazję wartości odstających.


Średnia ↑↑↑

1. Poznaliśmy już średnią, i medianę, AVERAGE i MEDIAN. Średnią zwykle oznacza się kreską nad literą oznaczającą zmienną losową, czyli, na przykład x z kreską u góry. Ta kreska jest kłopotliwa, ale da się zrobić: x̄ (kod #772). Można też pisać xśr .

2. Średnia to też zmienna losowa. (Co to "zmienna losowa"?)
Trzydzieścioro pierwszoklasistów z Ib miało wzrost jak widać. (Jasio był nieobecny, a on jest malutki.) W Ic byłoby inaczej. Dlaczego mierzyliśmy Ib? Bo akurat mieli WF. To jest losowość. A zmienna = liczba. To wszystko (prawie wszystko).
Jeżeli 30 liczb jest losowych (jeszcze błędy pomiarowe dochodzą, one przecież nie ustoją spokojnie) to jakaś funkcja tych liczb też będzie losowa. Ta funkcja to =SUM(A1:A30)/30 albo =AVERAGE(A1:A30) albo śr = 1/n*Σxi albo
.

2b. Suma, to suma. Było o tym trochę tutaj, w punkcie 3. Pojawiło się tam pojęcie splotu, to trudna sprawa. Otóż mamy tu sumę zmiennych normalnych. Jeżeli suma zmiennych równomiernych (rand) daje normalny, to suma zmiennych normalnych już na pewno ma rozkład normalny. PRAWDA. Średnia ma rozkład normalny, nawet jak x ma jakiś inny, ale najczęściej i tak ma normalny. U nas x ma rozkład N(133 cm, 17 cm) albo ogólnie N(μx, σx).
No to średnia ma rozkład N(μśr, σśr), tak możemy sobie napisać, ale jak te parametry mają się do tych wcześniejszych?
Sprawa jest poważna, skupmy się na istocie rzeczy, na sumowaniu. Powiedzmy, że x ma rozkład N(0,1), najprostszy. Pytanie: jaki będzie rozkład x+x (może lepiej zapisać to x1+x2, żeby komuś 2x nie przyszło do głowy, bo chodzi o =rand()+rand()).
Zacznijmy od
Wzór do przeklejenia mamy tu, w punkcie 3, jednak trzeba go przepisać, bo Wolfram nie czyta indeksów górnych:
(2*pi)^(-1/2)*e^(-x^2/2)
powyższe wklejamy jako pierwszą i drugą funkcję. To jest N(0,1). Wynik mamy pod Exact result:
(2*2*pi)^(-1/2)*e^(-x^2/(2*2)) Wyboldowana dwójka, to po prostu σ2, która się pojawiła na skutek sumowania (splatania).
(2*pi*σ2)^(-1/2)*e^(-x^2/(2*σ2)) - to jest N(0,σ), a nam w splocie ta σ wyszła równa 21/2 = 1.41.
Dlaczego to musi być takie skomplikowane? Nic nie wyjaśnimy zanim nie zajrzymy do czarnej skrzynki Wolfram. Weźmy arkusz.
Widzimy tu nasze dwa sumowane Gaussy. I jakiś trzeci malutki. Te dwie funkcje nie musiałyby być takie same. Rozważamy, po prostu najprostszy przypadek sumowania dwóch normalnych zmiennych losowych. Na obrazku akurat suma wynosi 2 (w A1 jest wpisywana ręcznie). No właśnie, co jest w arkuszu. W kolumnie B jest, co widać, od -6 do 6 co 0.2. C1: =NORMDIST(B1, 0,1,false) itd., D1: =NORMDIST(B1, A$1,1,false), E1: =C1*D1, F1: =sum(E1:E61)*A6 - jedna komórka. Jeżeli ktoś jest tu przypadkiem to: 1. Używamy Arkusza Googla, w Excelu byłoby niemal tak samo. 2. Do komórek wkopiowujemy to "=NORMDIST(B1, 0,1,false)" bez cudzysłowu. i przeciągamy w dół (61 wierszy). 3. Dolar w A$1 oznacza, że mimo przeciągania w dół, zawsze będzie to adres A$1, czyli A1. Aha, 4. w B1 wpisujemy -6, w B2 -5.8, zaznaczamy obie komórki , łapiemy za mały kwadracik, kursor staje się krzyżykiem, i ciągniemy w dół.
Ten niewymiarowy (nie ma pola 1), fioletowy Gauss jest iloczynem (=C1*D1) niebieskiego i czerwonego. Jego pole powierzchni (całka) liczona jest w F1. Teraz, gdy chcemy mieć x+x=2, pole to wynosi 0.104. I to właśnie jest wynikowa gęstość zmiennej y=x+x (y jak u Wolframa). Tak, punkt po punkcie możemy sprawdzać, czy poniższa formuła da to, co mamy w F1 (zamiast y wpisujemy A1)
=(4*pi())^(-1/2)*exp(-(A1^2)/4)
DLACZEGO tak? Dlaczego pole pod iloczynem? To jest tak. Chcemy, żeby suma była równa y. Jak to osiągnąć? Ano z pierwszego rozkładu dowolne x, ale z drugiego y-x, żeby x+(y-x)=y. Dowolne x oznacza, że musimy posumować, scałkować iloczyn gęstości prawdopodobieństwa. Iloczyn, bo ma być x i y-x. Łatwo się googlują ładne ruchome wizualizacje splotu.

2c. Zróbmy nowy punkt, choć to ciągle o sumie dwóch N(0,1). Otóż ta suma ma rozkład N(0,σ2=2). Wygodniej użyć wariancji σ2, żeby było, po prostu 2, a nie √2. Poza tym łatwiej sformułować tezę:
Przy sumowaniu dwóch zmiennych losowych normalnych (ale nie tylko) wartości oczekiwane i wariancje się sumują.
Skąd wiemy, że μsuma = μ1 + μ2 ? Z takiego samego rozumowania, jak to które doprowadziło do σsuma2 = σ12 + σ22 , tylko prostszego. Jak ma być inaczej, jak postawimy dziecko na dziecku? Wartość oczekiwana "wzrostu" będzie 133 cm + 133 cm = 266 cm.

2d. Tak więc, suma n zmiennych x podlegających rozkładowi N(μ, σ2) ma rozkład N(n*μ, n*σ2).
Drobna uwaga. Jaki rozkład ma: y=a*x+b? Taki ma: N(a*μ+b, (a*σ)2). Nie może być inaczej, przesuwanie rozkładu (b) nie ma wpływu na jego szerokość (σ). Natomiast zmiana skali tak. Jak dzieci zaczniemy mierzyć metrach, to σ=17cm zamieni się na σ=0.17m, to żadna statsytyka. Nowe μ=1.33m, a jak dzieci będą (przez pomyłkę) stały na metrowym pudle, to będzie μ=2.33m, a σ się nie zmieni (no chyba, że pudło się kołysze i wprowadza dodatkową zmienność, która wynosząc 10cm, czyli σpudło=0.1m, spowoduje, że w pomiarach wyjdzie σ=(0.172+0.12)1/2=0.2m). Zwróćmy uwagę, że wariancja nie ma takiego wymiaru jak x, μ, czy σ, ma kwadrat tego wymiaru. Jeżeli σ=17cm, to σ2=290cm2. Niby banał, a ile pomyłek. Pomijając już tę, że ktoś tu zarzuci, że 172=289. Ale to jest inna kwestia, raczej trudna.

2e. Tak więc suma n zmiennych podzielona przez n, czyli średnia, ma rozkład N(μ, σ2/n). Czyli σśr=σ/√n . Powoli: σ2 sumy = n*σ2, czyli σ sumy = √n*σ, czyli σ średniej = √n*σ/n = σ/√n, cbdo.
To fundamentalny wynik, ten obrazek nie jest za duży!


Odchylenie standardowe ↑↑↑

1. W temacie parametry rozkładów są definicje μ i σ. Wiemy już jak oceniać (estymować) μ mając serię pomiarów (serię obserwacji zmiennej losowej). Estymatorem wartości oczekiwanej, czyli μ jest średnia (mediana też). Co z σ?

2. Drugi parametr, σ2 zdefiniowany jest jako wartość oczekiwana kwadratu różnicy zmiennej losowej minus μ, czyli takiej wielkości: (x-μ)2. To może estymatorem tego parametru byłaby średnia z (x-xśr)2? To dobry pomysł. Choć ma pewną wadę. Lepiej byłoby użyć (x-μ), tyle, że μ nie znamy, a xśr sobie policzymy. Musimy pozostać przy (x-xśr)2. Za jaką cenę (czy to jest jakaś strata?)? Wzór jest taki: s2 = 1/(n-1)*Σ(x-μ)2. To nie jest średnia, sumę n składników dzielimy przez n-1! DLACZEGO? Bo tak jest dobrze. Gdybyśmy dzielili przez n, otrzymywalibyśmy wielkość nieco za małą. Dla n=2, 3, 4 to ma znaczenie. Dlaczego tak jest? To kara za użycie we wzorze wielkości obliczonej z danych, które tu pojawiają się drugi raz. Mamy o jeden "stopień swobody" mniej. To żaden kłopot, żadna strata, tak jest i już.

3. Estymatorem wariancji czyli wielkości σ2 jest zdefiniowane przed chwilą s2. To dobrze, że używa się innej litery, to JEST coś innego! Natomiast s=√s2 (że tak to napiszemy) nazywa się odchyleniem standardowym. Tu niestety drobna uwaga. s2 jest "nieobciążonym" estymatorem σ2, tzn. wartość oczekiwana estymatora równa się estymowanemu parametrowi. Niestety, nie dotyczy to s, które otrzymujemy z nieliniowego przekształcenia (pierwiastek) s2. W przemyśle (SPC) stosują odpowiednie korekty. Jednak jest to pewna przesada, gdyż dokładność s nie jest wysoka i upieranie się przy kolejnych cyfrach po przecinku nie ma sensu.

4. Zobaczmy to wszystko (średnią i odchylenie standardowe) w arkuszu. Wypełnijmy ze 30 wierszy i kolumny od A do Z najprostszą zmienną normalną N(0,1), czyli =NORM.S.INV(rand()).
W A32 policzmy średnią =AVERAGE(A1:A30), w A33 odchylenie standardowe =STDEV(A1:A30). Obie formuły przeciągamy do Z.
W A35 policzmy średnią średnich =AVERAGE(A32:Z32)
W B35 policzmy odchylenie standardowe średnich =STDEV(A32:Z32). Przeciągnijmy obie formuły z wiersza 35 do wiersza 36. Co tu się policzyło?
W A36 mamy średnią odchyleń standardowych.
W B36 mamy odchylenie standardowe odchyleń standardowych.
W A35 "powinno" wyjść 0, w A36 1, bo takie parametry są w NORM.S.INV(rand()), μ=0, σ=1. Gorzej z odchyleniem standardowym średnich. Ale to wiemy. Powinno wynosić 1/√30=0.183, zgadza się? Bardzo źle jest z odchyleniem standardowym odchyleń standardowych. To nawet trudno wygooglować. Powinno być 1/√(2*30)=0.129, zgadza się?
Jeszcze takie sprawdzenia, liczby w wierszu 32 nie powinny odchylać się od zera o dużo więcej niż ich σ=0.183.
Liczba w A35 nie powinna być odchylona od zera o dużo więcej niż jej σ=1/√(30*26)=0.036. Dlaczego 30*26? Czego mamy tu 26?

5. To jest moment aby powiedzieć wyraźnie, co robimy mając 30 pomiarów pewnej cechy (wzrost pierwszaków).
Liczymy średnią (xśr), liczymy odchylenie standardowe (s). To drugie dzielimy przez pierwiastek z n (u nas n=30).
Mamy dwie liczby 139.803 i 2.2728. Wszystko w porządku, tylko tych cyferek za dużo. Nie wolno pisać losowych cyferek i podawać je za prawdę. Niepewności pomiarowej nie podajemy z większą liczbą cyfr znaczących jak 2. Tu to jest omówione i jest programik, który nam niepotrzebne cyfry obetnie. W naszym przypadku ma być tak: 2.3. Z liczbą cyferek przy średniej sprawa jest prostsza. Piszemy tyle po kropce, ile jest w niepewności, a więc: 139.8. Końcowy wynik:

139.8 ± 2.3 cm

Jeszce dwa słowa o liczbie cyfr. Powiedzmy, że mamy wynik: 123456.78 ± 444.55. Nikt mi nie uwierzy, ale powinien on być tak zapisany 123460 ± 440.


Prawda statystyczna ↑↑↑

1. Zmierzyliśmy dzieci, wyszło 139.8 ± 2.3 cm. Dyrektor mówi: co takie małe, średnia krajowa jest 144 cm. To 144 też ma swój błąd, pewnie mały, bo dużo dzieci pomierzono, ale ma. O tym raczej się nie myśli.
Jaka jest prawda? Prawda statystyczna, ale to jedyna jaką możemy mieć. Taka jest.
Policzmy (144-139.8)/2.3 = 1.83. Stosując regułę 68–95–99.7 wiemy, że prawdopodobieństwo odchylenia większego od 2σ od średniej wynosi (100-95%)=5%=0.05=1/20. To już, owszem, dość małe prawdopodobieństwo, p=0.05 uważane jest często za graniczną wartość. Że coś, co jest mniej prawdopodobne, jest już podejrzane, jeśli wystąpiło. My mamy odchylenie tylko 1.8, mniej niż 2. Taka więc jest prawda: Nasze odchylenie od zadanej wartości (144) jest nieistotne statystycznie. Jest przypadkowe. NIE można mówić, że nasze dzieci są niższe. Gdyby mały Jasiek nie był akurat chory, moglibyśmy przekroczyć próg p=0.05, ale było jak było. Nie wolno jednego robić. Małego Jaśka wysłać do domu BO mamy pomiary wzrostu. To nieuczciwe. Jeżeli ktoś wie co to "brzytwa Ockhama", to takie działanie nazywa się miotłą Ockhama (Ockham's broom).
Zastosowaliśmy regułę 68–95–99.7. Można by też postąpić inaczej, policzyć =1-NORMDIST(144,139.8,2.3,true), wynik jest 0.034. To jest prawdopodobieństwo. Ono jest mniejsze od 0.05. Czyli nasze dzieci jednak SĄ mniejsze? Sprawa jest trochę pokomplikowana. Zacznijmy od narysowania reguły 68–95–99.7 (99.7 nie będzie widać).
Otóż na to 5% składają się dwa trójkąciki ciemno niebieskie. Gdyby one "szły" od 1.8 byłyby większe i razem miałyby 6.8% (2*0.034 w %). Wyraźnie powyżej 5%. Ale dlaczego brać oba? Poprawne rozumowanie jest takie. Jeżeli nasze dzieci zawsze by małe i dyrekcja wyłącznie pyta, czy ciągle są mniejsze. W ogóle nie słuchając, gdyby były wyższe. To wtedy bierzemy jeden trójkąt i mamy 3.4% < 5% i mamy te dzieci statystycznie, istotnie mniejsze. Natomiast, jeżeli pytają nas, czy dzieci mieszczą się w średniej krajowej, żeby tylko nie były za małe ANI za duże. To bierzemy oba ogony rozkładu i mamy 6.8% > 5%. Jest OK.


Test t ↑↑↑

1. W poprzednim temacie podzieliliśmy różnicę między średnią a zadaną wartością przez sśr (=2.3cm).
z = (xśr - x0) / sśr. To x0, to jakaś zadana wartość, że niby nasza średnia ma być blisko. Pozostałe dwie rzeczy po prawej stronie znaku równości, to obliczone przez nas estymatory. Wynik oznaczony jest przez z, o tym była już mowa. Z "powinno" mieć rozkład normalny N(0,1). Jeżeli rzeczywiście wartość oczekiwana x-ów, z których policzyliśmy xśr wynosi x0 (μ=x0), to xśr - x0 ma rozkład N(0,σ2/n). Jeżeli podzielilibyśmy (xśr - x0)/σ, to otrzymany iloraz jest z-tem; ma rozkład N(0,1). Ale my w mianowniku mamy s, nie σ! Na pewno s≠σ Bo s jest estymatorem, jest zmienną losową. Prawdopodobieństwo, że dwie liczby rzeczywiste będą przypadkiem równe wynosi 0, jak wiemy, znając pojęcie gęstości prawdopodobieństwa.

2. Wielkość t = (xśr - x0) / sśr ma pewien rozkład prawdopodobieństwa, ale nie jest to rozkład normalny, nazywa się on rozkładem t, lub t-Studenta. Dla dużych n, nie różni się on bardzo od normalnego N(0,1). Posiada μ =0, a σ2 = m/(m-2), gdzie m jest liczbą stopni swobody, m=n-2, więc u nas m=28. Czyli σ=1.04, o 4% większa od 1. Kształt rozkładu t też podobny do normalnego.
Policzmy =1-T.DIST((144-139.8)/2.3,28,true), otrzymamy 0.039. Poprzednio, używając NORMDIST otrzymaliśmy 0.034.
Abstrahując od wzrostu pierwszaków, sprawdźmy =T.DIST(2,A1,true)-T.DIST(-2,A1,true), przy A1=30. Wychodzi 0.945. Co my tu liczymy? Prawdopodobieństwo między -2 i 2, które, wg normalnego wynosi, jak wiemy 0.95 (dokładniej 0.9545). Zgadza się dość dobrze. Jak w A1 wpiszemy 100, będzie jeszcze lepiej: 0.952. Natomiast jeżeli mamy tylko pięć pomiarów, czyli n=5, czyli m(=A1)=3, to okazuje się, że p(-2<t<2)=0.86. Sprawdźmy to w arkuszu. Arkusz Google ma kolumn od A do Z, a wierszy 1000, dopóki nie poprosimy o więcej. Jak robimy eksperyment Monte Carlo, to z 1000 prób by sie przydało. Oznacza to, że nasze małe, 5-elementowe próbki trzeba wpisać w wierszu. Komórki A1:E1 wypełniamy znaną formułą =NORMSINV(rand()). Mamy więc μ=0 - nie będziemy musieli odejmować. Do F1 idzie formuła =AVERAGE(A1:E1)/STDEV(A1:E1)*2.236. Co to za tajemnicza liczba na końcu? To pierwiastek z 5, a dlaczego tu jest? Bo autor jest starej daty i nie może przełknąć, żeby biedny komputer musiał 1000 razy liczyć pierwiastek. To jest przyczyna, że tu jest liczba. Ale po co ten pierwiastek tam? O tym było dość dużo, dopiero co. Ok, tam musi być odchylenie standardowe średniej, a nie x-ów. Idziemy dalej, G1: =And(F1>-2,F1<2). Łapiemy i ciągniemy w dół. Oczywiście można to zrobić inaczej, z IQ>130. Teraz, H1: =countif(G1:G1000,true). Zgadza się z 0.86?

3. T.TEST. Ten punkt poświęcony będzie czytaniu ze zrozumieniem. Sprawę testu statystycznego zostawmy na później. Googlujemy "google spreadsheet t.test". Dostajemy "Google Support":
=T.TEST(A1:A4,B1:B4,2,1) możemy sobie od razu wpisać do naszego arkusza. Ponaciskajmy klawisz Delete. Co wychodzi? To samo, co przy =RAND(), liczba losowa (0, 1). Dlaczego? Bo w A i B mamy to samo (w sensie statystycznym). Wynik formuły T.TEST to p-wartość (p-value) - zostawmy to. Czytajmy Help. Napisane jest, że zwracana wartość jest "związana" z testem??? Żeby uniknąć hasła "p-value". Zdanie "Determines..." jest mniej więcej ok, ale... . Pytamy, czy "dwie" próbki pochodzą z tych samych "dwóch" rozkładów? Czyli pierwsza pochodzi z dwóch i druga z tych samych dwóch? Idźmy dalej. Syntax, tailes. Określamy tu liczbę ogonów rozkładu. Ale rozkład t ma oczywiście dwa ogony, lewy i prawy. Chodzi o test jednostronny lub dwustronny. To JEST trudne do wyjaśnienia. Jest o tym w końcówce poprzedniego tematu.
Typ testu. Są trzy. O co chodzi w "sparowanym". Robimy badania kliniczne, testujemy przez tydzień diety A i B, badamy spadek masy ciała. Jeżeli mamy mnóstwo chętnych (i dieta nie jest kosztowna), to połowie dajemy A, połowie B i po tygodniu robimy t.test typu 2, nie 1. Natomiast przy małej liczbie osób robimy inaczej, najpierw dajemy im A lub B (losowo), po tygodniu wysyłamy na urlop, a potem dajemy im te drugie diety. Na końcu mamy, dla każdej z osób dwa wyniki, dla A i dla B. To są te pary. (A1, B1), (A2, B2) itd. Zmienność między osobami jest duża i zaciemni nam obraz, dlatego lepiej, żeby każda osoba przeszła obie terapie (o ile to możliwe). Test bierze wtedy średnią z różnic, a nie różnicę średnich.
Typ 2 i 3, to dość prosta sprawa. Robiąc test interesujemy się wartościami oczekiwanymi, pytamy, czy μAB? Nie interesują nas σ. Natomiast, jeżeli WIEMY, że σAB, to tę wiedzę trzeba wykorzystać. Test będzie działał wtedy trochę lepiej. Jeżeli włączamy jakąś naszą wiedzę a priori, to dobrze jest ją sprawdzić. Jest osobny test na weryfikację hipotezy σAB, nazywa się FTEST. Rozważanie, czy większym błędem będzie użycie typu 2 przy (nieuniknionej) niepewności wiedzy o sigmach, czy też użycie typu 3 i zmarnowanie tej wiedzy (być może prawdziwej), musiałoby włączyć ocenę prawdopodobieństwa fałszywości wiedzy a priori. Natomiast, w arkuszu, wpisać raz 2, a potem 3 i zobaczyć wyniki, nie jest trudno. Wybrać właściwy jest trudno (patrz "miotła Ockhama", w poprzednim temacie).


Test statystyczny ↑↑↑

1. Rozważamy ofertę wycieczki do Neapolu.
Są dwa wyniki tego rozważania. Jedziemy, albo nie jedziemy. Jeden bit. Są też dwa możliwe błędy do popełnienia.
Błąd I rodzaju: pojechaliśmy, a nie było warto.
Błąd II rodzaju: nie pojechaliśmy, a ci, co pojechali mówią nam, że błąd.
Są też dwa przypadki słusznej decyzji, ogólnie tak to wygląda:

               ocena   ocena
rzeczywistość   ok     błąd II
rzeczywistość  błąd I   ok


Prawdopodobieństwo popełnienia błędu I rodzaju oznaczamy α, II rodzaju, β.
Koszt błędu II rodzaju jest zwykle wyższy. Nie pojechać i potem żałować... a pojechać i zastać strajk muzealników... (dlatego wycieczka była taka tania).

2. Mieliśmy hipotezę "zobaczyć Neapol i umrzeć", to była nasza wstępna hipoteza, nazwijmy ją H0. Była też alternatywna hipoteza, H1: nie warto.
Mamy zmienną losową normalną N(μ,σ). A przynajmniej tak nam się wydaje. Skupmy się na μ. Wydaje nam się (jesteśmy przekonani), że μ=μ0. To jest nasza hipoteza zerowa, H0. Nasze pierwszaki mają μ=μ0=144 cm. Ale może nie mają (są większe ALBO mniejsze - hipoteza alternatywna dwustronna, są mniejsze - hipoteza jednostronna).
Skądś wiemy (na pewno), że σ=1. Stawiamy hipotezę (najprostszą): μ0=0. Przeciwko hipotezie alternatywnej: μ0≠0. To brzmi jakoś niezbyt mądrze. Co, jeżeli μ0=0.0001? Nic specjalnego, test nie wykaże nam "prawdy". Przy małym odchyleniu od H0 mamy duże β, duże prawdopodobieństwo błędu II rodzaju (ale też jego szkodliwość nie jest wtedy wielka). Testujemy hipotezę, robimy pomiar (my ograniczymy się do rand()). Wychodzi nam 0.6. Jasne, że ok, mieści się nawet w przedziale (-1, 1). Nawet gdyby wyszło -1.2. Też musielibyśmy się cieszyć, że nie ma podstaw odrzucenia H0, bo dopiero prawdopodobieństwo wyjścia poza przedział (-2, 2) jest małe, wynosi 0.05. Czyli 2.1 kazałoby nam odrzucić H0. Najczęściej p=0.05 przyjmuje się jako graniczną wartość. Nazywa się ją poziomem istotności testu. To jest właśnie α. Arbitralna wartość.

3. W poprzednim temacie robiliśmy T.TEST. Funkcja ta zwracała liczbę, a przecież test kończy się jednobitowo: jechać / nie jechać. W zasadzie, dopiero to jest test:
=T.TEST(A1:A4,B1:B4,1,1)>0.05, prawda oznacza: nie odrzucać H0.

4a. W temacie Rozkład dwumianowy, mimochodem pojawiło się pytanie, ile razy trzeba by rzucać krzywą monetą (p=0.52, na przykład), żeby stwierdzić, że nie należy ona do uczciwych (odrzucić H0: p=0.5). Rozpatrzmy sprawę jakościowo, rysunkowo. Rzućmy tą monetą 10000 razy. Narysujmy sobie dwa rozkłady prawdopodobieństwa:
=BINOMDIST(A1,10000,0.5,false), i drugi, z 0.52 zamiast 0.5. Wiemy, że σ=(10000*0.5*0.5)1/2=50. Tak więc w A1 wpisujemy 5000-150=4850 i ze skokiem 5 ciągniemy w dół, do 5200+150. Otrzymujemy taki obrazek.
Widzimy, że przy N=10000 rzutach mamy szanse. Rozkłady są rozdzielone. Małe jest prawdopodobieństwo by "czerwona" moneta dała wynik, który zaakceptowalibyśmy, zakładając H0, czyli niebieski rozkład.

4b. Dochodzimy do sprawy prawdopodobieństwa błędu II rodzaju β. Z powyższego rysunku możemy je odczytać. Najpierw jednak musimy się umówić, co do testu. Powiedzmy, że odrzucimy H0 jeżeli dostaniemy wynik 5110 lub więcej (jakie wtedy będzie α? będzie sumą niebieskich słupków dla k>=5110). Jakie jest β, czyli prawdopodobieństwo, że wg czerwonego rozkładu otrzymamy jednak mniej niż 5110 orłów? To suma czerwonych słupków. Policzone. Mamy α i β, wiemy wszystko o teście statystycznym. Jak te rzeczy NA PRAWDĘ policzyć? Zmienić w formułach false na true, i wykonać jakieś odejmowania.
Oczywiście β będzie inne, dla innego p, my zbadaliśmy tylko przypadek p=0.52. Wykres 1-β w funkcji parametru rozkładu to funkcja mocy testu. Moc testu zależy od α. Im większe α, tym większa moc. Czyli: większe prawdopodobieństwo błędu I rodzaju oznacza mniejsze prawdopodobieństwo błędu II rodzaju. I odwrotnie. Skąd my to znamy? Albo albo. Optymalizacja α wymaga znajomości kosztów w obu przypadkach. W produkcji przemysłowej to jest do zrobienia. W nauce i w życiu codziennym, raczej nie. W SPC pojęcie mocy testu jest dobrze znane.

4c. Wykreślmy moc najprostszego testu. Najprostszy test jest taki:
H0: N(0, 1), H1: N(μ>0, 1). Wydaje nam się, że nasza zmienna losowa ma rozkład normalny o parametrach μ=0, σ=1, ale podejrzewamy, że jednak μ może być większe (test jednostronny, jednoogonowy). Poziom istotności testu ustalamy na α=0.05, standardowo. Tak więc, graniczną wartością będzie:
xgr =NORMSINV(1-0.05) = 1.65. Popatrzmy na ten rysunek. Bladą kreską narysowany tam jest (przypadkowo) przypadek α=0.1 (xgr=1.2). Wykres funkcji mocy, na osi poziomej ma μ, od zera, do, powiedzmy 5 (A1=0, A2=0.1, ..., Aileśtam=5). Na osi pionowej ma prawdopodobieństwo
1-β =1-NORMSDIST(1.65-A1).
Jest to, po prostu kawałek tej niebieskiej krzywej kumulacyjnej. Zwróćmy uwagę na wartość dla μ=0; wynosi ona α=0.05. Bo μ=0; oznacza, że H1=H0, a w takim wypadku 1-β=α. Prawda o Neapolu jest taka, że jest paskudny, ale poziom paskudności wynosi 0, czyli rzeczywistość = rzeczywistość. Za to przy paskudności=5 prawdopodobieństwo popełnienia błędu drugiego rodzaju, czyli oceny, że trzeba jechać, jest praktycznie zerowe. Wykres jest jasny, tylko jak wyznaczyć μ, przeglądając ofertę turystyczną?
W przemyśle to robią, i tak to może wyglądać (źródło):


Test χ2 ↑↑↑

0. Ta literka podnoszona do kwadratu, to nie X, to jest Chi. W Times, i większe, jest takie (przy okazji inne nasze): χ α β μ ρ σ Σ .

1. Wróćmy do histogramu. Tym razem zróbmy go jednak w arkuszu, żeby krok po kroku wszystko kontrolować.
1. A1:A300 =NORMSINV(rand()) - nasza próba losowa, dość duża.
2. B - śmieszna rzecz, choć dość naturalna. Pewne funkcje nie chcą mieć za argument niepewnych rzeczy w rodzaju RAND czy NOW, tzn. takich, co się ciągle zmieniają. Dlatego musimy do B przekopiować całe A, ALE wklejamy je "specjalnie", jako wartości (jest na to skrót: Ctrl+Shift+V). W końcu, gdyby to były porządne DANE, a nie nasze "modele"...
3. C1:C15 - wpiszmy takie liczby: -9, -3, -2.5, -2, ... 2.5, 3, 9. Te dziewiątki, to granice, poza które nic nam nie wyjdzie, a pozostałe liczby ... to granice przedziałów, w których będziemy zliczali "pomiary".
4. Do zliczania mamy taki programik:
function hist(x0,p0) {
 var p=new Array(), h=new Array(), n=x0.length, m=p0.length
 for(var i=0;i<m;i++) { p[i]=Number(p0[i]); h[i]=0}; h[0]="-"
  for(var i=0;ii<n;i++) { x=Number(x0[i])
   for(var j=1;ji<m;j++) { if (x>p[j-1] && xi<=p[j]) { h[j]++ } } }
 return h }

który kopiujemy do edytora wywołanego: Tools / Script editor. Zapisujemy program (Ctrl+S), po czym powinien być dostępny w arkuszu. Wywołujemy go w D1: =hist(B1:B300,C1:C15), powinniśmy otrzymać 14 liczb zliczeń, a na początku kreseczkę. Mamy histogram. Możemy go sobie narysować.

2. Test χ2. Co testujemy? Testujemy zgodność naszych danych z rozkładem N(0,1), lub jakimkolwiek innym, ale teraz ten. Zauważmy, że mamy konkretny rozkład, funkcja i parametry. Gdybyśmy jakieś inne dane chcieli przetestować na okoliczność zgodności z N(μ, σ), to wystarczyłoby zestandaryzować je z=(x-μ)/σ w kolumnie B, i wszystko dalej działa. Właściwie, to testować będziemy histogram. Test będzie się odwoływał tylko do kolumny D.
Idźmy dalej: E2 =(NORM.S.DIST(C2)-NORM.S.DIST(C1))*300 i przeciągamy do E15. Co to jest? To są, po prostu prawdopodobieństwa wg N(0,1) wpadnięcia do przedziałów (C1, C2) itd. Jeszcze pomnożone przez 300. Tak, że powinniśmy w E mieć liczby podobne jak w D. Liczby w E są wartościami oczekiwanymi zmiennych losowych z D.
Jaki rozkład prawdopodobieństwa mają całkowite liczby w D? Dwumianowy mają, p=NORM.S.DIST(C2)-NORM.S.DIST(C1), N=300. Sukcesem jest trafienie do danego przedziału. Zauważmy, że p jest małe, bo przedziałów jest trochę. Maksymalnie p=57.4/300=0.2. Powiedzmy, że to małe p. Test χ2 jest bardzo stary, dawniej nie tak sie to liczyło, stąd mnóstwo uproszczeń. Mimo wszystko dokończmy zapoznawanie się ze starym, dobrym statystycznym testem zgodności.
Następne uproszczenie jest takie, zamiast rozkładu dwumianowego używać będziemy normalnego, który jest dość podobny, o ile Np (kolumna E) nie jest zbyt małe, co najmniej 10 żeby było. Dlatego w F i G połączymy po 3 skrajne komórki. Z bidą jest. To akurat żadne oszustwo. Liczniejsze komórki, po prostu przepisujemy. Żal łączyć dalej, bo skrajne wartości mają wartość diagnostyczną.
Rozkład normalny zastępujący dwumianowy jest taki: N(μ=Np, σ2=Np). Drugi parametr piszemy jako σ, albo jako σ2, jak nam wygodnie. Byleby nie było nieporozumienia. Tak więc taka zmienna: z2=(k-Np)2/Np jest kwadratem zmiennej N(0,1). Dlaczego z2, a nie z? Bo interesuje nas tylko wielkość odstępstw niebieskich kwadracików od czerwonych kółeczek, a nie ich kierunek (góra, dół). Teraz możemy te z2 zsumować, i to właśnie jest zmienną losową χ2
χ2n jest sumą n kwadratów zmiennych normalnych unormowanych N(0,1).
H4: =(F4-G4)^2/G4 itd. H1: =sum(H4:H13) - to jest właśnie χ2n. Tylko ile wynosi n? Otóż n=9, o jeden mniej niż jest składników w H4:H13. Dlaczego o jeden mniej. Bo definicja χ2n jest taka, że to jest suma n NIEZALEŻNYCH zmiennych losowych. A nasze F4:F13 są powiązane ze sobą zależnością sum(F4:F13)=300.
wpiszmy jeszcze do H2 (wolne jest): =CHIDIST(H1,9). To jest wynik testu CHI2.

To jest p-value, w H2, niebieskie. Czyli wyszedł nam wynik (16.40), którego prawdopodobieństwo, wobec H0 (a ona JEST prawdziwa, bo to nie żadne dane, tylko model) wynosi 0.059. Nie za dużo, ale >0.05. Czyli jest OK, zgadza się. Tak wyszło. Najgorszy jest słupek histogramu (nie żaden słupek, tylko kwadracik - tak wyszło) dla przedziału (1.5, 2). Na wykresie on jest nad 2. Ale któryś musi być najmniej zgodny z gładkim rozkładem.
Gdyby wyszło 17, albo 20 p-value (wynik =CHIDIST(H1,9)) byłoby jeszcze mniejsze. A jak wpisać w miejsce argumentu H1 1, albo nawet 5, to otrzymujemy liczbę bardzo bliską 1, same dziewiątki. A przecież różne DISTribution nie tak wyglądają, tylko odwrotnie, rosną od 0 do 1 z lewej na prawą. Celowo tak to zrobili. Przy testowaniu zmiennej normalnej czasem odejmowaliśmy 1-NORMSDIST, bo interesował nas prawy ogon. Tu, to odejmowanie zawiera się już w funkcji, bo sens ma tylko za duże χ2. Za duża niezgodność. Za dobra zgodność, wskazywana przez małe χ2, jest, co najwyżej podejrzana. Że może ktoś manipulował danymi (miotła Ockhama) i przesadził.

3. Rozkład χ2. Oczywiście jest niesymetryczny, bo dziedzina (0, ∞), bo to suma kwadratów - skośność √(8/n). Zapamiętajmy wartość oczekiwaną, bo łatwo μ = n i przydaje się do szybkiej oceny wyniku (16.4 od razu coś duże było przy 9). Ciekawa rzecz, oczekiwana wartość sumy n kwadratów liczb N(0,1) wynosi n. To jest okazja do przetestowania Arkusza, ile obliczeń wytrzyma. Wpiszmy do B1:Z1000 =NORMSINV(rand())^2. Potem w A1 wpiszmy =SUM(B1:Z1000), co powinno dać 25000±160 (skąd to 160? patrz następne zdanie). Jeżelibyśmy jeszcze zapamiętali, że σ2=2n, to wiemy już o rozkładzie wszystko.


Regresja ↑↑↑

0. To kontynuacja, na polu statystyki, czyli analizy danych, tematu Zależność.

1. Zacznijmy jednak od współczynnika korelacji. Tym razem chodzi o estymator parametru ρ. Estymator ten zwykle oznaczany jest przez r, czasem R. Popatrzmy jeszcze raz na funkcję gęstości prawdopodobieństwa dwuwymiarowego rozkładu normalnego. Ten programik pokazuje ją z góry. Na osi poziomej jest x, na pionowej y, a f(x,y) ilustrowane jest kolorem. Parametr ρ zmienia sie tu płynnie miedzy -1 i 1. Zaczyna od ρ=0, wtedy jest kółko. Jaka jest różnica między kółkiem, a obrazem jaki mamy dla ρ=0.9 (który możemy sobie wyobrazić, wiedząc, że -1 i 1 to kreski ( \ / )). Taka jest, w prostych słowach. Jak jest kółko, to cztery ćwiartki kwadratu wykresu są tak samo zapełnione kolorem. A przy ρ=0.9 niemal całe niezerowe f(x,y) będzie w ćwiartkach (x>0, y>0) i (x<0, y<0). Co to oznacza? Że jak weźmiemy iloczyn x*y, to (prawie) zawsze będzie on dodatni.
Powtórzmy punkt 5 z Zależność, wpiszmy do A1: =NORMSINV(rand()), B1: =A1+0.1*NORMSINV(rand()) i przeciągnijmy obie komórki do wiersza 100. Następnie, gdzieś obok =CORREL(A1:A100,B1:B100). Wyjdzie dużo, 0.995. Nic dziwnego, suma iloczynów, same dodatnie, jak kilka ujemnych, to małe, bo z okolicy początku układu współrzędnych (nie zrobiliśmy wykresu?). Oczywiście tę sumę dzielimy przez n=100. Z dokładnie poprzedniego tematu wiemy, że suma kwadratów unormowanej zmiennej normalnej wynosi n. A przecież u nas B1 to niemal A1, czyli A1*B1≈A12. Drobny szczegół, nie wiemy jak CORREL się liczy. Porównajmy to z =SUM(C1)/100, po tym, jak do C wpiszemy =A1*B1. Nie wyjdzie dokładnie tak samo, bo MY WIEMY, że μx=0, μy=0, σx=1 i σy=1. Ale CORREL tego nie wie i zrobi tak:
r=Σ(xi-xśr)*(yi-yśr)/(sx*sy)
Pomodyfikujmy to 0.1 w B1 i patrzmy na r. Dodajmy też minus przed A1 w tej formule (dlaczego dodawanie - przed NORMSINV nie ma ŻADNEGO sensu?).

2. Regresja, dziwne słowo (jak Histogram). Chodzi o to, JAK y zależy od x. Zaczynamy oczywiście od najprostszej regresji - liniowej. A więc zakładamy, że y=ax+b. Mamy serię pomiarów (dwuwymiarowych):
(x1, y1), (x2, y2), ..., (xn, yn) i chcielibyśmy z niej "wydobyć" a i b.
Im więcej pada, tym lepiej idą parasole. Proste. Te dane pewnie z klimatu bez śnieżnego (egzotyka), miesiąc ignorujemy.
82 15
92.5 25
83.2 17
97.7 28
131.9 41
141.3 47
165.4 50
140 46
126.7 37

urwało się. Separatorem jest spacja. Kopiujemy i wklejamy do arkusza. To się DA zrobić.

Robimy wykres i dodajemy do niego "linię trendu" Edit, Customize, Series, Trendline, Label=Use Equation.
a = 0.432 1/mm, b = -17. Po prostu, mamy wzór na obliczanie ile zamówić parasoli wobec prognozy opadu x mm.
Jedno małe ALE, do arkusza. Mianowicie a i b to estymatory, mają swoją niepewność. Nie została ona podana.
Inna uwaga, subtelna, ale... względnie ważna. Taka regresja, jak i te dane, to nie jest przykład dwuwymiarowego rozkładu normalnego. Tu (w tego typu regresji) zakłada się, że xi są dokładnie znane, a tylko yi są zmiennymi losowymi (o rozkładzie normalnym). To się zgadza, jako tako z naszymi danymi. Opad jest mierzony dość dokładnie (na stacji meteo), a liczba klientów... losowa (ale od opadu zależna).
Wróćmy do kwestii niepewności "pomiarowych" a i b. Wpiszmy do kolumny z x'ami 1, 2, ..., 9, a za y wpiszmy =NORMSINV(rand()). Tu nie ma żadnej zależności. A przypadkowy wynik jest taki:
Punkty, jak punkty, wylosowane. Policzone osobno r=0.46, takie sobie, mało punktów. Ale JAK ON MOŻE wypisywać równanie regresji, z całkiem dokładnie podanymi parametrami a i b, jak tu żadnej zależności nie ma? Może to jest tak, że to MY WIEMY, że nie ma zależności, bośmy sami ten model zbudowali. Ale ON nie wie. Ale on ma r! Może wygooglować sobie (arkusz Google nie może?) "Table of Critical Values: Pearson Correlation" i odczytać tam, że dla n=0 i p=0.05 mamy r0=0.602 (dwuogonowo - bo przecież nie wiemy, czy korelacja miałaby być dodatnia, czy ujemna). No to wyraźnie nieistotne r wyszło (0.46<0.60), to NIE WOLNO rysować prostej regresji i wypisywać wzoru! Wiadomo jak to jest, to nie Arkusz jest głupi. On jest tylko dostosowany. Typowy użytkownik nie będzie sobie zaprzątał głowy, dlaczega czasem "linia trendu" się rysuje, a czasem nie. Ma się rysować i już.

2a. Bazując na Wikipedii zróbmy sobie to sami. Dane mamy w A (x) i B (y). Co musimy policzyć? Najpierw średnie. C1: =AVERAGE(A1:A9) i przeciągamy do D. Mamy średnie. Teraz E1: =(A1-C$1)*(B1-D$1) F1: =(A1-C$1)^2 G1: =(B1-D$1)^2 i przeciągamy E-G w dół.
Pod spodem, może w E11: =SUM(E1:E9) i przeciągamy w prawo. W C i D mamy trochę miejsca.
C3: =E11/F11 -to jest a, zgadza się z wykresem, z wzorem na linię trendu?
C4: =D1-C3*C1 -to jest b. Teraz niepewności tych a i b. Najpierw musimy policzyć reszty (ε), a właściwie ich kwadraty.
H1: =B1-C$4-C$3*A1 i przeciągamy w dół, a potem dociągamy w prawo G11. Teraz te niepewności i sa:
D3: =(H11/(9-2)/F11)^0.5 -to jest sa. Trochę głupawo wygląda, w formule, to 9-2. Dwa jest zawsze, a 9 to n. Ten nasz arkusz nie jest bardzo trudny do ougólnienia, do dowolnego n. Głównie wiersz 11 trzeba relokować. Żeby policzyć sb potrzebujemy (w kolumnie I) kwadratów A1^2, a potem ich sumy, która będzie w I11.
D4: =D3*(I11/9)^0.5
Policzmy coś bardziej sensownego niż brak zależności. Wstawmy do B =A1+NORMSINV(rand()), żeby była zależność.

Niebieskie liczby w arkuszu to a ± sa i niżej b ± sb te czerwone ± dopisano na bitmapie. Jak widać, a i b zgadzają się z tym co policzył wykres. Natomiast sa i sb nie sprawdzimy, bo nie ma z czym porównać. Musimy popatrzeć na wykres i ocenić czy a i b mogą się sensownie różnić od podanych wartości o swoje niepewności. To nie jest łatwe.

2b. A skądże to, jakże to, czemu tak gna? A co to to, co to to, kto to tak pcha?
Skądże te wzory? Z czego wynikły? Ta para gorąca, to zasada minimalizacji kwadratów odchyleń.
Patrzmy na powyższy wykres. Wartości x, celowo tak proste 1, 2, ..., 9, są zadane, a wartości y są wynikami pomiarów, nie są dokładne, są zmiennymi losowymi (teraz, to one już SĄ, czyli są "realizacjami" zmiennych losowych). Kropki nie leżą idealnie na dopasowanej prostej. Należałoby się dziwić, gdyby tak było.
εi = yi - ( a*xi + b ), to są odchylenia w kierunku pionowym (jak widać z wzoru) punktów od prostej. ε5 jest największe. Parametry prostej (a, b) są tak dobrane, żeby suma kwadratów tych odchyleń Σεi2 była jak najmniejsza. Dlaczego kwadraty? Prosta odpowiedź: żeby się dodatnie i ujemne εi nie odejmowały. Wtedy wszystko można by dopasować. Dlaczego nie moduły, abs, |εi|? Z kwadratami jest lepiej, one są bardziej czułe na duże odchylenia. Weźmy
ε1=1 i ε2=10 wtedy:
ε12=1 a ε22=100. Komentarz zbędny.
Użyte w naszym arkuszu wzory na a i b nawet dość łatwo jest wyprowadzić z zasady minimalizacji (pochodną liczymy) kwadratów odchyleń.

3. Regresja nieliniowa. Użyjmy zasady najmniejszych kwadratów w sposób bezpośredni, do trudniejszego, nieliniowego zagadnienia.
Nieliniowe, czyli jakie? Nie y=ax+b, to może y=ax2. Może nasze dane mają taką właśnie właściwość? Może x to odległość lampy od ściany, a y to natężenia światła (mierzone aparatem fotograficznym). Co prawda, wtedy mamy y=ax-2, ale przykład jest. Jednak parabola, a nawet hiperbola, to banał. Po prostu bierzemy √y, albo jeszcze √(1/y) i używamy naszego, gotowego arkusza - regresji liniowej.
BTW zadanie z fizyki, które wreszcie dostarczyło by nam PRAWDZIWYCH danych. Bierzemy ścianę, lampę i aparat fotograficzny. Lampa 30 cm od ściany. Ustawiamy taką przysłonę i ISO, żeby dostać krótki czas migawki, może 1/1000. Oddalamy lampę. Patrzymy jaki jest czas migawki. I tak z 10 przypadków. Mamy odległość i czas migawki. 1. Wykreślamy to w arkuszu (co jest zmienną niezależną x, a co zależną y?). 2. Wymyślamy jaka powinna być relacja. 3. Próbujemy wymyśloną funkcję dodać do wykresu.

3a. Weźmy sinus, tym bardziej, że sprawa jest całkiem praktyczna.
Widać sinus? Oś Ziemi kołysze się harmonicznie. Wcale się nie kołysze, tylko Słońce jest coraz to z innej strony. Owszem kołysze się też, ale tego nie widzimy, bo okres tego ruchu jest 25 tysięcy lat. Odczytujemy temperatury: -3 -1.5 1 7.5 13.5 17 18 17.5 13 8 3 -1.
Ile parametrów ma sinus? Cztery: skala i przesunięcie na osi x, skala i przesunięcie na osi y. Prosta ma dwa, zaproponowana wyżej parabola mogłaby mieć trzy, ale miała jeden (y=ax2).
y = by + ay*sin(ax*x - bx)
Cztery parametry do estymowania, to dużo. Bo zachodzi pewna, bardzo niekorzystna sytuacja. Jeszcze raz, te ax, bx, ay, by trzeba tak dobrać, żeby Σεi2 było możliwie najmniejsze. Zadanie jest jasne. To szukamy ax. Wymyślamy jakieś wartości i sprawdzamy sumę kwadratów. Jak dobrze wymyślamy, to dość szybko zaczniemy zacieśniać granice, osaczać optymalne ax. I dalej, bierzemy się za bx. Tylko zaraz, poprzednio jakieś bx musieliśmy mieć do wzoru wstawione. Obojętnie jakie? Nic z tych rzeczy. Sprawa wygląda tak, musimy JEDNOCZEŚNIE przeszukiwać wszystkie parametry. Zmiana każdego może popsuć optymalizację pozostałych. Jeżeli dla każdego parametru mamy 100 propozycji wartości, to Σεi2 trzeba policzyć 1004=1e8, sto milionów razy. Nie jest aż tak źle. Bo działając inteligentnie weźmiemy tylko po 10 wartości. Znajdziemy minimum. I przeszukujemy następne 104pozycji, ale już w dziesięciokrotnie zawężonych zakresach. I tak jeszcze raz, dwa razy. Obyśmy tylko w lokalne minimum nie wpadli. Widzimy dołek, szukamy dokładnie jego środka, a obok jest inne, głębszy dołek.
Jednak cztery parametry to za dużo. Kombinujemy (myślimy). Przecież ax jest oczywiste! Jedna rzecz tu jest pewna, ten sinus musi mieć okres 12 miesięcy. Czyli =SIN(2*PI()/12*A1), gdzie A1 to x=numer miesiąca. Dla A1=0 i A1=12 okresowa funkcja o okresie 12 daje nam to samo. Tak ma być. Policzmy średnią temperaturę, 7.75 °C. Dlaczego by to nie miało być by? Jak będzie źle, to to zobaczymy potem. Tak więc mamy już: =7.75+SIN(2*PI()/12*A1). Pozostałe parametry pozostawiamy do ustalenia, formuła jest wiec taka:
=7.75+C$1*SIN(2*PI()/12*(A1-C$2)), do C1 możemy wpisać 1, ale chyba lepiej będzie 10, a do C2 możemy spokojnie wpisać 0. Z tym C$2 trochę w ostatniej formule pokombinowano. Chodzi o to, żeby parametr mieć w miesiącach, a nie w radianach.
Wygląda to tak, dane, model i w C dwa parametry.
W E, jak widać są kwadraty różnic (ε2). W C3 ich średnia (średnia, suma bez różnicy, ma być minimum). W C3, żeby było blisko parametrów, nad którymi trzeba pracować. Oj trzeba, szczególnie nad przesunięciem w poziomie (C2). Amplituda się przypadkiem udała (10 °C).
Warto się tym pobawić, wynik jest piękny, dopasowanie jakby nawet za dobre, ale przecież dane tak autentycznie wyglądają.
Nie jest to pewnie optimum, ale jest dobrze: 10.7, 4.1, co daje średni kwadrat 0.22. Ponieważ √0.22 = 0.5 °C średni błąd naszego modelu jest pół stopnia, jakby nieduży. Co to jest model? Tu danych mamy niewiele, 12 punktów, a model wygląda tak: y=7.75+10.7*SIN(2*PI()/12*(x-4.1)), ma 4 parametry, zawsze to mniej niż 12. Jednak danych mogłoby być znacznie więcej, ale ważne jest to, że zamiast grupy liczb mamy właśnie "model", który coś nam mówi. Że temperatura w roku zmienia się sinusoidalnie, że przesunięcie w fazie jest ok. 4 miesiące, że amplituda jest ok. 11 °C, oraz, że okres rzeczywiście jest 12 miesięcy, bo krzywe są niemal doskonale nałożone na siebie (co zobaczy, kto zrobił).


Szeregi czasowe ↑↑↑

1. Pierwszy wykres w tym tekście, w temacie Chaos deterministyczny, to szereg czasowy. Bo kolejność liczb generowanych wg wzoru xi+1=4*xi*(1-xi) jest ważna. Wskazuje na to sama struktura wzoru; nowe x otrzymujemy z poprzedniego. Ostatni punkt ostatniego tematu, gdzie jest sinus i kolejne miesiące, opisuje szereg czasowy. Natomiast nieco wcześniejszy, pozornie podobny temat, ze sprzedażą parasoli w różnych miesiącach nie potrzebował pojęcia szeregu czasowego. Miesiące były tam właściwie tylko identyfikatorami przypadków. Nie wchodziły do analizy. Na wykresie rozrzutu zawierającym prostą regresji zupełnie nie wiadomo jaka jest kolejność punktów wg miesięcy. Podobnie było z wcześniejszym przykładem ze wzrostem dzieci. Konkretne wzrosty są w tabeli wpisane w jakiejś kolejności. Kolejne dziecko podchodzi i jest mierzone, a wynik zapisywany. Jaka jest kolejność dzieci, nie ma jednak żadnego znaczenia. Zupełnie inaczej jest tam, gdzie obserwujemy jakiś proces (np. w produkcji, SPC - Statistical Process Control). Gdzie coś się DZIEJE. Operator bierze z taśmy produkcyjnej 5 sztuk produktu do kontrolnych pomiarów. Ta piątka nie tworzy szeregu czasowego, bo mu się pomieszały w kieszeni. Ale mogły mu się pomieszać, bo one były sąsiednie, dzielił je czas, ale krótki. Jak za godzinę przyjdzie po następną piątkę, to te dwie piątki nie mogą mu się już pomieszać. Jeżeli chce, jak najszybciej wykryć ewentualne zmiany. Bo coś się dzieje.
Ważna uwaga! Kolejność, to jedno. Nie wolno pomieszać kolejności. Ale, ponadto, między próbkami/pomiarami muszą być jednakowe odstępy czasu. Nie wolno zgubić pomiaru. Jeżeli, z jakiegoś ważnego powodu nie mamy danej liczby, musimy ją jakoś zastąpić. Dobrym sposobem jest wpisać średnią z sąsiednich wartości. Trzeba tylko ten fakt odnotować.
Poniżej jest najważniejszy szereg czasowy świata (trochę obecnie zdominowany przez inny, ale tu jest geoinformatyka, nie medycyna).

2. Przyjrzyjmy się temperaturom z ostatnich 25 lat. Tu mamy je z rozdzielczością dobową. Dzięki uprzejmości The University of Dayton. Po przeliczeniu z °F na °C, utworzeniu porządnej daty i zaopiekowaniu braków danych (missing data) mamy czym się zająć.
Ten programik pozwala przejrzeć dane graficznie. Mamy tu niemal 10 tysięcy liczb (dni). Prawdziwych danych.
Import danych z Dayton Univ. albo wydłubanie ich z powyższego programu (muszą tam być), to ciekawe zadanie, jednak informatyczne, nie statystyczne. Poniżej jest wykres średnich rocznych temperatur, czyli średnich z 365 dni, z kolejnych lat znalezionych w omawianych danych.
W Warszawie też jest coraz cieplej.

3. Zróbmy model szeregu czasowego. To nic nowego, stare, dobre:
=NORMSINV(rand()). Wpiszmy tę formułę do A1:A128. Tak, 128 liczb. Może ktoś się domyśla dlaczego potęga dwójki. W sumie, co za różnica.

Co to jest? NIC, a nazywa się biały szum. Dlaczego biały? Zostawmy to na razie. Jak zrobić, żeby coś było? Zrobić łatwo, ale co ma być? Niech będzie to, co w najważniejszym szeregu czasowym świata! Wzrost! Mógłby też być spadek. Ogólna nazwa na oba, to TREND. Jak robimy trend rosnący?
=NORMSINV(rand())+row()/30. Co to jest row()? Jak to co? Numer wiersza. Dlaczego dzielimy go przez 30? Żeby trend nie zdominował całego szeregu, żeby szumu trochę zostało, żeby nie było zbyt łatwo z analizą. Po co dodano linię trendu? Tak sobie. Bo ma sens (na poprzednim wykresie nie miała, tam linią trendu jest oś y=0).

Co by innego. Mieliśmy już szereg czasowy, w poprzednim temacie, trochę przypadkiem. Sinusoidalny był. Te nasze warszawskie dane (z Dayton) też takie są. Wiemy już jak się robi sinus. Niech ma okres 64.
=NORMSINV(rand())+2*SIN(2*PI()/64*ROW()), mnożymy, go przez 2, żeby w szumie nam nie zniknął.

A zróbmy coś innego zamiast dodawania do szumu jakiegoś deterministycznego sygnału. Zróbmy tak:
A1: =NORMSINV(rand())     A2: =NORMSINV(rand())+A1    i przeciągamy do A128

Co to jest? Ni pies ni wydra, ni sinus ni trend. Tu jest AUTOKORELACJA. Korelację już generowaliśmy. Mieliśmy wtedy dwie zmienne. Gdzie y to był x z dodatkiem szumu. Tu, kolejne x to jest poprzednie x z większym lub mniejszym dodatkiem szumu. To już nie jest biały szum, wpada w czerwień (nie widać?).
- to razem daje biel. A jak słabnie prawa część robi się żółto (jak dawne żarówki) i czerwonawo. Po prawej stronie tej tęczy mamy szybką zmienność (krótkie fale), po lewej wolną zmienność (długie fale). Na ostatnim wykresie nie ma szybkiej zmienności (jest jej mało), jest wolna. Dołek po lewej trwa niemal ćwierć szeregu.
Jak obrócić kota ogonem? W A2 dać minus zamiast plusa (i przeciągnąć).

Okres tego dodawanego wyżej sinusa możemy zrobić jaki chcemy, możemy kilka różnych sinusów dodać na raz. Tak jak to się przecież jakoś robi w tych różnych elektronicznych, oczywiście cyfrowych keybordach. Temat częstotliwości jest tu ważny, stąd aż dwa obrazki do kojarzenia z tematem.
źródło

4. Wiemy już jak modelować różne cechy szeregów. Ale jak je wykrywać w nowym, nieznanym, zadanym szeregu czasowym? Oto jest pytanie.

4a. Zacznijmy od pozbywania się szumu. Może wtedy łatwiej zobaczymy na wykresie właściwy (ciekawy) sygnał.
źródło
Jako szeregu do wygładzania (coś takiego trzeba zrobić) użyjmy kombinacji szumu, sinusa i trendu.
=NORMSINV(rand())+row()/30+SIN(2*PI()/64*ROW())    może sinusa trochę mniej.
Jaki jest sposób na nierówności, na szum? Papier ścierny? Nie tu. Średnia jest sposobem. Zmniejsza szum √n razy. Jak mamy 9 pomiarów, to średnia ma odchylenie standardowe 3 razy mniejsze (niż pojedynczy pomiar).
W A1:A128 mamy szereg. Policzmy średnią z 9 punktów. W B5 wpisujemy =AVERAGE(A1:A9) i ciągniemy do B124 (uwaga, żadnych $). Dlaczego nie od A1 do A128?

Czerwona linia nie obejmuje całej długości szeregu, to drobiazg, coś byśmy tu wymyślili. Gorzej, że to wygładzenie nie taki bardzo gładkie. Moglibyśmy zrobić średnią ruchomą -tak to się nazywa, nie 9-punktową, a 20 punktową, albo i 3 punktową. Sprawdźmy wynik.

4b. Okno. On tam chyba ramę okna szlifuje? Zastosowaliśmy "okno" o długości 9, i jedziemy nim. Nagle, w pewnym momencie zahaczamy o ekstremalny punkt. Pojawi się w oknie nagle, będzie miał udział w średniej taki sam jak wtedy, gdy będzie w samym środku.
Tu mamy rzeczywiście ekstremalne maksimum. Czy tak powinno wyglądać jego wygładzanie? Trochę tak, ale dlaczego ten "postument" tak się ostro zaczyna i kończy? Nie mógłby gładko, jak dzwon (Gauss)?

4c. Średnia ważona. Pojawia nam się nowy temat, potrzebny na wielu różnych polach, nie tylko w szeregach czasowych.
Średnia to coś takiego xśr = Σ xi/n = x1/9 + x2/9 + ... +x9/9. Dlaczego tak śmiesznie zapisane? Tam gdzie x to wzrost dzieci z Ib, indeks przy x ma tylko formalne znaczenie. W szeregach jest inaczej. Tu x5 jest w środku, jest najważniejsze, a x1 i x9 ledwo, ledwo. Czyli jak? Jakoś tak: xśr = Σ wi*xi, gdzie wi to wagi, które wcale nie muszą być takie same =1/9. Może do b5 wpiszemy:
=0.02*A1+0.05*A2+0.11*A3+0.20*A4+0.24*A5+0.20*A6+0.11*A7+0.05*A8+0.02*A9 zapis ten nie wygląda specjalnie inteligentnie. Wagi są takie jak na tym obrazku:
Gauss, to, nie Gauss? Ważne jest jedno SUMA wag musi być równa 1! Wynik jest taki jak niżej. Tylko jedna uwaga. Poprzednio okno miało szerokość 9 i teraz, niby też. Jednak prostokąt o szerokości 9 jest zdecydowanie szerszy od tego dzwonu, który 9 ma dopiero razem z ogonami.
Wygładzanie na tym wykresie jest słabsze niż na poprzednim, ale jest ładniejsze. Przy czym, oczywiście, siłę wygładzania robimy jaką chcemy, tylko te wagi musimy umieć jakoś produkować.
Wagi robimy tak. Niech okno ma szerokość N=20. W którejś kolumnie, od góry: =EXP(-0.5*((ROW()-10.5)/3)^2) w którejś komórce liczymy sumę tych 20 liczb, kopiujemy ją i wklejamy do formuły =EXP(-0.5*((ROW()-10.5)/3)^2)/7.5138 Teraz suma powinna już być 1. Mamy 20 wag. Dlaczego w formule mamy 10.5 zamiast N/2? Żeby wagi były symetryczne. Muszą być? Nie muszą, ale to już inny temat. O otrzymanym właśnie gaussowskim oknie mówimy, że ma σ=3, bo przecież ma, popatrzmy na formułę.
To te wagi: 0.0009, 0.0024, 0.0058, 0.0127, 0.0248, 0.0432, 0.0674, 0.0940, 0.1175, 0.1313, 0.1313, 0.1175, 0.0940, 0.0674, 0.0432, 0.0248, 0.0127, 0.0058, 0.0024, 0.0009
Jeszcze sprawa końcówek, początku i końca szeregu. Przy N=20 dziesięć pierwszych i ostatnich punktów się nie policzy, bo jak? Będzie odwołanie do komórki A0, A-1 itd. Jak to zrobić. Nieistniejących komórek nie liczymy, proste. Tylko o jednym trzeba pamiętać. Suma wag ciągle musi wynosić 1. Mamy wagi od w1 do w20. Jak chcemy użyć, w pewnym momencie tylko od w5 do w20, to musimy te 16 wag tak przeliczyć, żeby ich suma też była 1. Takie rzeczy robimy już raczej programistycznie.

5. Wykrywanie trendu. Na tym wykresie mamy trend (rosnący) i, od razu podpowiedź, jak go wykrywać. Zrobić regresję, ale porządnie, w posiadanym arkuszu. Otrzymujemy parametr nachylenia a i jego niepewność sa. Jeżeli |a|>3*sa, to mamy trend.
Możemy też użyć testu t. Trzeba podzielić szereg na połowy, albo wziąć początek (może 1/3) i z końca tyle samo i zrobić test. Czy się różnią, co do średniej wartości.

5a. Wykrywanie autokorelacji. Trend, to też autokorelacja, sinus też, jakakolwiek prawidłowość. Na tym wykresie (nad tęczą) nie ma żadnej regularnej prawidłowości, ale jest silna autokorelacja. Na tym niżej, jest anty-autokorelacja, czyli ujemna autokorelacja - tego nie będzie dotyczyła metodo, do której przechodzimy.
Wróćmy w arkuszu do białego szumu do A1:A128 dajemy =NORMSINV(rand()) . Zostawmy obok szeregu wolną kolumnę, a gdzieś z boku policzmy =STDEV(A1:A128), powinno wyjść 1±1/√128. Byłoby dobrze gdybyśmy rozumieli ten ostatni zapis, to ważna rzecz (jeżeli nie najważniejsza w całym kursie statystyki).
Od B2 zróbmy tak: =(A2-A1)^2 i przeciągamy do B128, liczymy średnią =AVERAGE(B2:B128). Podnieśmy jeszcze do kwadratu =STDEV(A1:A128)^2 i podzielmy przez to średnią kwadratów różnic kolejnych pomiarów. Jak się zgubiliśmy, to policzmy gdzieś =AVERAGE(B2:B128)/STDEV(A1:A128)^2. Powinno wyjść dwa (oczywiście około 2). Okazuje się, że średnia kwadratów odchyleń od średniej wychodzi dwa razy mniejsza od średniej kwadratów różnic sąsiednich wartości. Do średniej bliżej niż do (przypadkowego) sąsiada. Niezależnie od μ i σ.
Zauważmy, że STDEV nie zależy od kolejności w danych, a suma kwadratów różnic sąsiednich wartości oczywiście tak. Spróbujmy do białego szumu dodać cokolwiek, ROW(), SIN, czy autokorelację. Parametr, który liczymy natychmiast będzie <2.
Musimy przyznać, że w tym co zrobiliśmy nie ma nic trudnego. Formułę na te kolejne różnice zrobiliśmy sami i nie była skomplikowana. STDEV umielibyśmy sami policzyć, chyba to nawet robiliśmy (jak nie, to warto). Mamy więc proste, nadające się do zrozumienia narzędzie, które działa. Nie zawsze tak jest.

5c. Dlaczego nasz szereg, którym tu się bawimy ma długość 128? To nie jest aż tak ważne, ale sygnalizuje osobny temat, który dobrze się googluje hasłem FFT, a omawiany jest tu.