Weblog Tomasza Przechlewskiego [Zdjęcie T. Przechlewskiego]


scrum
random image [Photo gallery]
Zestawienie tagów
1-wire | 18b20 | 1wire | 2140 | 3rz | adamowicz | alsamixer | amazon | anniversary | antypis | apache | api | applebaum | arm | armenia | astronomy | asus | atom.xml | awk | aws | bachotek | bakłażan | balcerowicz | balta | banan | bash | batumi | berlin | bibtex | bieszczady | biznes | blogger | blogging | blosxom | bme280 | bono | borne-sulinowo | breugel | bt747 | budapeszt | budyń | bursztyn | canon | cedewu | chello | chiller | chillerpl | chown | chujowetaśmy | ciasto | cmentarz | contour | coronavirus | covid19 | cron | css | csv | curl | cycling | d54250wykh | dbi | debian | dejavu | dhcp | dht22 | dia | docbook | dom | dp1500 | ds18b20 | dulkiewicz | dyndns | dynia | ebay | economy | ekonomia | elka | elm | emacs | emacs23 | english | ep | erasmus | erasmusplus | ess | eu | eurostat | excel | exif | exiftool | f11 | fc | fc11 | fc15 | fc29 | fc5 | fc8 | fedora | fedora21 | fenix | ffmpeg | finepix | firefox | flickr | folau | fontforge | fontspec | fonty | food | fop | foto | france | francja | fripp | froggit | fuczki | fuji | fuse | gammu | garmin | gawk | gazwyb | gdańsk | gdynia | gender | geo | geocoding | georgia | gft | git | github | gmail | gmaps | gnokii | gnus | google | googlecl | googleearth | googlemaps | gotowanie | gphoto | gphoto2 | gps | gpsbabel | gpsphoto | gpx | gpx-viewer | greasemonkey | gruzja | grzyby | haldaemon | handbrake | hhi | historia | history | hitler | holocaust | holokaust | hp1000se | hpmini | humour | iblue747 | ical | iiyama | ikea | imap | inkscape | inne | internet | j10i2 | javascript | jhead | k800i | kajak | kamera | karob | kleinertest | kml | kmobiletools | knuth | kociewie kołem | kod | kolibki | komorowski | konwersja | krutynia | kuchnia | kurski | latex | latex2rtf | latex3 | lcd | legend | lenny | lesund | lewactwo | lgbt-folly | liberation | linksys | linux | lisp | lisrel | litwa | lizbona | logika | ltr | lubowla | lwp | lwów | m2wś | malta | mapquest | mapsource | marchew | marimekko | marvell | math | mathjax | mazury | mbank | mediolan | mencoder | mevo | mh17 | michalak | michlmayr | microsoft | monitor | mp4box | mplayer | ms | msc | mssql | msw | mswindows | mtkbabel | museum | muzyka | mymaps | mysql | nafisa | nanopi | natbib | navin | nekrolog | neo | neopi | netbook | niemcy | niemieckie zbrodnie | nikon | nmea | nowazelandia | nuc | nxml | oauth | oauth2 | obituary | odessa | okular | olympus | ooffice | ooxml | opera | osm | otf | otftotfm | other | overclocking | ozbekiston | panoramio | paryż | pdf | pdfpages | pdftex | pdftk | pedophilia | perl | photo | photography | picasa | picasaweb | pim | pine | pis | pit | plotly | pls | plugin | po | podróże | pogoda | politics | polityka | polsat | portugalia | postęp | powerpoint | połtawa | prelink | problem | propaganda | pstoedit | putin | python | pywws | r | radio | random | raspberry | raspberry pi | raspberrypi | raspbian | refugees | relaxng | ridley | router | rower | rowery | rpi | rsync | rtf | ruby | rugby | rumunia | russia | rwc | rwc2007 | rwc2011 | rwc2019 | rzym | samba | sds011 | selenium | sem | sernik | sheevaplug | sienkiewicz | signature | sks | skype | skytraq | smoleńsk | sqlite | srtm | sshfs | ssl | staszek wawrykiewicz | statistics | stats | statystyka | stix | stretch | suwałki | svg | svn | swanetia | swornegacie | szwajcaria | słowacja | tbilisi | terrorism | tex | texgyre | texlive | thunderbird | tomato | totalnaopozycja | tourism | tramp | trang | transylwania | truetype | ttf | turcja | turkey | turystyka | tusk | tv | tv5monde | twitter | typetools | ubuntu | uchodźcy | udev | ue | ukraina | umap | unix | upc | updmap | ups | utf8 | uzbekistan | varia | video | vienna | virb edit | vostro | wammu | wdc | wdfs | weather | weathercloud | webcam | webdav | webscrapping | weewx | wh2080 | wiedeń | wikicommons | wilno | win10 | windows | windows8 | wine | wioślarstwo | word | wordpress | wrt54gl | ws1080 | wtyczka | wunderground | ww2 | www | wybory | wybory2015 | włochy | węgry | xemex | xetex | xft | xhtml | xine | xml | xmllint | xsd | xslt | xvidtune | youtube | yum | zakopane | zakupy | zdf | zdrowie | łeba | świdnica | żywność
Archiwum
O stronie
wykorzystywany jest blosxom plus następujące wtyczki: tagging, flatarchives, rss10, lastbuilddatexhtmlmime. Niektóre musiałem dopasować nieco do swoich potrzeb. Więcej o blosxom jest tutaj
Subskrypcja
RSS 1.0
Worldometer vs ECDC

Jak już pisałem danych nt COVID19 jest multum bo są traktowane jako treść promocyjna, przyciągająca klikających. Każda tuba medialna (gazeta/portal/telewizja) w szczególności publikuje dane nt.

Źródłem pierwotnym każdego wydają się być raporty narodowe (bo jak inaczej), ale ponieważ te raporty narodowe są składane w różny sposób, to ich połączenie w jedną bazę też różnie może wyglądać. Generalnie ci co takie bazy robią to albo przyznają się, że działają na hmmm niekonwencjonalnych źródłach (Twitter) albo nic nie piszą, skąd mają dane. Mają i już...

Wydaje się (chyba, że czegoś nie wiem), że ECDC, OWiD, CSSE oraz Worldometers (dalej WMs) są najpopularniejszymi źródłami danych nt COVID19 w przekroju międzynarodowym. (Nawiasem mówiąc: WHO nie publikuje danych -- publikuje raporty z danymi w formacie PDF. Wydobycie z nich danych jest nietrywialne i kosztowne, bo nie da się tego na 100% zautomatyzować. W rezultacie prawie nikt nie powołuje się na WHO jako źródło danych -- lekki szejm przyznajmy, bo niby ta organizacja jest od tego żeby m.in. zbierać i udostępniać informację n/t.) Taka drobna różnica na początek: ECDC, OWiD oraz CSSE to prawdziwe bazy: zarejestrowane z dzienną częstotliwością zgony, przypadki, testy i co tam jeszcze. OWiD kopiuje dane z ECDC, kiedyś kopiowało z WHO ale napisali że WHO zawierało liczne błędy i to ich skłoniło do korzystania z ECDC (0:2 dla WHO). WMs publikuje stan na, bazy jako takiej nie ma (przynajmniej publicznie albo nie potrafię jej odszukać na stronie). Można założyć że jak się ogląda stronę WMs z ,,notowaniami'' nt/ koronawirusa w dniu X o godzinie T to jest to stan na X:T. Nawiasem mówiąc tak jest wygodniej, ale jednocześnie komplikuje to sprawę w aspekcie: dzienna liczba przypadków chociażby z uwagi na różnice czasu (jak w PL kończy się dzień to na Fiji jest w połowie inny; inna sprawa, że wątpię żeby ktoś się tym przejmował). Niemniej WMs ma rubrykę "nowe przypadki", tyle że nie bardzo wiadomo co to znaczy...

No więc po tym przydługim wstępie do rzeczy: jak się mają dane z WMs względem ECDC? Jak wspomniałem, na stronie WMs nie ma bazy -- jest tabela z danymi ze stanem ,,na teraz''. ECDC z kolei publikuje bazę w postaci arkusza kalkulacyjnego. Ściągam dane codziennie. Ze strony WMs o 21:00 (koniec dnia, przynajmniej w PL) oraz o 23:00 ze strony ECDC. Dane te wyglądają jakoś tak (WMs, po konwersji HTML→CSV):

date;country;totalC;newC;totalD;newD;totalT
04040600;USA;277467;+306;7402;+10;830095

Stempel czasu jest ustalany w momencie pobrania danych. Na stronie WMs czas nie jest podany explicite (nie ma czegoś takiego jak np. dane aktualizowano o H:M). Czyli 04040600 to dane z 2020/04/04 z godziny 6:00.

Dane ECDC wyglądają jakoś tak:

date;id;country;newc;newd;totalc;totald
2020-04-04;US;United_States_of_America;32425;1104;277965;7157

NewC -- nowe przypadki (dzienne); NewD -- nowe zgodny; totalC -- przypadki łącznie; totalD -- zgony łącznie. Baza ECDC ma stempel czasu (dzień).

W przypadku PL wiem, że Ministerstwo Zdrowia (MinZ) publikuje dane generalnie o godzinie 10-coś-tam oraz o 17/18-coś-tam. (Czemu tak nie wiem). Patrząc na dane z WMs wiedzę, że o 21:00 publikują już dane aktualne na ten dzień, w tym sensie, że uwzględnią stan z ostatniego dziennego komunikatu MinZ (ale jakiego formalnie dnia te dane dotyczą, to już inna sprawa, bo ten dzień przecież się nie skończył :-)). Jeżeli chodzi o ECDC to dane pobrane w dniu X zawierają dane do dnia X-1, żeby było śmieszniej ECDC dla tego dnia przypisuje dane z komunikatu MinZ z dnia X-2. Czyli na przykładzie: arkusz pobrany o 23:00 dnia 24/04/2020 będzie miał ostatni wiersz datowany 23/04 ale dane w tym wierszu będą tymi które pojawiły się na stronie MinZ w dniu 22/04.

Uzbrojony o tę wiedzę dla wybranych 24 krajów wykreśliłem dane (z kwietnia) w wersji WMs oraz ECDC, w dwóch wariantach: z oryginalnymi stemplami czasowymi (górny wiersz) oraz ze stemplem skorygowanym przy założeniu że dane ECDC są 24H opóźnione (czyli dzień 23/04 tak naprawdę to dzień 22/04 itd). Te ,,skorygowane dane'' to dolny wiersz. Dla 90% krajów dane łącznie nakładają się czyli dane są identyczne (są wyjątki--ciekawe czemu). Dane dzienne to misz-masz, każda baza ma własną wersję, nie wiadomo która jest prawdziwa, tyle, że ECDC ma zawsze dane dzienne a WMs niekoniecznie (dla Japonii prawie zawsze ta kolumna była pusta)

Dane i komplet wykresów jest tutaj

