Rachunek prawdopodobieństwa i statystyka Home

Kontakt: a@adamwalanus.pl; temat maila: GIN. Unikamy załączników, a na pewno ich nie kompresujemy; nazwa pliku załącznika: IAD Nazwisko Imię numer albo inny identyfikator

Zadania

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ść     P-wartość    

Statystyka

Dlaczego statystyka, to nie matematyka?     Średnia     Odchylenie standardowe     Prawda statystyczna     Test t     Test statystyczny     Test χ²     Regresja     Szeregi czasowe     R    

Drugi stopień    Transformata Fouriera w R    Wnioskowanie Bayesowskie

Różne    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 zapisać formułę albo wartość. Ponieważ za chwilę będziemy zmieniać B1, musimy zapisywać (kolejno w C) wartości. Możemy nawet ręcznie je przepisywać. W zasadzie, to w C możemy od razu wpisać, od góry 0.01, 0.001, 0.0001. Czyli w C1 0.01, a w C2: =C1/10 i tę formułę ciągniemy w dół (jakieś 15 wierszy).

7c. Zmieniamy liczbę w B1 na 0.201. Wykonujemy 7b. Czyli w D2 wpisujemy numer wiersza, przy którym wykresy się "rozjeżdżają". I tak dalej. Dodajemy dalsze zera (0.2001), do oporu (jakiego rodzaju opór napotkamy?).

8a. Tworzymy wykres rozrzutu (x,y) z kolumn C i D. Będzie on nieciekawy, dlatego...

8b. W kolumnie E liczymy logarytm dziesiętny z C (wynik będzie prosty, żaden logarytm nie jest potrzebny), 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)! ) * p^k * (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/dm³.
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(-x²/2) czyli e do -x².
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(μ, σ²), 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 e^x, 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)! ) * p^k * (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:
μ = k₁ * p₁ + k₂ * p₂ + k₃ * p₃ + 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)! ) * p^k * (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 σ². 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:
σ² = (k₁ - μ)² * p₁ + (k₂ - μ)² * p₂ + (k₃ - μ)² * p₃ + 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: σ² = N*p*(1-p). Dla małych p będzie więc w przybliżeniu σ² = μ. To ciekawy i ważny wynik, przydatny w statystyce. Domyślamy się, wobec tego, że dla rozkładu Poissona dokładnie: σ² = μ = λ.
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-μ)²f(x).
Dla normalnego wychodzi (domyślamy się) σ² = σ².
Dla wykładniczego: σ² = T².

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ć.
σ² = 1/12, σ = 0.29. To jest ta dwunastka z sumy dwunastu rand'ów. Trzeba ich aż tyle, żeby osiągnąć σ². Bo wariancje (σ²) się sumują. Dlaczego?

