Rachunek prawdopodobieństwa i statystyka Home
Wykład 2 Wykład 3 Wykład 4 Wykład 5 Wykład 6 Wykład 7 Wykład 8
w1, 29 lut cz 08.00 s221
Zaczęliśmy od danych o temperaturze powietrza. (Martwi nas przecież globalne ocieplenie.)
Wzięliśmy dane, które kiedyś można było ściągnąć z Univ. of Dayton, ale zdaje się, że ich tam już nie ma.
Jest to zapis temperatur dobowych w Warszawie (Krakowa nigdy tam nie było) od 1 stycznia 1995 do roku 2020.
Gdyby ktoś chciał, to oryginalny plik txt jest tu, a same temperatury, już w °C i z interpolowanymi brakami danych są tu.
Na wykładzie dopasowaliśmy "na oko" do danych funkcję sinus (zrobiliśmy to w Excelu).
Sinus ten ma oczywiście 4 parametry, z których jeden (T = 356.25 dni) nie podlega dyskusji, a pozostałe 3 dopasowaliśmy, jako się rzekło "na oko".
Teraz trzeba wymyślić jakąś obiektywną metodę doboru parametrów. Tak, by dwa wykresy, danych i modelu, jak najlepiej do siebie pasowały.
Te wykresy są wykresami dwóch kolumn liczb, które powinny się jak najmniej (na ile to możliwe) różnić od siebie.
w2, 7 marca
Wspomniane wyżej dwie kolumny liczb są trochę widoczne niżej w kolumnach F i G. W F są pomiarowe dane, w G są wartości obliczone z modelu:
temp = tempav + amp * sin(( t - faza )*2*pi/365.25) + a * ( t - 4632 )
temp - temperatura modelowa; tempav - średnia (average) temperatura obliczona dla 24, pełnych lat; amp - amplituda sinusa [°C];
faza - przesuniecie fazowe [dni]; a - współczynnik kierunkowy prostej wzrostu temperatury [°C/dzień]; t - czas [dni];
liczba 4632[dni] to środek szeregu czasowego, odejmując go dostajemy prostą, która nie zmienia średniej temperatury, a jedynie daje jej przyrost lub spadek w czasie.
Obiektywną (w przeciwieństwie do porównywania "na oko" pomarańczowego modelu z granatowymi danymi) metodą oceny zgodności modelu z danymi
jest metoda najmniejszych kwadratów (Least Squares). Parametry tak trzeba dobrać, żeby suma kwadratów odchyleń danych od modelu
{Kolumna N= (F1 - G1) ^ 2 } była jak najmniejsza.
Poniżej makro (program) VisualBasic przepisujące wybrany parametr z M do H i następnie przepisujące sumę kwadratów z L1 do N Sub ls() 'Least Squares For i = 1 To 11 Cells(7, 8) = Cells(i, 13) Cells(i, 14) = Cells(1, 12) Next i End Sub
Paraboliczny wykres po lewej, na dole to sumy kwadratów w zależności od parametru a. Niestety minimum mamy przy a = 9.6 [°C/100 lat], co oznacza przyrost temperatury w Warszawie o 1°C co 10 lat. Trzeba to będzie jeszcze sprawdzić w programie STATISTICA, który będzie dopasowywał wszystkie parametry (3 lub 4) jednocześnie.
w3, 21 marca
Sprawdziliśmy w STATISTICA pesymistyczny (dla młodych mieszkańców Planety) wynik otrzymany w Excelu. Zastosowano metodę estymacji nieliniowej. Wzrost temperatury w Warszawie 0.96°C/10lat potwierdził się z dokładnością do dwóch cyfr znaczących (podobnie pozostałe trzy parametry).
Badaliśmy chaos deterministyczny; patrz np. tutaj.
Badaliśmy funkcję Excela LOS(), oraz runif(1000) w R. Przetworzyliśmy liczby losowe równomierne na rzut monetą albo kostką do gry.
Sprawdziliśmy, że średnia jest modelem danych minimalizującym sumę kwadratów odchyleń modelu od danych. Zobaczyliśmy też do czego może przydać się inny model (inny estymator wartości oczekiwanej - cokolwiek to znaczy), mianowicie mediana. Jest ona przydatna, gdy w danych występują błędy grube.
w4, 4 kwietnia
Ze względu na zamówienie z ćwiczeń omówiliśmy rozkład normalny.
Z danych o temperaturach warszawskich wzięliśmy różnice R = dane - model, które powinny mieć rozkład normalny.
STATISTICA zrobiła histogram i domyślnie dopasowała krzywą Gaussa.
Rozkład ma trochę za dużo dużych ujemnych wartości, jest lewostronnie skośny. To cecha pogody; częściej w zimie zdarzają się bardzo zimne dni, niż w lecie bardzo ciepłe.
Rozkład normalny ma dwa parametry: μ i σ, patrz tutaj.
Jakie w naszych danych mamy μ dowiadujemy się (w przybliżeniu) obliczając średnią. Parametr skali σ znajdujemy obliczjąc odchylenie standardowe,
patrz tutaj.
W R robiliśmy następujące rzeczy: > x=rnorm(1e6) > sum(x>3) [1] 1293 > 1-pnorm(3) [1] 0.001349898 Jak widać, bardzo ładnie się zgadza.
w5, 18 kwietnia
Był rozkład dwumianowy, czyli jakie jest prawdopodobieństwo otrzymania w N próbach
k sukcesów, jeżeli prawdopodobieństwo tego ostatniego w pojedynczej próbie (rzucie monetą) jest p.
Było porównianie rozkładu dwumianowego z normalnym, a więc oczywiście centralne twierdzenie graniczne rachunku prawdopodobieństwa
i kwestia zmiennych dyskretnych i ciągłych.
Jak widać jest to sprawa ważna, ilustrowana na www błędnym rysunkiem kolorowym
(na wykresie po prawej, punkty są po prostu połączone kreskami) i poprawnym czarno-białym.
w6, 25 kwietnia
Próbowaliśmy wykonać test statystyczny, czy przypadkiem ten wynik głosowania nie jest wynikiem rzucania monetą.
Okazało się, że nie. Prawdopodobieństwo (p-wartość) otrzymania aktualnego wyniku przy założeniu hipotezy zerowej H0: p=0.5 wychodzi niepoliczalnie małe.
W Excelu i w R. Mamy tu 10.6σ (σ = (N*p*(1-p))^0.5), więc p-value będzie rzędu 1e-20.
Zaczęliśmy też myśleć o relacji teoretycznego rozkładu do liczb losowych z niego generowanych:
> sum(rbinom(10000,10,0.5)==0)
[1] 9
> dbinom(0,10,0.5)*10000
[1] 9.765625
Jak widać, bardzo ładnie się zgadza. (Jaki rozkład ma liczba generowana przez pierwszą linijkę?)
w7, 16 maja
Temat: odchylenie standardowe średniej
> sd(replicate(10000,mean(rnorm(100))))
[1] 0.1000706
Obliczyliśmy tu odchylenie standardowe (sd) z serii 10000 średnich obliczonych, każda ze stu zmiennych normalnych N(0,1).
Wynik zgadza się z oczekiwaniami, jest zgodny z wzorem σśr = σ/√n, gdzie n=100, σ=1.
Tak więc, jeżeli mamy serię n liczb o rozkładzie normalnym i obliczyliśmy z nich średnią, to dokładność tej
średniej (σśr) obliczamy wg powyższego wzoru.
Natomiast do wyznaczenia liczby cyfr znaczących błędu średniej (σśr) używamy poniższej funkcji R
z drugim parametrem równym 2 (oczywiście nie musimy używać aż R, wystarczy rzecz zrozumieć).
> signif(123.4567, 4)
[1] 123.5
Następnie w średniej obcinamy tyle cyfr po przecinku,
żeby zostało ich tyle samo co w "błędzie" (nawet jeżeli nie ma żadnego przecinka (ani kropki dziesiętnej)).
Powinno wyjść jakoś tak:
139.8 ± 2.3 cm
w8, 23 maja
Dziś była korelacja (i zależność). Jako przykład wzięliśmy, jak zwykle dobowe temperatury warszawskie z 25 lat. Wzięliśmy 25 temperatur z 1 stycznia i tyleż samo z 2 stycznia. Jak się okazało współczynnik korelacji wynosił tylko 0.77. Oznacza to, że temperatura 2 stycznia danego roku nie musi być aż tak bardzo taka sama jak temperatura dnia poprzedniego. Jednak r=0.77 przy n=25 jest istotna statystycznie. Natomiast korelacja 1 stycznia z 7 stycznia wynosiła tylko 0.37 i była nieco poniżej istotności (na poziomie α=0.05). Korelacji 1 lipca z 1 stycznia nie było żadnej. Konkretnie r=-0.2, co nic nie znaczy, gdyż przy tak małej liczbie punktów (n=25) nie jest istotne statystycznie.
Robiliśmy wykresy rozrzutu, a nawet, w Statistice macierz korelacji, również w postaci macierzy wykresów.