Poniżej kilka wybranych krajów:





url | Sun, 26/04/2020 05:17 | tagi: , ,
Wględne tempo wzrostu (koronowirusa)

Financial Times zamieścił wykres wględnego tempa wzrostu (rate of growth) czyli procentu liczonego jako liczba-nowych / liczba-ogółem-z-okresu-poprzedniego x 100%. Na wykresie wględnego tempa wzrostu zachorowań na COVID19 wszystkim spada: Every day the Covid-19 virus is infecting an increasing number of people, but the rate of growth in cases in some of the worst-hit countries is starting to slow. Powyższe Czerscy przetłumaczyli jako m.in. trend dotyczy niemal wszystkich krajów rozwiniętych. [he, he... Rozwiniętych pod względem liczby chorych, pewnie chcieli uściślić, ale się nie zmieściło]


Spróbowałem narysować taki wykres samodzielnie:

library("dplyr")
library("ggplot2")
library("ggpubr")
##
surl <- "https://www.ecdc.europa.eu/en/publications-data/download-todays-data-geographic-distribution-covid-19-cases-worldwide"
today <- Sys.Date()
tt<- format(today, "%d/%m/%Y")

#d <- read.csv("covid19_C.csv", sep = ';',  header=T, na.string="NA", stringsAsFactors=FALSE);
d <- read.csv("covid19_C.csv", sep = ';',  header=T, na.string="NA", 
   colClasses = c('factor', 'factor', 'factor', 'character', 'character', 'numeric', 'numeric'));

d$newc <- as.numeric(d$newc)
d$newd <- as.numeric(d$newd)

Zwykłe read_csv skutkowało tym, że newc/newd nie były liczbami całkowitymi, tylko czynnikami. Z kolei dodanie colClasses kończyło się błędem. W końcu stanęło na tym, że czytam dane w kolumnach newc/newd zadeklarowanych jako napisy a potem konwertuję na liczby. Czy to jest prawidłowa strategia to ja nie wiem...

Kolejny problem: kolumny newc/newd zawierają NA, wykorzystywana później funkcja cumsum z pakietu dplyr, obliczająca szereg kumulowany nie działa poprawnie jeżeli szereg zawiera NA. Zamieniam od razu NA na zero. Alternatywnie można korzystać z funkcji replace_na (pakiet dplyr):

# change NA to 0
d[is.na(d)] = 0

# Alternatywnie replace_na
#d %>% replace_na(list(newc = 0, newd=0)) %>%
#  mutate( cc = cumsum(newc), dd=cumsum(newd))

Ograniczam się tylko do danych dla wybranych krajów, nie starszych niż 16 luty 2020:

d <- d %>% filter(as.Date(date, format="%Y-%m-%d") > "2020-02-15") %>% as.data.frame
str(d)

last.obs <- last(d$date)
c1 <- c('IT', 'DE', 'ES', 'UK', 'FR')
d1 <- d %>% filter (id %in% c1) %>% as.data.frame

str(d1)

Obliczam wartości skumulowane (d zawiera już skumulowane wartości, ale obliczone Perlem tak nawiasem mówiąc):

t1 <- d1 %>% group_by(id) %>%  summarise(cc = sum(newc, na.rm=T), dd=sum(newd, na.rm=T))

t1c <- d %>% group_by(id) %>%  mutate(cum_cc = cumsum(newc), cum_dd = cumsum(newd)) %>% 
  filter (id %in% c1) %>%
  filter(as.Date(date, format="%Y-%m-%d") > "2020-02-15") %>% as.data.frame

  str(t1c)

Wykres wartości skumulowanych:

pc1c <- ggplot(t1c, aes(x= as.Date(date, format="%Y-%m-%d"), y=cum_cc)) + 
  geom_line(aes(group = id, color = id), size=.8) +
  xlab(label="") +
  theme(plot.subtitle=element_text(size=8, hjust=0, color="black")) +
  ggtitle(sprintf("COVID19: total confirmed cases (%s)", last.obs), subtitle=sprintf("%s", surl)) 

ggsave(plot=pc1c, "Covid19_1c.png", width=15)

Kolumny cum_lcc/cum_ldd zawierają wartości z kolumny cum_cc/cum_dd ale opóźnione o jeden okres (funkcja lag):

## 
t1c <- t1c %>% group_by(id) %>% mutate(cum_lcc = lag(cum_cc)) %>% as.data.frame
t1c <- t1c %>% group_by(id) %>% mutate(cum_ldd = lag(cum_dd)) %>% as.data.frame

t1c$gr_cc <- t1c$newc / (t1c$cum_lcc + 0.01) * 100
t1c$gr_dd <- t1c$newd / (t1c$cum_ldd + 0.01) * 100

## Początkowo wartości mogą być ogromne zatem
## zamień na NA jeżeli gr_cc/dd > 90
t1c$gr_cc[ (t1c$gr_cc > 90) ] <- NA
t1c$gr_dd[ (t1c$gr_dd > 90) ] <- NA

Wykres tempa wzrostu:

pc1c_gr <- ggplot(t1c, aes(x= as.Date(date, format="%Y-%m-%d"), y=gr_cc,  colour = id, group=id )) + 
  ##geom_line(aes(group = id, color = id), size=.8) +
  geom_smooth(method = "loess", se=FALSE) +
  xlab(label="") +
  theme(plot.subtitle=element_text(size=8, hjust=0, color="black")) +
  ggtitle(sprintf("COVID19: confirmed cases growth rate (smoothed)"), 
      subtitle=sprintf("%s", surl)) 

ggsave(plot=pc1c_gr, "Covid19_1g.png", width=15)




To samo co wyżej tylko dla PL/CZ/SK/HU:

c2 <- c('PL', 'CZ', 'SK', 'HU')

t2c <- d %>% group_by(id) %>%  mutate(cum_cc = cumsum(newc), cum_dd = cumsum(newd)) %>% 
  filter (id %in% c2) %>%
  filter(as.Date(date, format="%Y-%m-%d") > "2020-02-15") %>% as.data.frame

##str(t2c)
t2c.PL <- t2c %>% filter (id == "PL") %>% as.data.frame
t2c.PL
head(t2c.PL, n=200)

pc2c <- ggplot(t2c, aes(x= as.Date(date, format="%Y-%m-%d"), y=cum_cc)) + 
  geom_line(aes(group = id, color = id), size=.8) +
  xlab(label="") +
  theme(plot.subtitle=element_text(size=8, hjust=0, color="black")) +
  ggtitle(sprintf("COVID19: total confirmed cases (%s)", last.obs), subtitle=sprintf("Total: %s\n%s", lab1c, surl)) 

ggsave(plot=pc2c, "Covid19_2c.png", width=15)

t2c <- t2c %>% group_by(id) %>% mutate(cum_lcc = lag(cum_cc)) %>% as.data.frame
t2c <- t2c %>% group_by(id) %>% mutate(cum_ldd = lag(cum_dd)) %>% as.data.frame

t2c$gr_cc <- t2c$newc / (t2c$cum_lcc + 0.01) * 100
t2c$gr_dd <- t2c$newd / (t2c$cum_ldd + 0.01) * 100

## zamień na NA jeżeli gr_cc/dd > 90
t2c$gr_cc[ (t2c$gr_cc > 90) ] <- NA
t2c$gr_dd[ (t2c$gr_dd > 90) ] <- NA

t2c.PL <- t2c %>% filter (id == "PL") %>% as.data.frame
t2c.PL

pc2c_gr <- ggplot(t2c, aes(x= as.Date(date, format="%Y-%m-%d"), y=gr_cc,  colour = id, group=id )) + 
  ##geom_line(aes(group = id, color = id), size=.8) +
  geom_smooth(method = "loess", se=FALSE) +
  xlab(label="") +
  theme(plot.subtitle=element_text(size=8, hjust=0, color="black")) +
  ggtitle(sprintf("COVID19: confirmed cases growth rate (smoothed)"), 
      subtitle=sprintf("%s", surl)) 

ggsave(plot=pc2c_gr, "Covid19_2g.png", width=15)

Koniec

url | Tue, 31/03/2020 05:16 | tagi: , ,
Dane Eurostatu nt zgonów/urodzeń

Trzeba coś robić w czasie kwarantanny

## https://b-rodrigues.github.io/modern_R/
## https://gist.github.com/imartinezl/2dc230f33604d5fb729fa139535cd0b3
library("eurostat")
library("dplyr")
library("ggplot2")
library("ggpubr")
## 
options(scipen=1000000)
dformat <- "%Y-%m-%d"

eu28 <- c("AT", "BE", "BG", "HR", "CY", "CZ", "DK",
   "EE", "FI", "FR", "DE", "EL", "HU", "IE", 
   "IT", "LT", "LU", "LV", "MT", "NL", "PL", 
   "PT", "RO", "SK", "SI", "ES", "SE")
eu6 <- c("DE", "FR", "IT", "ES", "PL")

### Demo_mor/ Mortality monthly ### ### ###
dm <- get_eurostat(id="demo_mmonth", time_format = "num");
dm$date <- sprintf ("%s-%s-01", dm$time, substr(dm$month, 2, 3))
str(dm)

## There are 12 moths + TOTAL + UNKN
dm_month <- levels(dm$month)
dm_month

## Only new data
dm28  <- dm %>% filter (geo %in% eu28 & as.Date(date) > "1999-12-31")
str(dm28)
levels(dm28$geo) 

## Limit to DE/FR/IT/ES/PL:
dm6  <- dm28 %>% filter (geo %in% eu6)
str(dm6)
levels(dm6$geo) 

pd1 <- ggplot(dm6, aes(x= as.Date(date, format="%Y-%m-%d"), y=values)) + 
 geom_line(aes(group = geo, color = geo), size=.4) +
 xlab(label="") +
 ##scale_x_date(date_breaks = "3 months", date_labels = "%y%m") +
 scale_x_date(date_breaks = "6 months",
   date_labels = "%m\n%y", position="bottom") +
 theme(plot.subtitle=element_text(size=8, hjust=0, color="black")) +
 ggtitle("Deaths", subtitle="https://ec.europa.eu/eurostat/data/database (demo_mmonth)")

## Newer data
dm6  <- dm6 %>% filter (as.Date(date) > "2009-12-31")

pd2 <- ggplot(dm6, aes(x= as.Date(date, format="%Y-%m-%d"), y=values)) + 
 geom_line(aes(group = geo, color = geo), size=.4) +
 xlab(label="") +
 scale_x_date(date_breaks = "3 months", date_labels = "%m\n%y", position="bottom") +
 theme(plot.subtitle=element_text(size=8, hjust=0, color="black")) +
 ggtitle("Deaths", subtitle="https://ec.europa.eu/eurostat/data/database (demo_mmonth)")

ggsave(plot=pd1, file="mort_eu_L.png", width=12)
ggsave(plot=pd2, file="mort_eu_S.png", width=12)
## Live births (demo_fmonth) ### ### ###