3. Czy są jeszcze jakieś ogólne parametry rozkładów? Dlaczego nie?
σ² = Σ(k_i - μ)² * p_i może jakiś sens miałyby jeszcze potęgi 3 i 4? Σ(k_i-μ)³*p_i. 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 σ² = 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(-x²/2) to jeden wymiar (patrz),
f(y)=(2*pi)^-1/2*exp(-y²/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(-(x²+y²)/2) Poniżej dołączamy 4 należyte parametry:
f(x,y)=1/(2*pi*σ_x*σ_y)*exp(-(((x-μ_x)/σ_x)²+((y-μ_y)/σ_y)²)/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*σ_x*σ_y*(1-ρ²)^1/2)*exp(-1/(2*(1-ρ²))*(((x-μ_x)/σ_x)²+((y-μ_y)/σ_y)²)-2*ρ*(x-μ_x)/σ_x*(y-μ_y)/σ_y))
Jeżeli wyrzucić cztery parametry, które już znamy, będzie krócej.
f(x,y)=exp(-(x²+y²-2*ρ*x*y)/(2*(1-ρ²)))/(2*pi*(1-ρ²)^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 (ρ_ij=ρ_ji), 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.

P-wartość (p-value) ↑↑↑

Wykonujemy RAND(), wychodzi 0.324, 0.892, 0.277. A nagle widzimy 1.123!!! NIEMOŻLIWE. Oczywiście popatrzyliśmy nie na tę komórkę. Ale co jeżeli? Wpisujemy =RAND() (powiedzmy, pierwszy raz dzisiaj - to ważne) i naciskamy Enter. Widzimy 0.01326. Jak reagujemy? W ogóle nie reagujemy bo liczba jest między 0 i 1, czyli OK? Dziwi nas, że jakoś taka mała. Podobnie, zdziwiło by nas 0.99324, nie mówiąc już o 0.000156, czy 0.9999847. Tak, zdziwienie jest właściwą reakcją. Co prawda, prawdopodobieństwo trafienie w przedział (0.5, 0.51) jest takie samo jak w przedział (0, 0.01), ale ten drugi (podobnie jak (0.99, 1)) jest wyjątkowy, jest skrajny. Cała statystyka (patrz niżej) polega na wypatrywaniu, czy p-value wyszło nam mniejsze od 0.05 (albo większe od 0.95). A tu nagle, głupie RAND() daje nam taki wynik. Owszem, takie myślenie nie wygląda na poważne. Ale jest. Po prostu, jesteśmy na granicy rachunku prawdopodobieństwa i statystyki.
Weźmy rozkład normalny =NORMSINV(RAND()). Wpisujemy, Enter, a tu 3.456. Jakby dużo, możliwe, ale dużo. Możliwe jest też 10, 100 itd. (oczywiście nie w sensie funkcji arkusza, ta ma ograniczenia, teoretycznie MOŻLIWE). Jednak dostając wynik 5.678 od razu wiemy, że błąd. Gdzie jest granica? Ktoś, kto pamięta zasadę 68–95–99.7 wie jakie wartości są sensowne. Ale jak ktoś nie pamięta, albo jak to nie rozkład normalny tylko jakiś inny, albo w ogóle wynik jakiegoś skomplikowanego rozumowania statystycznego? DLATEGO nie wartość, tylko p-wartość się podaje. P - prawdopodobieństwo, jest między 0 i 1, i od razu wszystko wiadomo. Podobnie RAND, ale to, w jakimś sensie przypadek. Jak RAND pomnożymy przez 10cm i dodamy 122cm, to wynik nie będzie już przypominał żadnego prawdopodobieństwa. Natomiast mnożenie prawdopodobieństwa przez 10cm (tym bardziej przez 10 metrów to żart) nie ma żadnego sensu. Chyba, żeby głębiej pomyśleć... Prawdopodobieństwo czasem jest interpretowane jako część. Część kulek jest zielona, 30%. A "poproszę część jabłka, 50%" ma sens. Poniższy temat: statystyka, pełen jest takich zawiłości.
Jak przerobić "normalne" 3.456 na p-value? Bardzo prosto: =NORMSDIST(3.345). Jaki wynik? Taki: 0.9996, dużo, blisko 1.
Z p-value jest tak, że godne uwagi jest jak jest <0.05. Ten przypadek: =NORMSDIST(-3.345) by to spełniał. Sprawa dwóch ogonów rozkładu jest dość prosta. Wydaje się, że nie ma sensu dzielić włosa na czworo (na dwa). Trudniejsze jest coś innego. Tam mamy kwestię binarną, prawa strona, czy lewa strona (a może obie?). Generowaliśmy liczbę losową normalną (standardową) i wyszło nam 3.456. Tych cyferek było 15 (3.45678901234567), arkusz nie wszystkie nam pokazał (i tak za dużo pokazuje, chce się pokazać jako bardzo dokładny). Jakie jest (teoretycznie) prawdopodobieństwo wylosowanie dokładnej liczby z rozkładu ciągłego? Zero jest. Jakie jest prawdopodobieństwo wylosowania w arkuszu 3.45600000000000 albo 3.45678901234567? Około 1e-15. Ale nie to jest interesujące. Co jest interesujące? Jakie jest prawdopodobieństwo otrzymania 3.456 lub liczby jeszcze bardziej ekstremalnej (dziwnej), czyli P(z>3.456), czyli =1-NORMSDIST(3.345). Nie jest dziwne otrzymanie 1.1, -1.9, 0.55. Otrzymanie 3.345 jest dziwne i wymaga wyjaśnienia. Otrzymanie p-value<0.05 wymaga wyjaśnienia. Naukowcy się cieszą, bo piszą artykuł. Dlaczego akurat 0.05? To jest bardzo dyskusyjne, ale czegoś trzeba się trzymać, na początku.

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*Σx_i 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 x₁+x₂, ż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 σ², która się pojawiła na skutek sumowania (splatania).
(2*pi*σ²)^(-1/2)*e^(-x^2/(2*σ²)) - to jest N(0,σ), a nam w splocie ta σ wyszła równa 2^1/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). Wygodniej użyć wariancji σ², ż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 = μ₁ + μ₂ ? Z takiego samego rozumowania, jak to które doprowadziło do σ_suma² = σ₁² + σ₂² , 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(μ, σ²) ma rozkład N(n*μ, n*σ²).
Drobna uwaga. Jaki rozkład ma: y=a*x+b? Taki ma: N(a*μ+b, (a*σ)²). 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.17²+0.1²)^1/2=0.2m). Zwróćmy uwagę, że wariancja nie ma takiego wymiaru jak x, μ, czy σ, ma kwadrat tego wymiaru. Jeżeli σ=17cm, to σ²=290cm². Niby banał, a ile pomyłek. Pomijając już tę, że ktoś tu zarzuci, że 17²=289. Ale to jest inna kwestia, raczej trudna.

2e. Tak więc suma n zmiennych podzielona przez n, czyli średnia, ma rozkład N(μ, σ²/n). Czyli σ_śr=σ/√n . Powoli: σ² sumy = n*σ², 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, σ² zdefiniowany jest jako wartość oczekiwana kwadratu różnicy zmiennej losowej minus μ, czyli takiej wielkości: (x-μ)². To może estymatorem tego parametru byłaby średnia z (x-x_śr)²? 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)². Za jaką cenę (czy to jest jakaś strata?)? Wzór jest taki: s² = 1/(n-1)*Σ(x-μ)². 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 σ² jest zdefiniowane przed chwilą s². To dobrze, że używa się innej litery, to JEST coś innego! Natomiast s=√s² (że tak to napiszemy) nazywa się odchyleniem standardowym. Tu niestety drobna uwaga. s² jest "nieobciążonym" estymatorem σ², tzn. wartość oczekiwana estymatora równa się estymowanemu parametrowi. Niestety, nie dotyczy to s, które otrzymujemy z nieliniowego przekształcenia (pierwiastek) s². 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 - x₀) / s_śr. To x₀, 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 x₀ (μ=x₀), to x_śr - x₀ ma rozkład N(0,σ²/n). Jeżeli podzielilibyśmy (x_śr - x₀)/σ, 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 - x₀) / 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 σ² = 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 μ_A=μ_B? Nie interesują nas σ. Natomiast, jeżeli WIEMY, że σ_A=σ_B, 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 σ_A=σ_B, 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).

4. Przykład wg rozwiązania zadania 3 przez JF.