dm <- get_eurostat(id="demo_fmonth", time_format = "num");
dm$date <- sprintf ("%s-%s-01", dm$time, substr(dm$month, 2, 3))
str(dm)

## There are 12 moths + TOTAL + UNKN
dm_month <- levels(dm$month)
dm_month

dm28  <- dm %>% filter (geo %in% eu28 & as.Date(date) > "1999-12-31")
str(dm28)
levels(dm28$geo) 

dm6  <- dm28 %>% filter (geo %in% eu6)
str(dm6)
levels(dm6$geo) 

pd1 <- ggplot(dm6, aes(x= as.Date(date, format="%Y-%m-%d"), y=values)) + 
 geom_line(aes(group = geo, color = geo), size=.4) +
 xlab(label="") +
 ##scale_x_date(date_breaks = "3 months", date_labels = "%y%m") +
 scale_x_date(date_breaks = "6 months", date_labels = "%m\n%y", position="bottom") +
 theme(plot.subtitle=element_text(size=8, hjust=0, color="black")) +
 ggtitle("Births", subtitle="https://ec.europa.eu/eurostat/data/database (demo_fmonth)")

##
dm6  <- dm6 %>% filter (as.Date(date) > "2009-12-31")

pd2 <- ggplot(dm6, aes(x= as.Date(date, format="%Y-%m-%d"), y=values)) + 
 geom_line(aes(group = geo, color = geo), size=.4) +
 xlab(label="") +
 scale_x_date(date_breaks = "3 months", date_labels = "%m\n%y", position="bottom") +
 theme(plot.subtitle=element_text(size=8, hjust=0, color="black")) +
 ggtitle("Births", subtitle="https://ec.europa.eu/eurostat/data/database (demo_fmonth)")

ggsave(plot=pd1, file="birt_eu_L.png", width=12)
ggsave(plot=pd2, file="birt_eu_S.png", width=12)
## Population (only) yearly ### ### ##
## Population change - Demographic balance and crude rates at national level (demo_gind)
dp <- get_eurostat(id="demo_gind", time_format = "num");
dp$date <- sprintf ("%s-01-01", dp$time)
str(dp)
dp_indic_dic <-  get_eurostat_dic("indic_de")

dp_indic_dic
dp28  <- dp %>% filter (geo %in% eu28 & time > 1999 & indic_de == "JAN")

str(dp28)
dp6  <- dp28 %>% filter (geo %in% eu6)

pdp1 <- ggplot(dp6, aes(x= as.Date(date, format="%Y-%m-%d"), y=values)) + 
        geom_line(aes(group = geo, color = geo), size=.4) +
        xlab(label="") +
        ##scale_x_date(date_breaks = "3 months", date_labels = "%y%m") +
        ##scale_x_date(date_breaks = "6 months", date_labels = "%m\n%y", position="bottom") +
        theme(plot.subtitle=element_text(size=8, hjust=0, color="black")) +
        ggtitle("Population", subtitle="https://ec.europa.eu/eurostat/data/database (demo_fmonth)")

ggsave(plot=pdp1, file="pdp1", width=12)

url | Wed, 25/03/2020 08:51 | tagi: , , ,
Dzienne dane dot. wypadków drogowych w Polsce

Na stronie http://policja.pl/pol/form/1,Informacja-dzienna.html udostępniane są dzienne dane dotyczące liczby interwencji, zatrzymanych na gorącym uczynku, zatrzymanych poszukiwanych, pijanych kierujących, wypadków, zabitych w wypadkach, rannych w wypadkach.

Ściągam wszystkie dane:

#!/bin/bash

rm pp.html

for ((i=0; i <= 274; i++)) do 
  if [ ! -f ${i}.html ] ; then
    curl -o ${i}.html "http://policja.pl/pol/form/1,Informacja-dzienna.html?page=${i}" ; 
    grep 'data-label' ${i}.html >> pp.html
    sleep 6
  else 
    grep 'data-label' ${i}.html >> pp.html
    echo done
  fi

done

zamieniam prostymi skryptami na plik CSV, który ma następującą strukturę:

data;interwencje;zng;zp;znk;wypadki;zabici;ranni
2008-12-01;NA;873;344;447;135;1;1

okazuje się że liczba interwencji jest podawana od roku 2018, wcześniej nie była. Nic to wstawiamy NA.

Na przyszłość dane będą aktualizowane w ten sposób, że codziennie (przez odpowiedni wpis w pliku crontab) będzie pobierany plik http://policja.pl/pol/form/1,Informacja-dzienna.html:

#!/usr/bin/perl
use LWP::Simple;

$PP="http://policja.pl/pol/form/1,Informacja-dzienna.html";
$PPBase="pp.csv";

$content = get("$PP");

$content =~ s/\r//g; # dla pewności usuń

@content = split (/\n/, $content);

foreach (@content) { chomp();
  unless ($_ =~ m/data-label=/ ) { next }

  if ($_ =~ m/Data statystyki/ ) { $d = clean($_); }
  elsif ($_ =~ m/Interwencje/ )  { $i = clean($_); }
  elsif ($_ =~ m/Zatrzymani na g/ ) { $zg = clean($_); }
  elsif ($_ =~ m/Zatrzymani p/ ) { $zp = clean($_); }
  elsif ($_ =~ m/Zatrzymani n/ ) { $zn = clean($_); }
  elsif ($_ =~ m/Wypadki d/ ) { $w = clean($_);  }
  elsif ($_ =~ m/Zabici/ )  { $z = clean($_);  }
  elsif ($_ =~ m/Ranni/ ) { $r = clean($_);
    $l = "$d;$i;$zg;$zp;$zn;$w;$z;$r";
    $last_line = "$l"; $last_date = "$d";
    ## pierwszy wpis powinien zawierać dane dotyczące kolejnego dnia
    ## więc po pobraniu pierwszego można zakończyć
    last;
 }
}

### read the database
open (PP, "<$PPBase") || die "cannot open $PPBase/r!\n" ;

while (<PP>) { chomp(); $line = $_; @tmp = split /;/, $line; }

close(PP);

### append the database (if new record)
open (PP, ">>$PPBase") || die "cannot open $PPBase/w!\n" ;

unless ("$tmp[0]" eq "$last_date") { print PP "$last_line\n" }
else {print STDERR "nic nowego nie widzę!\n"}

close(PP);

sub clean  {
 my $s = shift;
 $s =~ s/<[^<>]*>//g;
 $s =~ s/[ \t]//g;

 return ($s);
}

Zaktualizowana baza jest wysyłana na githuba. Tutaj jest: https://github.com/hrpunio/Nafisa/tree/master/PP

Agregacja do danych tygodniowych okazała się nietrywialna

Niektóra lata zaczynają się od tygodnia numer 0 a inne od 1. Okazuje się, że tak ma być (https://en.wikipedia.org/wiki/ISO_week_date#First_week):

If 1 January is on a Monday, Tuesday, Wednesday or Thursday, it is in W01. If it is on a Friday, it is part of W53 of the previous year. If it is on a Saturday, it is part of the last week of the previous year which is numbered W52 in a common year and W53 in a leap year. If it is on a Sunday, it is part of W52 of the previous year.

Nie bawię się w subtelności tylko tygodnie o numerze zero dodaję do tygodnia z poprzedniego roku.

Sprawdzam czy jest OK i się okazuje że niektóre tygodnie mają 8 dni. W plikach html są błędy:

Błędne daty 2019-10-30 winno być 2019-09-30; podobnie błędne 2019-03-28 (winno być 2019-02-28), 2018-11-01 (2018-12-01), 2018-12-01 (2017-12-01), 2016-04-30 (2016-03-30), 2009-08-31 (2009-07-31). Powtórzone daty: 2016-03-10, 2010-07-25, 2010-01-10 (zdublowane/różne/arbitralnie usuwamy drugi) Ponadto brak danych z następujących dni: 2015-12-04--2015-12-07, 2015-04-17--2015-04-20, 2014-10-02--2014-10-05, 2014-01-23 i jeszcze paru innych (nie chcialo mi się poprawiać starych.)

Teraz jest OK, plik ppw.csv ma nast strukturę:

rok;nrt;interwencje;in;zng;zngn;zp;zpn;znk;znkn;wypadki;wn;zabici;zn;ranni;rn;d1;d7 coś co się kończy na `n' to liczba tego co jest w kolumnie poprzedniej, np zn to liczba dni tygodnia dla kolumny zabici. Generalnie kolumny kończące się na `n' zawierają 7 :-) Kolumna d1 to pierwszy dzień tygodnia a kolumna d7 ostatni.

maxY <- max (d$zabici)
pz <- ggplot(d, aes(x= as.factor(nrt), y=zabici )) + 
 geom_bar(fill="steelblue", stat="identity")  +
 xlab(label="") +
 ggtitle("Wypadki/zabici (Polska/2020)", subtitle="policja.pl/pol/form/1,Informacja-dzienna.html") 

W sumie agregacja jest niepotrzebna, bo można ją zrobić na poziomie R używając funkcji stat_summary:

pw <- ggplot(d, aes(x= week, y=wypadki)) + 
 stat_summary(fun.y = sum, geom="bar", fill="steelblue") +
 scale_x_date( labels = date_format("%y/%m"), breaks = "2 months") +
 xlab(label="") +
 #theme(plot.subtitle=element_text(size=8, hjust=0, color="black")) +
 ggtitle("Wypadki (Polska/2018--2020)", subtitle="policja.pl/pol/form/1,Informacja-dzienna.html") 

albo najpierw agregując dane a potem wykreślając wykres szeregu zagregowanego. Drugi sposób pozwala na przykład na dodanie linii oznaczających poziomy zagregowanego zjawiska/etykiety słupków w sposób `inteligentny'. Dodajemy etykiety (z numerem tygodnia) tylko dla słupków poniżej/powyżej Q1/Q3:

## agregowanie do danych tygodniowych kolumn ranni, zabici, wypadki
dw <- d %>% group_by ( YrWeek) %>% summarise_at ( vars(ranni,zabici,wypadki), sum )

## Obliczanie mediany i kwartyli
median.zw <- median(dw$zabici)
q1.zw <- quantile(dw$zabici, 1/4) 
q3.zw <- quantile(dw$zabici, 3/4) 

## YrWeekZZ to numer tygodnia, dla tygodni w których liczba wypadków jest typowa
## numer jest pusty (żeby nie był drukowany); taki trick może jest lepszy
dw$YrWeekZZ <- substr(dw$YrWeek,4,5)
dw$YrWeekZZ[ (dw$zabici > q1.zw) & (dw$zabici < q3.zw) ] <- ""

pz.2 <- ggplot(dw, aes(x= YrWeek, y=zabici)) + 
 geom_bar(stat="identity", fill = "steelblue") +
 geom_text(data=dw, aes(label=sprintf("%s", YrWeekZZ), x=YrWeek, y= zabici), vjust=-0.9, size=3 ) +
 geom_hline(yintercept=median.zw, linetype="solid", color = "violet", size=.4) +
 geom_hline(yintercept=q3.zw, linetype="solid", color = "red", size=.4) +
 geom_hline(yintercept=q1.zw, linetype="solid", color = "red", size=.4) +
 xlab(label="rok/tydzień") +
 ylab(label="zabici") +
 scale_x_discrete(breaks=c("18/01", "18/10", "18/20",  "18/30", "18/40",
          "19/01", "19/10", "19/20",  "19/30", "19/40", "20/01", "20/10"))  +
          #  labels = c("/18/01", "18/10", "18/20", "")) ## tutaj niepotrzebne
 ggtitle("Wypadki/zabici (Polska/2018--2020)", 
  subtitle="Linie poziomie: q1/me/q3 (źródło: policja.pl/pol/form/1,Informacja-dzienna.html)") 

url | Wed, 25/03/2020 07:29 | tagi: , ,
Mevo: koncentracja rowerów na stacjach

29 czerwca był upał i pierwszy raz w życiu zobaczyłem stację MEVO literalnie zawaloną rowerami. Była 13:00. Po dojechaniu do domu sprawdziłem, że rowerów było tam aż 24. To mnie zainspirowało do sprawdzenia jak wygląda koncentracja rowerów na stacjach. Na szybko zmajstrowałem skrypt wypisujący ile jest rowerów na stacjach (rowery), udział w całości zaparkowanych w danym momencie (udzial) oraz udział w całości zaparkowanych w danym momencie w tym konkretnym mieście (ludzial): Dla 29 czerwca 2019, godz 13:00 okazało się że:

stacja;wspolrzedne;miasto;rowery;udzial;ludzial
S11358;18.57196808,54.40258423;Gdańsk;55;6.782;16.566
S11069;18.59038998,54.42794681;Gdańsk;24;2.959;7.229
S10100;18.57092997,54.44534256;Sopot;24;2.959;35.294
S11357;18.63640600,54.38637008;Gdańsk;18;2.219;5.422
S12007;18.53351997,54.49663118;Gdynia;16;1.973;7.843
S11186;18.57639103,54.39859534;Gdańsk;15;1.850;4.518
S10121;18.56255831,54.45392821;Sopot;13;1.603;19.118
S12126;18.47267507,54.54728035;Gdynia;13;1.603;6.373
S12053;18.54832346,54.51633152;Gdynia;13;1.603;6.373
S12054;18.54687652,54.51960321;Gdynia;12;1.480;5.882
S12033;18.56357898,54.48005340;Gdynia;10;1.233;4.902

Czyli na stacji 11358 było 55 rowerów co stanowiło 6,782% wszystkich zaparkowanych w systemie MEVO albo 16,566% zaparkowanych w Gdańsku. Dla Sopotu było nawet jeszcze lepiej bo na stacjach 10100/10121 było 24+13 (37 rowerów) ale było to 35,294 + 19.118, tj. prawie 55% wszystkich zaparkowanych w Sopocie (o godzinie 13:00). Dokładnie było 68 rowerów w 28 miejscach wtedy (27 stacji i jeden luźny bajk). Na 11 stacjach nie było nic a na 10 jeden rower.

Jest taka prosta miara koncentracji, co się nazywa w języku HH-Index, albo po polsku Wskaźnik Herfindahla-Hirschmana. Jest on nawet stosowany w USA do mierzenia koncentracji na rynku. Formuła jest banalnie prosta: dla $N$ wartości $x_i$ ($i=1...N$) sumujących się do 100 (czyli udziałów w całości), HHI liczy się jako: $\sum_{i=1}^N x_i^2$. Łatwo sprawdzić że HHI < 10000. Interpretacja jest taka, że HHI < 1000 wskazuje na słabą koncentrację, 1000 < HHI < 1800 umiarkowaną, a wartości większe od 1800 na dużą. BTW HHI da Sopotu o 13:00 (29/7/2019) wynosiło około 1850...

No to ja policzyłem HHI dla MEVO. Udziały były chwilowe, tj. $x_i = r_i/r_t$, dzie $r_t$ łączna liczba zaparkowanych rowerów na wszystkich stacjach w mieście $M$ (w danym momencie); no a $x_i$ to oczywiście liczba rowerów na stacji $i$. Potem uśredniłem, tj. wszystkie $HHI_i$ z godziny $h$ zsumowałem i podzieliłem przez liczbę pomiarów w tej godzinie (zwykle przez 30, bo pomiar jest co 2 minuty). Policzyłem oddzielnie dla GD/GA/Sopot/Tczew dla dni pracujących oraz dla świąt, sobót i niedziel osobno...

Wyniki dla maj--lipiec na wykresie (obok). Ciekawostkowo koncentracja w GD jest zaskakująco inna niż w GA. Teoretycznie im więcej stacji tym wartość HHI powinna być mniejsza, a tak nie jest: w GD wartość HHI wynosi w szczycie około 500, a w GA tylko 300. Szczyt wypada tak 9--11 zresztą. Można sobie wyobrazić, że użytkownicy jadą do pracy i zostawiając rowery przez biurami ogołacają stacje poza centrum a zapełniają te w rodzaju stacji o numerze 11358 (Gdańsk/Oliva Business Centre). W niedziele i święta do pracy nie jeżdżą to koncentracja jest mniejsza. Ma sens, ale nie w GA gdzie akurat w święta jest większa, wprawdzie chwilowo (w znaczeniu, że szybko rośnie ale potem równie szybko spada), ale jednak. Jakby w GA masowo jeździli gdzieś koło południa, a potem wracali z powrotem prawie że od razu... Mniejsza liczba stacji/rowerów w GA niż GD może powodować że wartość HHI ,,łatwiej'' rośnie...

W Sopocie i Tczewie jest znacząco mniej rowerów niż w GD/GA więc nic dziwnego że wartość HHI jest też dużo większa. BTW przeciętne dzienne wartości HHI są następujące (GD/GA/Sopot/Tczew): 170/222/1255/1143 (pon-piątek) oraz 93/225/1169/1200 (niedziele-soboty-święta). W Sopocie (poniedziałek--piątek) zmiany HHI są jeszcze inne niż w GD/GA. Amplituda jest mniejsza, a jedyny wyraźny dół jest rano około 5 a nie w godzinach 21--5 jak na przykład w GD. W niedziele i święta jest tradycyjnie jeden szczyt koło południa, a oprócz tego około 19--20 oraz 2--3 rano.

Skrypty i plik CSV z danymi jest tradycyjnie w archiwum GitHub

Na koniec przypomnienie, że prezes obiecał na 18.08 +4000 rowerów w systemie. Na koniec lipca wykazywane jest póki co: +1500 (niektóre w Warszawie albo w Cedrach Wlk.). Z tej puli bajków (numerów bajków?) codziennie pojawia się w pliku locations.js +1350, jeździ zaś mniej, około 1200 (reszta stoi). Na 100 procent wszyscy już wiedzą, że ni-chu-chu nie będzie 4000 rowerów nawet na 31 września, ale media zachowują w tej sprawie zgodne milczenie. Nikt nikogo nie pyta, a zwłaszcza prezesa. Nie ma sprawy...

url | Thu, 01/08/2019 12:38 | tagi: , , , ,
Wybory 2018. Różnica w liczbie mandatów do Sejmików

Powyborczo. Primo: szkoda, że PiSom te kamery nie wyszły byłoby się z czego pośmiać oglądając debili ze świeczkami (innego pożytku z zainstalowania -- przy założeniu #1kamera na jedną komisję -- 27tys kamer nie widzę).

Ale do rzeczy: dane pobrane z PKW (na Wikipedii za 2014 mają dokładnie takie same, za 2018 nie sprawdzałem)

require(ggplot2)

#d <- read.csv("mandaty.csv", sep = ';',  header=T, na.string="NA");
# Albo po prostu bo danych mało
# https://www.datamentor.io/r-programming/data-frame/
x <- data.frame("komitet" = c("PIS", "PO", "PSL", "SLD", "INNI"),
   "y2018" = c(254,194,70,15,19),
   "y2014" = c(171,179,157,28,20) )

# różnica w liczbie uzyskanych mandatów
d$diff <- d$y2018 - d$y2014 

ggplot(d, aes(x= komitet, y=diff, fill=komitet)) +
  geom_bar(stat="identity") +
  scale_fill_manual("legend",
    values = c("PIS" = "#421C52", "PO" = "blue",
    "PSL" = "green", "SLD" = "red", "INNI" = "pink")) +
    geom_text(aes(x=komitet, y=diff, label=diff),
    hjust=0, vjust=-0.25, size=3.5) +
ggtitle ("Mandaty sejmików wojewódzkich 2018--2014 (zmiana)")

BTW, jeżeli protokoły komisje obwodowe wysłały (zapewne elektronicznie) do PKW góra w poniedziałek (w mojej już poniedziałek-rano okleili kopią drzwi), to co niezawisłe Hermelińskie robiły we wtorek i środę? Niestety tego prostego pytania żaden z tzw. dziennikarzy (aka specjalistów od pierdołowatych njusów czyli #pierdokontentu) nie zadał.

A mnie ono ciekawi.

BTW2: ten wpis jest 500 w blogu. Wychodzi jakieś 45/rok średnio (z tendencją spadkową).

url | Fri, 26/10/2018 05:47 | tagi: , , ,
Wybory wójtów/burmistrzów/prezydentów

Analiza eksploracyjna wyborów wójtów/burmistrzów/prezydentów. W PL wybiera się radnych w wyborach do rad powiatów/rad gmin (oba ciała IMO zbędne), radnych sejmików wojewódzkich oraz uwaga: wójtów/burmistrzów/prezydentów (WBP) na poziomie gmin. O ile wybory sejmików kierują się tym samym mechanizmem co wybory sejmowe to wybory WBP są większościowe -- każdy może wystartować i wygrać. Do tego taki WBP ma dużą władzę więc warto być WBP. Takich wyborów w PL jest 2477 -- tyle ile gmin. W zależności od statusu gminy w jednych wybiera się wójta a w innych prezydenta czy burmistrza. Mówiąc konkretnie wójtów jest 1547, burmistrzów 823 a prezydentów 107. Poniższa tabela zestawia dane dotyczące kandydatów w wyborach 2018/2014/2010

  Rok      N   1KN   1K%   2KN   2K%   >4N    >4%     śr
  ------------------------------------------------------
  2018  6965   329  13,30  849  34,30  262  10,58   2,81
  2014  8019   245   9,90  666  26,90  471  19,00   3,25
  2010  7776   303  12,20  683  27.57  430  17,36   3,14

Jak na moje to kandydatów za dużo to się nie zgłasza, do tego (w tym roku/w tych wyborach) w 13,30% gmin jest jeden, a w 34,30% dwóch (co daje co najwyżej dwóch w prawie połowie wyborów WBP). Do tego tendencja jest jakby nie w tę stronę co trzeba: mniej kandydatów ogółem, więcej gmin z małą liczbą kandydatów, mniej gmin z dużą liczbą kandydatów. Można podsumować że demokracja na lokalnym poziomie słabnie...

Ilustruje to wykres krzywych gęstości liczby kandydatów na urząd WBP (dla każdego roku oddzielna krzywa).

## ramka g ma następującą strukturę: razem;teryt;rok
g$r <- as.factor(g$rok)
p <- ggplot(g, aes(x=razem, color=r)) + geom_density() +
labs(title="Krzywa gęstości liczby kandydatów na urząd wójta/burmistrza/prezydenta",
x="Liczba kandydatów", 
y = "Gęstość", color="Rok")

Dane są tutaj

url | Thu, 18/10/2018 17:31 | tagi: , , ,
Koniec pobierania danych wyborczych

Dobrnąłem w końcu do finału pobierając ostatecznie ze strony PKW dane dotyczące siedmiu wyborów, które odbyły się w latach: 2015, 2014 (samorządowe), 2011, 2010 (samorządowe), 2007, 2006 (samorządowe), 2005.

Wyniki wcześniejszych wyborów nie są już dostępne na poziomie komisji obwodowych (a przynajmniej ja nie potrafię takowych odszukać). Protokoły z wyborów z 2006 roku też nie były dostępne, ale udało się je w części odtworzyć ze stron z wynikami kandydatów (zawierającymi liczbę głosów oddanych na kandydata w poszczególnych komisjach obwodowych).

Dla każdych wyborów wykreśliłem histogram poparcia dla mainstreamowych partii: PSL, PO, PiS oraz SLD. Zgodnie z oczekiwaniami rozkłady poparcia są jednomodalne, prawostronnie symetryczne, ale z dwoma wyjątkami: rozkład poparcia dla PO jest bimodalny i ta tendencja wydaje się stała. Rozkład poparcia dla PSL z roku 2014 (cud nad urną) różni się -- na zasadzie znajdź element niepasujący do pozostałych -- od sześciu pozostałych rozkładów poparcia dla tej partii (czemu to już inna historia).

Dane są tutaj

url | Thu, 11/10/2018 08:31 | tagi: , , ,
Wybory 2014 i jeszcze więcej rozkładów

Rozkład odsetka głosów nieważnych (definiowanego jako głosy nieważne / (głosy ważne + nieważne)) w wyborach samorządowych w 2014. Pierwszy histogram dotyczy całej Polski (27455 komisji), drugi województwa pomorskiego (1856) a trzeci Mazowieckiego (3574).

#!/usr/bin/Rscript
# Skrypt wykreśla histogramy dla danych z pliku ws2014_komisje.csv
# (więcej: https://github.com/hrpunio/Data/tree/master/ws2014_pobranie_2018)
#
par(ps=6,cex=1,cex.axis=1,cex.lab=1,cex.main=1.2)
komisje <- read.csv("ws2014_komisje.csv", sep = ';',
       header=T, na.string="NA");

komisje$ogn <- komisje$glosyNiewazne  / (komisje$glosy + komisje$glosyNiewazne) * 100;

summary(komisje$glosyNiewazne); fivenum(komisje$glosyNiewazne);
sX <- summary(komisje$ogn);
sF <- fivenum(komisje$ogn);
sV <- sd(komisje$ogn, na.rm=TRUE)
skewness <- 3 * (sX[["Mean"]] - sX[["Median"]])/sV

summary_label <- sprintf ("Śr = %.1f\nMe = %.1f\nq1 = %.1f\nq3 = %.1f\nW = %.2f", 
  sX[["Mean"]], sX[["Median"]], sX[["1st Qu."]], sX[["3rd Qu."]], skewness)

## ##
kpN <- seq(0, 100, by=2);
kpX <- c(0, 10,20,30,40,50,60,70,80,90, 100);
nn <- nrow(komisje)

h <- hist(komisje$ogn, breaks=kpN, freq=TRUE,
   col="orange", main=sprintf ("Rozkład odsetka głosów nieważnych\nPolska ogółem %i komisji", nn), 
   ylab="%", xlab="% nieważne", labels=F, xaxt='n' )
   axis(side=1, at=kpN, cex.axis=2, cex.lab=2)
   posX <- .5 * max(h$counts)
text(80, posX, summary_label, cex=1.4, adj=c(0,1))

## ##
komisje$woj <- substr(komisje$teryt, start=1, stop=2)

komisjeW <- subset (komisje, woj == "22"); ## pomorskie
nn <- nrow(komisjeW)
sX <- summary(komisjeW$ogn); sF <- fivenum(komisjeW$ogn);
sV <- sd(komisjeW$ogn, na.rm=TRUE)
skewness <- 3 * (sX[["Mean"]] - sX[["Median"]])/sV

summary_label <- sprintf ("Śr = %.1f\nMe = %.1f\nq1 = %.1f\nq3 = %.1f\nW = %.2f", 
  sX[["Mean"]], sX[["Median"]], sX[["1st Qu."]], sX[["3rd Qu."]], skewness)

h <- hist(komisjeW$ogn, breaks=kpN, freq=TRUE,
   col="orange", main=sprintf("Rozkład odsetka głosów nieważnych\nPomorskie %i komisji", nn), 
   ylab="%", xlab="% nieważne", labels=T, xaxt='n' )
   axis(side=1, at=kpX, cex.axis=2, cex.lab=2)
   posX <- .5 * max(h$counts)
text(80, posX, summary_label, cex=1.4, adj=c(0,1))

komisjeW <- subset (komisje, woj == "14"); ## mazowieckie
nn <- nrow(komisjeW)
sX <- summary(komisjeW$ogn); sF <- fivenum(komisjeW$ogn);
sV <- sd(komisjeW$ogn, na.rm=TRUE)
skewness <- 3 * (sX[["Mean"]] - sX[["Median"]])/sV

summary_label <- sprintf ("Śr = %.1f\nMe = %.1f\nq1 = %.1f\nq3 = %.1f\nW = %.2f", 
  sX[["Mean"]], sX[["Median"]], sX[["1st Qu."]], sX[["3rd Qu."]], skewness)

h <- hist(komisjeW$ogn, breaks=kpN, freq=TRUE,
   col="orange", main=sprintf("Rozkład odsetka głosów nieważnych\nMazowieckie %i komisji", nn), 
   ylab="%", xlab="% nieważne", labels=T, xaxt='n' )
   axis(side=1, at=kpX, cex.axis=2, cex.lab=2)
   posX <- .5 * max(h$counts)
text(80, posX, summary_label, cex=1.4, adj=c(0,1))

Wyniki są takie oto (indywidualne wykresy tutaj: #01 #02 #03):

Rozkłady odsetka poparcia dla PSL/PiS/PO w wyborach samorządowych w 2014 w całej Polsce, w miastach/poza miastami oraz w poszczególnych województwach. Poniższy skrypt generuje łącznie 60 wykresów słupkowych:

#!/usr/bin/Rscript
# Skrypt wykreślna różnego rodzaju histogramy dla danych z pliku ws2014_komitety_by_komisja_T.csv
# (więcej: https://github.com/hrpunio/Data/tree/master/ws2014_pobranie_2018)
#
showVotes <- function(df, x, co, region, N, minN) {
   ## showVotes = wykreśla histogram dla województwa (region)
   kN <- nrow(df)
   sX <- summary(df[[x]], na.rm=TRUE);
   sV <- sd(df[[x]], na.rm=TRUE)
   ## współczynnik skośności Pearsona
   skewness <- 3 * (sX[["Mean"]] - sX[["Median"]])/sV

   summary_label <- sprintf ("Śr = %.1f\nMe = %.1f\nq1 = %.1f\nq3 = %.1f\nS = %.1f\nW = %.2f", 
     sX[["Mean"]], sX[["Median"]],
     sX[["1st Qu."]], sX[["3rd Qu."]], sV, skewness)

   if (minN < 1) {
   t <- sprintf("Rozkład głosów na %s\n%s ogółem %d komisji", co, region, kN ) } 
   else { t <- sprintf("Rozkład głosów za %s\n%s ogółem %d komisji (N>%d)", co, region, kN, minN ) } 

   h <- hist(df[[x]], breaks=kpN, freq=TRUE, col="orange", main=t, 
     ylab="%", xlab="% poparcia", labels=F, xaxt='n' )
     axis(side=1, at=kpN, cex.axis=2, cex.lab=2)
   ## pozycja tekstu zawierającego statystyki opisowe
   posX <- .5 * max(h$counts)
   text(80, posX, summary_label, cex=1.4, adj=c(0,1))
}

## Wczytanie danych; obliczenie podst. statystyk:
komisje <- read.csv("ws2014_komitety_by_komisja_T.csv", 
   sep = ';', header=T, na.string="NA");

komisje$ogn <- komisje$glosyNiewazne  / (komisje$glosy 
   + komisje$glosyNiewazne) * 100;

summary(komisje$PSL); summary(komisje$PiS); summary(komisje$PO);
fivenum(komisje$PSLp); fivenum(komisje$PiSp); fivenum(komisje$POp);

## ## ###
par(ps=6,cex=1,cex.axis=1,cex.lab=1,cex.main=1.2)
kpN <- seq(0, 100, by=2);
kpX <- c(0, 10,20,30,40,50,60,70,80,90, 100);
kN <- nrow(komisje)
region <- "Polska"
minTurnout <- 0

## cała Polska:
showVotes(komisje, "PSLp", "PSL", region, kN, minTurnout);
showVotes(komisje, "PiSp", "PiS", region, kN, minTurnout);
showVotes(komisje, "POp",  "PO",  region, kN, minTurnout);

## Cała Polska (bez małych komisji):
## ( późniejszych analizach pomijane są małe komisje)
minTurnout <- 49
komisje <- subset (komisje, glosyLK > minTurnout); 
kN <- nrow(komisje)

showVotes(komisje, "PSLp", "PSL", region, kN, minTurnout);
showVotes(komisje, "PiSp", "PiS", region, kN, minTurnout);
showVotes(komisje, "POp",  "PO",  region, kN, minTurnout);

## Typ gminy U/R (U=gmina miejska ; R=inna niż miejska)
komisjeW <- subset (komisje, typ == "U"); 
kN <- nrow(komisjeW)
region <- "Polska/g.miejskie"
showVotes(komisjeW, "PSLp", "PSL", region, kN, minTurnout);
showVotes(komisjeW, "PiSp", "PiS", region, kN, minTurnout);
showVotes(komisjeW, "POp",  "PO",  region, kN, minTurnout);

komisjeW <- subset (komisje, typ == "R"); 
kN <- nrow(komisjeW)
region <- "Polska/g.niemiejskie"
showVotes(komisjeW, "PSLp", "PSL", region, kN, minTurnout);
showVotes(komisjeW, "PiSp", "PiS", region, kN, minTurnout);
showVotes(komisjeW, "POp",  "PO",  region, kN, minTurnout);

## woj = dwucyfrowy kod teryt województwa:
komisje$woj <- substr(komisje$teryt, start=1, stop=2)

cN <- c("dolnośląskie", "dolnośląskie", "kujawsko-pomorskie",
 "lubelskie", "lubuskie", "łódzkie", "małopolskie", "mazowieckie",
 "opolskie", "podkarpackie", "podlaskie", "pomorskie", "śląskie",
 "świętokrzyskie", "warmińsko-mazurskie", "wielkopolskie",
 "zachodniopomorskie");
cW <- c("02", "04", "06", "08", "10", "12", "14", "16", "18",
 "20", "22", "24", "26", "28", "30", "32");

## wszystkie województwa po kolei:
for (w in 1:16) {
  wojS <- cW[w]
  ###region <- cN[w];
  region <- sprintf ("%s (%s)", cN[w], wojS);

  komisjeW <- subset (komisje, woj == wojS); ##

  showVotes(komisjeW, "PSLp", "PSL", region, kN, minTurnout);
  showVotes(komisjeW, "PiSp", "PiS", region, kN, minTurnout);
  showVotes(komisjeW, "POp",  "PO",  region, kN, minTurnout);
}
## ## koniec

Dla całej Polski wyniki są następujące:

Indywidualne wykresy zaś tutaj: #01 #02 #03 #04 #05 #06 #07 #08 #09 #10 #11 #12 #13 #14 #15 #16 #17 #18 #19 #20 #21 #22 #23 #24 #25 #26 #27 #28 #29 #30 #31 #32 #33 #34 #35 #36 #37 #38 #39 #40 #41 #42 #43 #44 #45 #46 #47 #48 #49 #50 #51 #52 #53 #54 #55 #56 #57 #58 #59 #60):

url | Tue, 02/10/2018 17:08 | tagi: , , ,
Żuławy w Koło 2018

Się tak rozochociłem, że zapisałem się na Żuławy w Koło 2018 (99 PLN). W ramach przygotowań, że tak powiem taktycznych postanowiłem rozpoznać możliwości przeciwnika:-) Konkretnie ustalić jak jechali ci co się zapisali, a co już startowali w ŻwK w roku 2017 albo 2016. Zadanie zatem polega na odszukaniu na liście zgłoszeń tych co się zapisali na edycję 2018 i jednocześnie ukończyli ŻwK w latach 2016/2017. Oczywiście nie ręcznie, tylko automatem:

## poniższe ściąga plik z listą zapisanych
wget 'http://www.czasomierzyk.pl/zapisy2016/zulawywkolo/index.php?akcja=lista' -O ZwK2018.out

Plik HTML ma tak prostą strukturę, że jego zamiana (za pomocą wyrażeń regularnych) na CSV jest banalna. Jak już mam ten plik CSV, to porównuję go do połączonych wyników z lat 2017/2016 (też w formacie CSV). Skrypt mam co porównuje pliki CSV:

perl join_csvs.pl -fn1 ZwK201809190908.csv  -fs1 1,2 -fn2 ZwK16_17.csv -fs2 1,2

Porównuje pliki ZwK201809190908.csv oraz ZwK16_17.csv, w oparciu o (wspólną) wartości dla kolumn nr 1 oraz nr 2 (w tym przypadku są to kolumny zawierające nazwisko i imię). Innymi słowy fs1 c1,c2..., to klucz główny, a fs2 c1,c2, to klucz obcy. Skrypt wypisuje połączone wiersze odpowiadające tym wierszom dla, których klucz główny = klucz obcy. Na dziś (19 września) takich wierszy wypisał 55, (na 104 zgłoszenia na dystansie 140km), ale pomijam tych co startowali kiedyś na najkrótszym dystansie lub tych, którzy startowali wprawdzie na najkrótszym, ale mieli średnią mniejszą niż 24kmh (odpada w ten sposób 10 zostaje 45). Na koniec plik jest zapodawany do prostego skryptu rysującego wykres słupkowy:

z <- read.csv("ZwK_2018_vs_2017.csv", sep = ';',
   header=T, na.string="NA", dec=".");

s140 <- summary(z$speed)

z <- subset (z, ( speed > 16.0 )); ## bez maruderów

# wykres słupkowy
h <- hist(z$speed,
  breaks=c(18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35),
  freq=TRUE,
  col="orange",
  main="Dystans: 140 (biorący udział w latach 2017-16)",
  xlab="Prędkość średnia w latach 2017--16 [kmh]",ylab="L.kolarzy",
  labels=T, xaxt='n' )
  axis(side=1, at=c(18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35))
  text(38, 37, summary_label, cex = .8, adj=c(1,1) )

Jak widać paru ludków w okolicach 30kmh jest. Będzie za kim jechać.

url | Wed, 19/09/2018 12:54 | tagi: , , ,
Rysowanie profilu wysokości w R

Ze śladu GPX prostym skryptem wyciągam co trzeba tworząc plik CSV o następującej zawartości (nazwy kolumn: data-czas,wysokość,prędkość,dystans przebyty):

daytime;ele;speed;dist  

Teraz poniższym skryptem rysuję profile wysokości (wysokość/prędkość vs czas oraz wysokość/prędkość vs dystans)

library(reshape)
require(ggplot2)

graphWd <- 6
graphHt <- 5

args = commandArgs(trailingOnly = TRUE);

if (length(args)==0) { stop("Podaj nazwę pliku CSV", call.=FALSE) }

fileBase <- gsub(".csv", "", args[1]);
outFile1 <- paste (fileBase, "_1.pdf", sep = "");
outFile2 <- paste (fileBase, "_2.pdf", sep = "");

what <- args[2];

# http://stackoverflow.com/questions/7381455/filtering-a-data-frame-by-values-in-a-column
d <- read.csv(args[1], sep = ';',  header=T, na.string="NA");
coeff <- median(d$ele)/median(d$speed)
d$speed <- d$speed * coeff


p1 <- ggplot(d, aes(x = as.POSIXct(daytime, format="%Y-%m-%dT%H:%M:%SZ"))) +
  geom_line(aes(y = ele, colour = 'wysokość', group = 1), size=1.5) +
  geom_line(aes(y = speed, colour = 'prędkość', group = 1), size=.5) +
  stat_smooth(aes(y=speed, x=as.POSIXct(daytime, format="%Y-%m-%dT%H:%M:%SZ"), colour ='prędkość wygładzona')) +
  ylab(label="Wysokość [mnpm]") +
  xlab(label="czas") +
  scale_y_continuous( sec.axis = sec_axis(name="Prędkość [kmh]",  ~./ coeff)) +
  labs(colour = paste( what )) +
  theme(legend.position="top") +
  theme(legend.text=element_text(size=12));
p1
ggsave(file=outFile1, width=graphWd, height=graphHt )

p2 <- ggplot(d, aes(x = dist)) +
  geom_line(aes(y = ele, colour = 'wysokość', group = 1), size=1.5) +
  geom_line(aes(y = speed, colour = 'prędkość', group = 1), size=.5) +
  ##geom_smooth() +
  stat_smooth(aes(y=speed, x=dist, colour ='prędkość wygładzona')) +
  ylab(label="Wysokość [mnpm]") +
  xlab(label="dystans") +
  scale_y_continuous( sec.axis = sec_axis(name="Prędkość [kmh]",  ~./ coeff)) +
  labs(colour = paste( what )) +
  theme(legend.position="top") +
  theme(legend.text=element_text(size=12));
p2

ps <- stat_smooth(aes(y=speed, x=dist));

ggsave(file=outFile2, width=graphWd, height=graphHt )

Teraz na koniec ciekawostka. Mój smartfon produkuje pliki GPX z superdokładnym stemplem czasu np. 2018-08-23T04:52:43.168Z, na czym wysypuje się R. Po prostu usuwam część po kropce dziesiętnej oraz samą kropkę (tj. .168Z) i działa.

url | Fri, 31/08/2018 07:40 | tagi: , ,
time plot tygodniowej liczby twitów

Załóżmy, że plik CSV zawiera liczbę opublikowanych twitów (dane tygodniowe). Problem: przedstawić szereg w postaci przebiegu czasowego (time plot). Taki skrypt R wymyśliłem do zrealizowania tego zadania:

require(ggplot2)

args <- commandArgs(TRUE)
ttname <- args[1];
file <- paste(ttname, ".csv", sep="")
filePDF <- paste(ttname, ".pdf", sep="")

d <- read.csv(file, sep = ';',  header=T, na.string="NA", );
## Plik CSV jest postaci:

##str(d)

## wiersze 1,2 + ostatni są nietypowe (usuwamy)
rows2remove <- c(1, 2, nrow(d));
d <- d[ -rows2remove, ];

## szacujemy prosty model trendu
lm <- lm(data=d, posts ~ no ); summary(lm)
posts.stats <- fivenum(d$posts);
posts.mean <- mean(d$posts);
sumCs <- summary(d$posts);

otherc <- coef(lm);
# W tytule średnia/mediana i równanie trendu
title <- sprintf ("Weekly for %s # me/av = %.1f/%.1f (y = %.2f x + %.1f)", 
  ttname, sumCs["Median"], sumCs["Mean"], otherc[2], otherc[1] );

##str(d$no)
## Oś x-ów jest czasowa
## Skróć yyyy-mm-dd do yy/mmdd
d$date <- sub("-", "/", d$date) ## zmienia tylko pierwszy rr-mm-dd
d$date <- sub("-", "", d$date) ## usuwa mm-dd
d$date <- gsub("^20", "", d$date) ## usuwa 20 z numeru roku 2018 -> 18
weeks <- length(d$no);
## https://stackoverflow.com/questions/5237557/extract-every-nth-element-of-a-vector
## Na skali pokaż do 20 element /dodaj ostatni `na pałę' (najwyżej zajdą na siebie)
## możnaby to zrobić bardziej inteligentnie ale nie mam czasu
scaleBreaks <- d$no[c(seq(1, weeks, 20), weeks)];
scaleLabs <- d$date[c(seq(1, weeks, 20), weeks)];

ggplot(d, aes(x = no, y = posts)) +
  geom_line() +
  ggtitle(title) +
  ylab(label="#") +
  xlab(label=sprintf("time (yy/mmdd) n=%d", weeks )) +
  scale_x_continuous(breaks=scaleBreaks, labels=scaleLabs) +
  geom_smooth(method = "lm")

ggsave(file=filePDF)  

url | Wed, 21/02/2018 13:56 | tagi: , ,
Żuławy wKoło 2017 -- podsumowanie



Podsumowanie wyników dla lat 2015--2017

z16 <- read.csv("wyniki_zulawy_2016_D.csv", sep = ';',  header=T, na.string="NA", dec=",");
aggregate (z16$time, list(Numer = z16$dist), summary)
z16$year <- 2016;

z15 <- read.csv("wyniki_zulawy_2015_D.csv", sep = ';',  header=T, na.string="NA", dec=",");
aggregate (z15$time, list(Numer = z15$dist), summary)
z15$year <- 2015;

z17 <- read.csv("wyniki_zulawy_2017_D.csv", sep = ';',  header=T, na.string="NA", dec=".");
aggregate (z17$time, list(Numer = z17$dist), summary)
z17$year <- 2017;

zz15 <- z15[, c("dist", "kmH", "time", "year")];
zz16 <- z16[, c("dist", "kmH", "time", "year")];
zz17  <- z17[, c("dist", "kmH", "time", "year")];

zz <- rbind (zz15, zz16, zz17);

## tylko dystans 140
zz140 <- subset (zz, ( dist == 140 ));
sum140 <- aggregate (zz140$kmH, list(Numer = zz140$year), summary)

boxplot (kmH ~ year, zz140, ylab = "Śr.prędkość [kmh]", col = "yellow", main="140km" )

## tylko dystans 55
zz75 <- subset (zz, ( dist > 60 & dist < 90 ));
sum75 <- aggregate (zz75$kmH, list(Numer = zz75$year), summary)
sum75
boxplot (kmH ~ year, zz75, ylab = "Śr.prędkość [kmh]", col = "yellow", main="80/75km" )

## tylko dystans 55
zz55 <- subset (zz, ( dist < 60 ));
sum55 <- aggregate (zz55$kmH, list(Numer = zz55$year), summary)
sum55
xl <- paste ("średnie 2015=", sum55$x[1,4], "kmh   2016=",
  sum55$x[2,4], "kmh   2017=", sum55$x[3,4], " kmh")

  boxplot (kmH ~ year, zz55, xlab = xl,
  ylab = "Śr.prędkość [kmh]", col = "yellow", main="55km" )

A ja (numer 418) byłem 70 w kategorii 140 km, z czasem 5:42:03 co dało 24,56 kmh przeciętną. Do pierwszego bufetu się spinałem, potem już nie...

url | Thu, 28/09/2017 05:17 | tagi: , ,
Chancellor Merkel victory for a visual person

Change in number of seats won by party (AfD is brown of course regardless official party colors :-):-)