Dane są umieszczone w wierszach, a nie w kolumnach. Wynika to z tego, że są to tylko dwa lata (dla potrzeb zadania) z całego szeregu lat, chyba do 2016 (ciekawiej byłoby wziąć ostatnie lata, kiedy już byliśmy na świecie (w 1951 to nawet autora tego tekstu nie było)). Najlepiej byłoby wziąć lata 2016 i 1951, wtedy można by się spodziewać zaobserwowania ocieplenia. To, czego się spodziewamy a priori, czyli zanim zaczniemy analizę, jest bardzo ważne w statystyce.
Zasadniczą cechą tych danych jest to, że lipiec jest ciepły, a styczeń zimny. I nie należy tych miesięcy mieszać. Mamy tu dane sparowane, ostatni parametr formuły musi być 1.
Spróbujmy mimo wszystko zignorować wiedzę o lecie i zimie i policzmy najzwyklejsze średnie i odchylenia standardowe. To te brązowe liczby po prawej. W kolumnie P jest kolumna O podzielona przez √12, czyli jest to odchylenie standardowe średniej. Tu drobna uwaga. Powyższa uwaga jest nieprawdziwa. Byłaby ok, gdybyśmy mieli do czynienia z rozkładem normalnym, ale nie mamy. Kto nie wierzy, niech wrzuci powyższe dane do histogramu (-1.30 0.60 1.70 8.80 12.40 16.90 17.90 19.00 15.00 5.80 6.60 1.30 -0.30 -1.50 -3.50 10.10 11.60 15.60 18.50 18.90 11.80 7.60 1.90 -1.40 po przerobieniu ich na kolumnę). Pomijamy nieścisłości. Różnica średnich jest 1.3°C a niepewności pomiarowe ich obu są ponad 2. Wrzućmy do Wolframa nasze rozkłady:
exp(-((x-8.73)/2.12)^2/2)    i    exp(-((x-7.44)/2.35)^2)

To są rozkłady średnich. Czy one się różnią? Niby tak. Ale wyobraźmy sobie, że losujemy z tych rozkładów po jednej liczbie (PO JEDNEJ). Prawdopodobieństwo, że ta z prawego (niebieskiego) rozkładu będzie większa, owszem, jest większe od 0.5, może wynosi nawet 0.6 (ile dokładnie?).
Zupełnie fałszywe podejście do danych, ignorujące sezonowość, wskazuje, że lata 1951 i 1952 się nie różnią. Gdybyśmy nie znali sparowanego testu t, moglibyśmy jeszcze tak zrobić. Do komórki B4 wpisać formułę =B2>B3 i przeciągnąć w prawo. Otrzymujemy 4 fałsze i 8 prawd. Wyraźnie więcej prawd, ale czy istotnie więcej? Odpowiedź: =ROZKŁ.DWUM(4;12;0.5;PRAWDA) daje wartość 0.19 - nieistotnie. Nawet tylko 3 fałsze jeszcze są nieistotne (0.07>0.05). Test oparty o prostą kwalifikację większe/mniejsze jest za słaby.
Na szczęście mamy test t Studenta. Wynik testu sparowanego (do czego nie mamy żadnych wątpliwości) i, przy tym jednostronnego (do dalszej dyskusji) został pokolorowany na czerwono bo jest <0.05!
Dalsza dyskusja: Do testu jednostronnego mamy prawo gdy wiemy (posiadamy wiedzę), że, jeżeli już, to 1952 może być tylko zimniejszy. Jakby był cieplejszy, to nas to nie interesuje. Jeszcze raz. Pytamy testu (tego testa): Czy 1952 jest zimniejszy od 1951, czy nie jest. W takim wypadku czerwono. ALE przecież my nie mamy takiej wiedzy. Jeżeli byśmy już mieli, to o ociepleniu, a nie ochłodzeniu. Kropka. Pytamy więc, czy 1951 różni się od 1952, czy nie, w którąkolwiek stronę. Test dwustronny. P-value teraz jest dwa razy większe, już nie jest mniejsze od 0.05. Nie ma podstaw do utrzymywania hipotezy, że się różnią. Dlaczego liczba w G6 jest dwa razy większa od A6? Bo ogony rozkładu t są symetryczne i dwa są dwa razy cięższe niż jeden.
O niższych wierszach tego arkusza szkoda gadać. Może tylko zwróćmy uwagę na piękne słowo "homoscedastyczna", żeby go przypadkiem nie zapamiętać (już za późno?). Mieszanie lipca ze styczniem nie daje żadnych szans testowi, wyniki kompletnie nieistotne. P-value to, mniej więcej prawdopodobieństwo otrzymania otrzymanego wyniku przy założeniu braku różnicy między 1951 i 1952.
Pamiętajmy jednak na przyszłość, że sytuacja wymagająca testu sparowanego jest raczej rzadkością. Najczęściej mamy pięć zielonych jabłuszek i siedem czerwonych, i po zważeniu każdego pytamy, czy rzeczywiście czerwone są cięższe (jeden ogon)? A może pytamy: czy ciężar zależy od koloru (dwa ogony)?

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ą H₀. Była też alternatywna hipoteza, H₁: nie warto.
Mamy zmienną losową normalną N(μ,σ). A przynajmniej tak nam się wydaje. Skupmy się na μ. Wydaje nam się (jesteśmy przekonani), że μ=μ₀. To jest nasza hipoteza zerowa, H₀. Nasze pierwszaki mają μ=μ₀=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. Przeciwko hipotezie alternatywnej: μ₀≠0. To brzmi jakoś niezbyt mądrze. Co, jeżeli μ₀=0.0001? Nic specjalnego, test nie wykaże nam "prawdy". Przy małym odchyleniu od H₀ 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 H₀, bo dopiero prawdopodobieństwo wyjścia poza przedział (-2, 2) jest małe, wynosi 0.05. Czyli 2.1 kazałoby nam odrzucić H₀. 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ć H₀.

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ć H₀: 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 H₀, 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 H₀ 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:
H₀: N(0, 1), H₁: 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:
x_gr =NORMSINV(1-0.05) = 1.65. Popatrzmy na ten rysunek. Bladą kreską narysowany tam jest (przypadkowo) przypadek α=0.1 (x_gr=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 H₁=H₀, 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 χ² ↑↑↑

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 χ². 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 χ² 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, σ²=Np). Drugi parametr piszemy jako σ, albo jako σ², jak nam wygodnie. Byleby nie było nieporozumienia. Tak więc taka zmienna: z²=(k-Np)²/Np jest kwadratem zmiennej N(0,1). Dlaczego z², 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 z² zsumować, i to właśnie jest zmienną losową χ²
χ²_n 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 χ²_n. Tylko ile wynosi n? Otóż n=9, o jeden mniej niż jest składników w H4:H13. Dlaczego o jeden mniej. Bo definicja χ²_n 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 CHI².