library(ggplot2)

df <- read.csv("de.csv", sep = ';', header=T, na.string="NA");

ggplot(df, aes(x=party, y=diff, fill=party )) +
geom_bar(stat="identity") +
geom_text(aes(label=diff), vjust=-0.5) +
labs(x = "", y="change") +

ggtitle("German elections results (#MP change)") +

## AfD is brown regardless official party colors :-)
scale_fill_manual(values=c("#8B4513", "#56B4E9",
"yellow", "green", "red", "#ff6666") )
url | Mon, 25/09/2017 08:43 | tagi: , ,
Przed Żuławy w Koło 2017

Jutro planuję przejechać 140km biorąc udział w imprezie pn. Żuławy wKoło 2017. Niby Żuławy a profil trasy sugeruje jakieś istotne wzniesienie w okolicach 25--40km:

Uważnie przyjrzenie się liczbom (zwłaszcza na osi OY) pozwala stwierdzić, że jest to złudzenie, wynikające z różnicy w jednostkach miary obu osi (kilometry vs metry). W rzeczywistości góra tam jest symboliczna o czym można się przekonać robiąc wykres nachyleń. Żeby pozbyć się przypadkowych błędów związanych z niedokładnością pomiaru oryginalne 673 punktowe dane zostały zmienione na 111 punktowe (uśrednienie minimum 1 km) lub 62 punktowe (uśrednienie minimum 2 km).

Przy czym uśrednienie minimum $x$ oznacza obliczenie nachylenia dla najkrótszego odcinka kolejnych $n$ punktów z oryginalnego śladu GPX, który będzie dłuższy niż $x$.

Skrypty R/dane są tutaj. Oryginale ślady GPX/TCX skopiowane ze strony ŻwK są tutaj.

url | Sat, 23/09/2017 17:02 | tagi: , , ,
Czytelnictwo prasy

Punktem wyjścia są dane ze strony ZKDP (w formacie Excel.) Ponieważ pobieram je od pewnego czasu mam tego więcej niż jest na ww. stronie bo od stycznia 2015 do lipca 2017, czyli 31 plików. Ręczna konwersja byłaby zatem ciut za bardzo czasochłonna.

  for i in *.xls do
    oocalc --headless --convert-to csv $i ;
    # albo ssconvert -v $i `basename $i .xls`.csv ;
    done
  # Wyciągam dane dotyczące sprzedaży ogółem dla SE
  grep 'Super Ex' *.csv | awk -F ',' '{print $7} ' > se_sales.csv
  # Analogicznie dla innych tytułów

Uwaga: program ssconvert znajduje się w pakiecie gnumeric, oocalc to oczywiście składni Libre/OpenOffice.

Wielkości sprzedaży dla trzech najpoczytniejszych tytułów pokazują wykresy liniowe (pierwszy w tys egz. a drugi w procentach nakładu ze stycznia 2015 r.)