To jest p-value, w H2, niebieskie. Czyli wyszedł nam wynik (16.40), którego prawdopodobieństwo, wobec H₀ (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 χ². Za duża niezgodność. Za dobra zgodność, wskazywana przez małe χ², jest, co najwyżej podejrzana. Że może ktoś manipulował danymi (miotła Ockhama) i przesadził.

3. Rozkład χ². 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 σ²=2n, to wiemy już o rozkładzie wszystko.

4. Przykład (wg MB). W tym pliku tekstowym są następujące liczby: 19354922 12815475 8675982 6381966 5733259 ... 3 3 3 2 2 2 2 2 1 1 1. Liczb tych jest niemal 30k i są to liczby ludności miast w USA. Doskonale nadają się one do sprawdzenie pewnej dziwnej reguły, która mówi, że prawdopodobieństwo, że losowe miasto (i wiele innych liczb o wielu rzędach wielkości) ma liczbę ludności zaczynającą się na cyfrę k wynosi P(k)=log₁₀(1+1/k). Czyli najczęściej jedynka. I to są, chyba jednak przez przypadek, dwa pierwsze przypadki, i ostatnie. Te ostatnie już nie przez przypadek. Lekką bezmyślnością jest chyba jednak włączanie jedno-osobowych miast do analizy, może od dwucyfrowych należałoby zacząć? Ale nie ma to znaczenia, te kilka sztuk przy dziesiątkach tysięcy. Potem dość częste powinny być 2, 3, a 9 - najmniej prawdopodobne. Taka Excelowa formuła =LEWY(A1;1) wydobywa pierwszy znak. Potem trzeba policzyć ile, w tych 30 tysiącach przypadków mamy jedynek, dwójek, ..., dziewiątek. Uwaga, w oryginalnych danych, gdyby ktoś tam zaglądał, są miasta o liczbie ludności 0, odrzucamy je. Test χ² będzie miał 8 stopni swobody. P-value wychodzi niemal idealne; 0.45. Jeżeli liczba miast to n, k=1,2,.. 9, a n_k to liczba wystąpień cyfry k, na początku, oraz p_k prawdopodobieństwo (policzone z wzoru z logarytmem), to χ²=Σ(n_k-n*p_k)²/(n*p_k)

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≈A1². 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=Σ(x_i-x_śr)*(y_i-y_śr)/(s_x*s_y)
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):
(x₁, y₁), (x₂, y₂), ..., (x_n, y_n) 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 x_i są dokładnie znane, a tylko y_i 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 r₀=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 s_a:
D3: =(H11/(9-2)/F11)^0.5 -to jest s_a. 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 uogólnienia, do dowolnego n. Głównie wiersz 11 trzeba relokować. Żeby policzyć s_b 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 ± s_a i niżej b ± s_b te czerwone ± dopisano na bitmapie. Jak widać, a i b zgadzają się z tym co policzył wykres. Natomiast s_a i s_b 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. Przykład. Dla odmiany w Excelu, wg JM.
Poza oczywistymi rzeczami, w arkuszu mamy następujące formuły:
C18 =(B3-A$18)*(C3-B$18)   D18 =(B3-A$18)^2    E18 =(C3-B$18)^2 w wierszu 27 są sumy od 18 do 25.
F18 =C27/D27   G18 =B18-F18*A18 jak widać, zgadza się z równaniem linii trendu. Teraz to, czego nie daje nam linia trendu. Niepewności otrzymanych liczb a i b, czyli prawdziwa statystyka, rzetelna wiedza.
Reszty, czyli odchylenia punktów od prostej (w kierunku pionowym), porównajmy je z wykresem. H18 =C3-G$18-F$18*B3 popatrzmy na sumę reszt w H27, wynosi zero. Ta suma musi być zerowa. Zauważmy jednak, że zerowanie się sumy reszt nie jest wystarczającym warunkiem regresji, ta prosta mogłaby być dowolnie nachylona, nawet odwrotnie, a przesuwając ja w górę lub w dół uzyskalibyśmy sumę reszt zerową. Jaki jest warunek regresji, patrz następny punkt. W komórce H27 figuruje 4E-14. Wykładnik -15, -16 od razu wskazuje nam, że mamy zero (IEEE 754), tu jest trochę więcej, ale też wiele operacji arytmetycznych złożyło się na H27. W J mamy kwadraty x'ów (Netflixa), i na dole ich sumę. WYNIK:
K18 =(I27/6/D27)^0.5   L18 =K18*(J27/8)^0.5  Mamy n=8 punktów danych; 6=n-2.
Wynik końcowy regresji: HBO = (0.36 ± 0.10) * Netflix + (103.0 ± 8.1)
I co z tego? Ile będzie miało HBO, jak Netfix będzie miał 200mln użytkowników? Dopóki nie mieliśmy niepewności, było łatwo. HBO = 0.36 * 200 + 103.0 = 175. Oczywiście "na pewno". Nie wiemy, jak obliczyć niepewność tego 175. Są na to metody, ale nie wszystko na raz. Zróbmy tak, w końcu arkusz wiele ułatwia. Policzmy 4 wersje ++, +-, -+, --. Wynik jest taki 146, 162, 186, 202. Proponuje się niniejszym podsumować to tak: 175 ± 20 mln. Proszę nie dopytywać dlaczego akurat ±20. Zasada pesymizmu w ocenie danych nakazywałaby może raczej ±25.
Mętne? To instalujemy darmowy program R, i w nim:
x=c(23.5,33.3,44.4,54.5,70.8,89.0,110.6,139.3)
y=c(93,114,127,138,131,134,142,147)
d=data.frame(y,x)
r=lm(d)
x1=data.frame(x=c(200))
predict(r, newdata = x1, interval = "confidence")
(lm - linear model, reszta jako tako się tłumaczy, c - concatenate - łączy liczby w wektor)
Wynik predykcji: 174.5 141.2 207.8, wartość i przedział ufności. Najpewniej będzie 174.5, a z prawdopodobieństwem 0.95 prawda będzie w przedziale od 141 do 208.