Sprzedaż w tys egz.

Sprzedaż w % poziomu ze stycznia 2015
library(ggplot2)
library(reshape2)

df <- read.csv("newspaper_sales_2015-17.csv", sep = ';',
               header=T, na.string="NA");

meltdf <- melt(df,id="month")

ggplot(meltdf,aes(x=month, y=value, colour=variable, group=variable)) +
  geom_line() +
  ylab(label="sales [ths]") +
  theme(legend.title=element_blank()) +
  scale_x_discrete (breaks=c("2015-01-01", "2015-06-01",
     "2016-01-01", "2016-06-01", "2017-01-01",  "2017-06-01"),
  labels=c("2015-01", "2015-06", "2016-01", "2016-06",
     "2017-01", "2017-06")  )

# https://stackoverflow.com/questions/10085806/extracting-specific-columns-from-a-data-frame
obs <- df[,c("month")]

normalize <- function(x) { return (x /x[1] * 100 )  }
dfN <- as.data.frame(lapply(df[-1], normalize))

# https://stackoverflow.com/questions/10150579/adding-a-column-to-a-data-frame
dfN["month"] <- obs

str(dfN)

meltdf <- melt(dfN,id="month")

# https://www.r-bloggers.com/what-is-a-linear-trend-by-the-way/
pN <- ggplot(meltdf,
 aes(x=month, y=value, colour=variable, group=variable)) + geom_line() +
 ylab(label="sales [ths]") +
 theme(legend.title=element_blank()) +
 stat_smooth(method = "lm", se=F) +
  scale_x_discrete (breaks=c("2015-01-01", "2015-06-01",
     "2016-01-01", "2016-06-01", "2017-01-01",  "2017-06-01"),
  labels=c("2015-01", "2015-06", "2016-01",
  "2016-06", "2017-01", "2017-06")  )

  pN

Spadek widoczny na wykresach można określić liczbowo na przykład szacując linię trendu:

# Trend liniowy
# http://t-redactyl.io/blog/2016/05/creating-plots-in-r-using-ggplot2-part-11-linear-regression-plots.html

# http://r-statistics.co/Time-Series-Analysis-With-R.html
seq = c (1:nrow(dfN))
dfN["trend"] <- seq

trendL.gw <- lm(data=dfN, gw ~ trend )
trendL.fakt <- lm(data=dfN, fakt ~ trend )
trendL.se <- lm(data=dfN, se ~ trend )

trendL.gw
trendL.fakt
trendL.se

Współczynniki trendu dla GW, Faktu i SE są odpowiednio równe -1.114 -0.6415 -0.4301, co należy interpretować następująco: przeciętnie z miesiąca na miesiąc nakład spada o 1,11%, 0,64% oraz 0,43% nakładu ze stycznia 2015 r., odpowiednio dla GW, Faktu i SuperExpresu.

Dane i wyniki są tutaj

url | Mon, 11/09/2017 17:33 | tagi: , , ,
Żuławy w koło 2016

Żuławy w koło to maraton rowerowy (czyli przejazd rowerem na dłuższym dystansie -- nie mylić z wyścigiem) organizowany od paru lat na Żuławach jak nazwa wskazuje. Sprawdziłem jak ta impreza wyglądała pod kątem prędkości w roku 2016. W tym celu ze strony Wyniki żUŁAWY wKOŁO 2016 ściągnąłem stosowny plik PDF z danymi, który następnie skonwertowałem do pliku w formacie XLS (Excel) wykorzystując konwerter on-line tajemniczej firmy convertio.pl. Tajemniczej w tym sensie, że nie znalazłem informacji kto i po co tą usługę świadczy.

Konwersja (do formatu CSV) -- jak to zwykle konwersja -- nie poszła na 100% poprawnie i wymagała jeszcze circa 30 minutowej ręcznej obróbki. Być może zresztą są lepsze konwertery, ale problem był z gatunku banalnych i wolałem stracić 30 minut na poprawianiu wyników konwersji niż 2 godziny na ustalaniu, który z konwerterów on-line konwertuje ten konkretny plik PDF (w miarę) bezbłędnie.

Po konwersji wypadało by sprawdzić (chociaż zgrubnie) czy wszystko jest OK.

## Czy każdy wiersz zawieraja 9 pól (powinien)
$ awk -F ';' 'NF != 9 {print NR, NF}' wyniki_zulawy_2016S.csv

## Ilu było uczestników na dystansie 140km?
$ awk -F ';' '$7 ==140 {print $0}' wyniki_zulawy_2016S.csv | wc -l
133

## Ilu było wszystkich (winno być 567 + 1 nagłówek)
$ cat wyniki_zulawy_2016S.csv | wc -l
568 # ok!

Przykładowy wykres pudełkowy

Do analizy statystycznej wykorzystano wykres pudełkowy (porównanie wyników na różnych dystansach) oraz histogram (rozkład średnich prędkości na dystansie 140km). BTW gdyby ktoś nie wiedział co to jest wykres pudełkowy to wyjaśnienie jest na rysunku obok. Objaśnienie: Me, $Q_1$, $Q_3$ to odpowiednio mediana i kwartyle. Dolna/górna krawędź prostokąta wyznacza zatem rozstęp kwartylny (IQR). Wąsy ($W_L$/$W_U$) są wyznaczane jako 150% wartości rozstępu kwartylnego. Wartości leżące poza ,,wąsami'' (nietypowe) są oznaczane kółkami.

Ww. wykresy wygenerowano następującym skryptem:

#
co <- "Żuławy wKoło 2016"
#
z <- read.csv("wyniki_zulawy_2016_C.csv", sep = ';',
  header=T, na.string="NA", dec=",");

aggregate (z$meanv, list(Numer = z$dist), fivenum)

boxplot (meanv ~ dist, z, xlab = "Dystans [km]",
    ylab = "Śr.prędkość [kmh]", col = "yellow", main=co )

## tylko dystans 140
z140 <- subset (z, ( dist == 140 ));

## statystyki zbiorcze
s140 <- summary(z140$meanv)
names(s140)

summary_label <- paste (sep='', "Średnia = ", s140[["Mean"]], 
  "\nMediana = ", s140[["Median"]],
  "\nQ1 = ", s140[["1st Qu."]],  "\nQ3 = ", s140[["3rd Qu."]],
  "\n\nMax = ", s140[["Max."]] )
# drukuje wartości kolumny meanv
# z140$meanv
# drukuje wartości statystyk zbiorczych
s140

# wykres słupkowy
h <- hist(z140$meanv, breaks=c(14,18,22,26,30,34,38), freq=TRUE, 
   col="orange", 
   main=paste (co, "[140km]"), # tytuł
   xlab="Prędkość [kmh]",ylab="L.kolarzy", labels=T, xaxt='n' )
# xaxt usuwa domyślną oś 
# axis definiuje lepiej oś OX
axis(side=1, at=c(14,18,22,26,30,34,38))
text(38, 37, summary_label, cex = .8, adj=c(1,1) )

Dane i wyniki są tutaj

url | Mon, 11/09/2017 08:10 | tagi: , , ,
Rozkład komisji obwodowych według liczby oddanych głosów

Rozkład komisji obwodowych według liczby oddanych głosów (na podstawie szczegółowych wyników wyborów do Sejmu RP, pobranych ze strony PKW -- por. Web scrapping protokołów wyborczych ze strony PKW):

komisje <- read.csv("komisje_glosy_razem.csv", sep = ';',  header=T, na.string="NA");
str(komisje);

hist(komisje$glosy, breaks=seq(0, 3200, by=25), col="orange",
     freq=TRUE,main="Komisje wg liczby oddanych głosów",
     xlab="# głosów",ylab="# komisji (N = 27859)" )

mtext(text="https://github.com/hrpunio/Data/tree/master/sejm", 4, cex=0.7)
text(3200,100, "Me = 495\nQ1 = 265\nQ3 = 782...", 2, cex=0.7,  adj=c(0,0));

fivenum(komisje$glosy);

quantile(komisje$glosy, c(.10));
quantile(komisje$glosy, c(.05));
quantile(komisje$glosy, c(.90));

url | Sat, 21/11/2015 19:23 | tagi: , , , ,
Wiek posłów ósmej kadencji Sejmu RP


Na stronie www.sejm.gov.pl już dziś pojawiły się strony o nowowybranych posłach 8 kadencji. Strony można ściągnąć na przykład takim oto prostym skryptem basha:

#!/bin/bash
# Przykładowy URL: http://www.sejm.gov.pl/Sejm8.nsf/posel.xsp?id=002&type=A
padtowidth=3
for ((i=1;i<=460;i++)) ; do
  ## parametr id w URLu ma wartość 001--460
  ## za pomocą printf/tricku z padtowidth dodajemy wiodące zera:
  POSEL=`printf "%0*d\n" $padtowidth $i`
  wget 'http://www.sejm.gov.pl/Sejm8.nsf/posel.xsp?id='$POSEL'&type=A'\
     -O $POSEL.html
done

Na stronach na razie jest niewiele informacji, ale jest data urodzenia, liczba zdobytych głosów oraz okręg wyborczy z którego poseł został wybrany. Za pomocą prostych skryptów Perla można wydłubać te dane, dodać informacje o wieku/płci i zapisać w pliku CSV:

imnz;rokur;wiek;klub;miejsce;okreg;glosy;plec
Adam Abramowicz;1961-03-10;54;PiS;NA;7 Chełm;10500;M
Andrzej Adamczyk;1959-01-04;56;PiS;NA;13 Kraków;18514;M
...

Jak wygląda struktura wiekowa w poszczególnych klubach? (na poniższym wydruku symbole x.1, x.2, x.3, x.4 oraz x.5, to odpowiednio: wartość minimalna, pierwszy kwartyl, mediana, trzeci kwartyl oraz wartość maksymalna)

p <- read.csv("Sejm_8_u.csv", sep = ';',  header=T, na.string="NA");
boxplot (wiek ~ klub, p, xlab="Klub", ylab="Wiek", col='yellow')

aggregate (p$wiek, list(Klub = p$klub), fivenum)
aggregate (p$wiek, list(Klub = p$klub), na.rm=TRUE, mean)

A jak wyglądała średnia wieku w poszczególnych kadencjach Sejmu?

p <- read.csv("Sejm1-8.csv", sep = ';',  header=T, na.string="NA");
boxplot (wiek ~ kadencja, p, xlab = "Kadencja", ylab = "Wiek", col='yellow')

aggregate (p$wiek, list(Kadencja = p$kadencja), fivenum)

 Kadencja  x.1  x.2  x.3  x.4  x.5
1     1991 22.0 37.0 43.0 49.0 70.0
2     1993 24.0 39.0 45.0 50.0 74.0
3     1997 23.0 40.5 46.0 51.0 72.0
4     2001 26.0 43.0 49.0 54.0 78.0
5     2005 23.0 41.0 47.0 53.0 67.0
6     2007 22.0 41.0 48.0 54.0 78.0
7     2011 22.0 42.0 50.0 56.0 73.0
8     2015 23.0 41.5 51.0 59.0 77.0

aggregate (p$wiek, list(Kadencja = p$kadencja), na.rm=TRUE, mean)

  Kadencja        x
1     1991 43.19438
2     1993 45.21535
3     1997 46.42500
4     2001 48.28221
5     2005 46.55230
6     2007 47.32948
7     2011 48.86739
8     2015 49.74783

Dane pobrane ze strony http://www.sejm.gov.pl/Sejm8.nsf/poslowie.xsp?type=A są dostępne tutaj.

url | Thu, 12/11/2015 23:36 | tagi: , , , ,
There is a cold summer this year in Sopot
Temperature in Sopot in July 2015
Temp in Sopot/July 2015

The following CSV (on-demand generated from raw data with simple Perl script) file contains outdoor temperature registred every hour starting from 2010 (with DS18B20 sensor):

dayhr;No;y2010;y2011;y2012;y2013;y2014;y2015;day30
d070100;001;14.6;17.5;14.9;10.1;12.9;12.2;0
d070101;002;13.4;16.7;14.1;10.1;12.8;12.5;3600
d070102;003;12.8;16.3;14.3;10.2;12.7;12.1;7200

dayhr is a label and day30 denotes number of seconds from the beginning od the period (for the first observation day30 equals 0, for the second 3600 etc.) The chart was produced with the following custom R script:

require(ggplot2)
library(scales)
number_ticks <- function(n) {function(limits) pretty(limits, n)}

d <- read.csv("july-by-day.csv", sep = ';',  header=T, na.string="NA");

datestart <- ISOdate(2015, 7, 1, tz = "");
d30 <- datestart + d$day30;
d[,"d30"]  <- d30;

ggplot(d, aes(x = d30)) +
  geom_line(aes(y = y2015, colour = 'y2015'), size=.3) +
  geom_line(aes(y = y2014, colour = 'y2014'), size=.3) +
  geom_smooth(aes(y = y2015, colour = 'y2015'), size=1) +
  geom_smooth(aes(y = y2014, colour = "y2014"), size=1) +
  ylab(label="Temp [C]") +
  xlab(label="Observation") +
  scale_y_continuous(breaks=number_ticks(15)) +
  scale_x_datetime(breaks = date_breaks("5 days")) +
  theme(legend.title=element_blank()) +
  ggtitle("Temperature in July in Sopot") +
  theme(legend.position=c(.6, .9)) +
  theme(legend.text=element_text(size=12));

ggsave(file="Temp-M7-2015.pdf", width=15, height=8)

url | Fri, 31/07/2015 14:34 | tagi: , ,
Aksjomat Balcerowicza: im większe wpływy związków zawodowych, tym mniej miejsc pracy

TU density vs GDP

TU density vs emp. rate

TU density vs unemp. rate

Kontunuując minianalizę rozpoczętą w poprzednim wpisie, a dotyczącą zależności pomiędzy zatrudnieniem a uzwiązkowieniem (w związku ze śmiałą tezą L. Balcerowicza, że taka zależność istnieje i jest ujemna):

require(ggplot2)

## https://stats.oecd.org/Index.aspx?DataSetCode=UN_DEN
## http://stats.oecd.org/Index.aspx?DatasetCode=STLABOUR
## employment rate Q42012
d <- read.csv("union_density_and_gdp.csv", sep = ';',  header=T, na.string="NA");

## tu.density = ratio of  wage and salary earners
## that are trade union members, divided by the total number of wage and salary earners:
## gdppc = GDP per capita
ggplot(d, aes(d$tu.density, d$gdppc)) + geom_point() +
  geom_text(aes(label=d$iso),size=2.0, vjust=-0.35)  +
  xlab("TU density (%)") + ylab("GDPpc (tys USD)") +
  scale_colour_discrete(name="") +
  geom_smooth(method="lm", se=T, size=2)

lm <- lm(data=d, gdppc ~ tu.density ); summary(lm);

## employment rate vs tu.density:
ggplot(d, aes(d$tu.density,d$emprate)) + geom_point() +
  geom_text(aes(label=d$iso),size=2.0, vjust=-0.35)  +
  xlab("TU density (%)") + ylab("Empolyment rate (%)") +
  scale_colour_discrete(name="") +
  geom_smooth(method="lm", se=T, size=2);

lm <- lm(data=d, emprate ~ tu.density ); summary(lm);

## youth unemployment rate vs tu.density:
## http://www.oecd-ilibrary.org/employment/youth-unemployment-rate_20752342-table2
ggplot(d, aes(d$tu.density,d$yur)) + geom_point() +
  geom_text(aes(label=d$iso),size=2.0, vjust=-0.35)  +
  xlab("TU density (%)") + ylab("Youth unempolyment rate (%)") +
  scale_colour_discrete(name="") +
  geom_smooth(method="lm", se=T, size=2);

lm <- lm(data=d, yur ~ tu.density ); summary(lm)

Prosta regresja daje następujące rezultaty: zależność #1 pomiędzy GDP per capita a Trade Union Density jest słabo dodatnia (to już wiemy); zależność #2 pomiędzy współczynnikiem zatrudnienia a Trade Union Density też jest słabo dodatnia; zależność #3 pomiędzy stopą bezrobocia w grupie wiekowej 15--24 lat a Trade Union Density jest wprawdzie ujemna, ale statystycznie nieistotna (współczynnik $R^2$ do tego równy 1,4%).

Jak to wygląda graficznie widać na wykresach obok.

Zbiór danych jest do pobrania tutaj.

BTW: do konwersji pliku PDF na JPG wykorzystano:

convert -density 150 Rplots.pdf Rplots_%02d.png

Uwaga na koniec: zapis method="lm" jest bardziej poprawny niż method=lm zastosowany w poprzednim wpisie.

url | Tue, 26/05/2015 18:20 | tagi: , , , ,