2c. 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 = y_i - ( a*x_i + b ), to są odchylenia w kierunku pionowym (jak widać z wzoru) punktów od prostej. ε₅ jest największe. Parametry prostej (a, b) są tak dobrane, żeby suma kwadratów tych odchyleń Σε_i² 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 i ε₂=10 wtedy:
ε₁²=1 a ε₂²=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=ax². 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=ax²).
y = b_y + a_y*sin(a_x*x - b_x)
Cztery parametry do estymowania, to dużo. Bo zachodzi pewna, bardzo niekorzystna sytuacja. Jeszcze raz, te a_x, b_x, a_y, b_y trzeba tak dobrać, żeby Σε_i² było możliwie najmniejsze. Zadanie jest jasne. To szukamy a_x. Wymyślamy jakieś wartości i sprawdzamy sumę kwadratów. Jak dobrze wymyślamy, to dość szybko zaczniemy zacieśniać granice, osaczać optymalne a_x. I dalej, bierzemy się za b_x. Tylko zaraz, poprzednio jakieś b_x 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 Σε_i² trzeba policzyć 100⁴=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 10⁴pozycji, 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ż a_x 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ć b_y? 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 (ε²). 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 x_i+1=4*x_i*(1-x_i) 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 = Σ x_i/n = x₁/9 + x₂/9 + ... +x₉/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 x₅ jest w środku, jest najważniejsze, a x₁ i x₉ ledwo, ledwo. Czyli jak? Jakoś tak: x_śr = Σ w_i*x_i, gdzie w_i 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 w₁ do w₂₀. Jak chcemy użyć, w pewnym momencie tylko od w₅ do w₂₀, 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ść s_a. Jeżeli |a|>3*s_a, 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.

R ↑↑↑

Arkusze kalkulacyjne nadają się bardzo dobrze do statystyki. Mają wiele możliwości i są wizualne. Jeżeli dodać do tego Visual Basic w Excelu czy Open Office, albo Java Script w Arkuszach Googla, to możemy zrobić wiele. Po co nam więc R?

Może lepsza byłaby Statistica, lub inny tego typu komercyjny program? Te programy są drogie i są komercyjne, zorientowane na użytkownika. Działają tak, żeby użytkownikowi się podobało, żeby przedłużył licencję. Nie, żeby lepiej rozumiał, co robi. R jest darmowy, jest oszczędny w słowach, matematycznie abstrakcyjny. Na początku może wydawać się trudny.
Zainstalujmy sobie R. Warto! Nie zajmie wiele pamięci. A jeżeli nas zajmie...
Zdefiniujmy zmienną x i wpiszmy do niej wartość. Możemy z nią potem coś zrobić.

Poruszanie się w konsoli różni się od tego, co robi się na dotykowych ekranach. Trudno. A może dobrze?
Wynik jest 10, zgadza się. Ale po co indeks [1]? Dla porządku numerujemy. Rzadko mamy tylko jedną liczbę.
Więcej liczb (wektor) wpisuje się tak: c(1,2,3,4,5,4,3,2,1) Enter. Oczywiście musimy nasz wektor zapisać, piszmy tak: x=c(1:100) Enter x Enter.
Tak wygląda to numerowanie dla porządku.
Rozbawić kogoś może wynik wpisania x*x Enter - kończymy z tym przypomnieniem. Zróbmy też tak: plot(x*x).ta kropka kończy zdanie, nie formułę
Idźmy do statystyki: x=rnorm(100), mamy 100 liczb normalnych unormowanych. Możemy je zobaczyć. Albo coś takiego zobaczyć: 0.1*x+10.5.
Statystyka wygląda tak: mean(x) oraz sd(x). Policzmy sobie też mean(rnorm(10000)), wynik nie powinien się różnić od zera dużo bardziej niż o 0.01. Dlaczego? A dlaczego to zrobiliśmy? Żeby zobaczyć różnicę z arkuszem, w którym musielibyśmy to dziesięć tysięcy liczb mieć (widzieć). No to jeszcze hist(x) albo boxplot(x). No to jeszcze korelacja: r=0.9 i dalej y=r*x+(1-r^2)*rnorm(100) i plot(x,y) oraz cor(x,y). Teraz zróbmy jednak tabelkę (frame): f=data.frame(y,x), zobaczmy ją sobie; f Enter. Potrzebna nam jest po to, żeby zrobić regresję (linear model): reg=lm(f) i summary(reg). Wynik (poniżej w wersji nieco skróconej) jest okazją do powtórki ze statystyki.

Residuals: Min 1Q Median 3Q Max -0.45877 -0.12909 0.00571 0.14770 0.53439 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -0.02794 0.01916 -1.458 0.148 y 1.04916 0.02503 41.922 <2e-16

Co to są reszty (Residuals)? Mamy model (lm) naszych danych (f). Ten model, to funkcja liniowa y1=1.049*x-0.028 Liczba cyfr jest tu taka, jaka powinna być! Zmienne zależna została tu nazwana y1, żeby można było tę formułę wpisać do konsoli, by potem napisać y-y1 i zobaczyć owe reszty. Można też wykreślić dane (y) z predykcją wg modelu (y1):
plot(x,y,ylim=c(-3,3)); par(new = TRUE); plot(x,y1,col='red',ylim=c(-3,3))
dbając o taką sama skalę dla obu wykresów na osi y (zakres zmiennej x jest taki sam).
Jesteśmy przy Residuals. One są ważne, te reszty nie powinny wykazywać żadnej regularności, zależności od x. Powinny mieć rozkład normalny o wartości oczekiwanej równej 0. Dlatego tyle na ich temat. Sens Min i Max jest jasny. Mediana, to mediana, wartość środkowa, jest blisko zera. 1Q i 3Q, to kwartyle, coś jak mediana (=2Q), tylko ćwiartki. Otóż 1Q to taka wartość, że 1/4 wartości jest mniejsza.
Współczynniki. Intercept to rzędna początkowa, czyli b w takim zapisie: y=ax+b. Litera "y" oznacza a. Oba te parametry mają swoje odchylenie standardowe. Podane jest też p-value, tu oznaczone precyzyjnie jako Pr(>|t|), gdzie t to oczywiście zmienna t-Studenta. Dlaczego p-value jest tak niskie? Tak istotnie różny od zera jest współczynnik nachylenia a? Bo korelacja w danych była wysoka. Powtórzmy wszystko dla r=0.2, w R powtarzanie jest łatwe.

Używane wyżej komendy:
x=5    x=c(1,2,3,4,5,4,3,2,1)    x=c(1:100)    plot(x*x)    x=rnorm(100)    x1=0.1*x+10.5    mean(x)    sd(x)    mean(rnorm(10000))    hist(x)    boxplot(x)
r=0.9; y=r*x+(1-r^2)*rnorm(100); plot(x,y)
f=data.frame(y,x); reg=lm(f); summary(reg)
y1=x*summary(reg)$coefficients[2,1]+summary(reg)$coefficients[1,1]    (pod ,2] mamy błędy, pod 4 p-value)
plot(x,y,ylim=c(-3,3)); par(new = TRUE); plot(x,y1,col='red',ylim=c(-3,3))