>> wybierz styl >> es :: ns :: bs

Weblog Tomasza Przechlewskiego [Zdjęcie T. Przechlewskiego] [[Ikona]]


scrum
random image [Photo gallery]
Zestawienie tagów
1-wire | 18b20 | 1wire | 2140 | 3rz | alsamixer | amazon | anniversary | antypis | apache | api | arm | armenia | astronomy | asus | atom.xml | awk | aws | bakłażan | balcerowicz | balta | bash | berlin | bibtex | bieszczady | biznes | blogger | blogging | blosxom | borne-sulinowo | breugel | bt747 | budapeszt | canon | cedewu | chello | chown | chujowetaśmy | cmentarz | contour | cron | css | csv | curl | d54250wykh | debian | dejavu | dhcp | dht22 | dia | docbook | dom | ds18b20 | dyndns | dynia | ebay | economy | ekonomia | elka | elm | emacs | emacs23 | english | ess | eu | excel | exif | exiftool | f11 | fc | fc11 | fc15 | fc5 | fc8 | fedora | fedora21 | fenix | ffmpeg | finepix | firefox | flickr | fontforge | fontspec | fonty | fop | foto | france | francja | fripp | fuczki | fuji | fuse | gammu | garmin | gawk | gazwyb | gdańsk | gdynia | gender | geo | georgia | gft | git | github | gmail | gnokii | gnus | google | googlecl | googleearth | googlemaps | gphoto | gphoto2 | gps | gpsbabel | gpsphoto | gpx | gpx-viewer | greasemonkey | gruzja | grzyby | haldaemon | handbrake | historia | history | hitler | holocaust | holokaust | hpmini | humour | iblue747 | ical | iiyama | ikea | imap | inkscape | inne | internet | j10i2 | javascript | jhead | k800i | kajak | kamera | kleinertest | kml | kmobiletools | knuth | kod | kolibki | komorowski | konwersja | krutynia | kuchnia | kurski | latex | latex2rtf | latex3 | lcd | legend | lenny | lesund | lewactwo | liberation | linux | lisp | lisrel | litwa | logika | ltr | lubowla | lwp | m2wś | mapsource | marvell | math | mathjax | mazury | mbank | mediolan | mencoder | mh17 | michalak | microsoft | monitor | mp4box | mplayer | ms | msc | msw | mtkbabel | museum | muzyka | mymaps | mysql | nanopi | natbib | navin | neo | neopi | netbook | niemcy | niemieckie zbrodnie | nikon | nmea | nowazelandia | nuc | nxml | oauth | oauth2 | obituary | okular | olympus | ooffice | ooxml | opera | otf | otftotfm | other | overclocking | panoramio | pdf | pdfpages | pdftex | pdftk | perl | photo | photography | picasa | picasaweb | pim | pine | pit | plotly | pls | plugin | po | politics | polityka | polsat | postęp | powerpoint | prelink | problem | propaganda | pstoedit | putin | python | r | radio | random | raspberry pi | refugees | relaxng | ridley | router | rower | rowery | rpi | rsync | rtf | ruby | rugby | russia | rwc | rwc2007 | rwc2011 | rzym | samba | sem | sheevaplug | sienkiewicz | signature | sks | skype | skytraq | smoleńsk | sqlite | srtm | ssl | statistics | stats | statystyka | stix | suwałki | svg | svn | swornegacie | szwajcaria | słowacja | tbilisi | terrorism | tex | texgyre | texlive | thunderbird | tomato | tourism | tramp | trang | truetype | ttf | turystyka | tusk | tv | tv5monde | twitter | typetools | ubuntu | uchodźcy | udev | umap | unix | upc | updmap | ups | utf8 | varia | video | vienna | virb edit | vostro | wammu | wdc | wdfs | webcam | webdav | wh2080 | wiedeń | wikicommons | wilno | windows | windows8 | wine | wioślarstwo | word | wordpress | wrt54gl | ws1080 | wtyczka | ww2 | www | wybory | wybory2015 | włochy | węgry | xemex | xetex | xft | xhtml | xine | xml | xmllint | xsd | xslt | xvidtune | youtube | yum | zakopane | zakupy | zdf | łeba | świdnica
Pobrania via google: [[Ikona]]
Archiwum
Inne blogi
N. Walsh | Morten H. Frederiksen | B. Clementson | prawo.vagla.pl | F. Hecker | M. Olson | J. Tennison | J. Clark | M. Nottingham | M. Shuttleworth | T. Isakowicz-Zalewski | J. Anglim | José A. Ortega Ruiz Modern Perl
Inne tematyczne
Ashwin Amanna | wiesia.nets.pl | Wojt | rwm.org.pl | DataBlog | Revolutions | Learning R | A. Gelman | C. Nel | J. Vogelgesang | ubl.xml.org/ | J.D. Long |
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
Ż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: , , , ,
Im większe wpływy związków zawodowych, tym mniej miejsc pracy
TU density vs GDP (OECD countries)
TU density vs GDP

Pan profesor Balcerowicz na finiszu kampanii prezydenckiej baaardzo mocno się zaangażował po stronie rządzącego układu, a to zaangażowanie przejawia się m.in. wzmożonym wypisaniem na Twitterze różnych mniej lub bardziej mądrych (zwykle mniej) sloganów (aka farmazonów). Np. "S" już poparła Dudę, który zabiega o poparcie OPZZ -- to zła wiadomość dla młodych. Im większe wpływy ZZ w państwie, tym mniej miejsc pracy..

Co szkodzi sprawdzić empirycznie tezę profesora? Pobrałem zatem ze strony stats.oecd.org dane dotyczące Trade Union Density (ratio of wage and salary earners that are trade union members, divided by the total number of wage and salary earners tj. udział fundusza płac związkowców do płac ogółem). A ze strony List of OECD countries by GDP per capita dane dotyczące GDP per capita (jakoś nie mogłem szybko odszukać tych liczb na stronie stats.oecd.org -- zakładam, że na wikipedia.org przepisano je bez błędów:-) Dane są z roku 2012.

require(ggplot2)

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

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)

Jak widać na wykresie Polska jest piąta od końca wśród krajów OECD pod względem GDP na głowę i szósta od końca jeżeli chodzi o wielkość uzwiązkowienia. Czepianie się związków w tej sytuacji (12,5% uzwiązkowienia w PL, podczas gdy przykładowo w Niemczech jest to 41.9%, a w Danii 67.2%) ma wszystkie znamiona obsesji podobnej przykładowo do popularnego wśród Palikotowców i innych antyklerykałów poglądu, iż jakoby Kościół Katolicki jest praprzyczyną wszelkiego zła (przynajmniej w PL).

Dodatkowo prosta regresja daje następujący rezultat: GDP = 0,25 tu.density + 30,5435, czyli 1% wzrost uzwiązkowienia daje 0,25 tys wzrostu GDP na głowę (dokładnie odwrotnie niż twierdzi Balcerowicz). Współczynnik przy zmiennej tu.density jest nawet istotny statystycznie (na poziomie 0,05) ale $R^2$ jest faktycznie bardzo marne -- 13%.

Zbiór danych jest do pobrania tutaj.

url | Sun, 17/05/2015 14:52 | tagi: , , , ,
Data sets for statistical education

IMHO there is acute lack of interesting data sets for statistical education purposes, particularly cross-sectional in nature and not related to biology/medicine/astronomy. The domains I am looking for are management, demography and economics but sport is OK too as many students are interested in sport.

I agree there is a lot of data published (mainly) by various government agencies but its almost 100% time series data. There are many test datasets too (like UC Irvine Machine Learning Repository), but they are usually boring, often categorical and short:-)

To remedy a little the lack of good data I compiled myself a few databases and I intend to add more in the future. At present my repository consists of 4 datasets concerning: -- 1100 or so German WW2 U-boots; -- 2200 or so Titanic passangers/crew members -- 1200 or so Rugby Union players who took part in RWC 2011 and RWC 2007 tournaments.

The repository is available at my github account.

url | Tue, 08/11/2011 19:52 | tagi: , , ,
Introduction to Data Technologies

Przypadkowo znalazłem całkiem grubą książkę nt. Introduction to Data Technologies. Na super rewelację toto nie wygląda, ale w wolnej chwili obejrzę dokładniej...

url | Thu, 25/11/2010 20:13 | tagi: ,
Metody ilościowe w R. Aplikacje ekonomiczne i finansowe

Metody ilościowe w R, Aplikacje ekonomiczne i finansowe, to książka napisana przez zespół autorów z Uniwersytetu Warszawskiego. Na okładce podpisani są Katarzyna Kopczewska, Tomasz Kopczewski oraz Piotr Wójcik. Na odwrocie strony tytułowej do tego dopisano także: Michała Ejdysa, Piotra Szymańskiego i Jakuba Świniarskiego. W paginie każdej strony parzystej pp. Kopczewska, Kopczewski i Piotr Wójcik są zaś określani zaś jako redaktorzy. Wydawca to CeDeWu.pl [pełna nazwa: CeDeWu wydawnictwa fachowe, która zakrawa na żart, o czym niżej]. Redaktor -- jeżeli w ogóle w tym CeDeWu, ktoś taki jest -- się nie podpisał.

Książka jest po prostu redakcyjnie fatalna a jej czytanie wręcz boli. Jak pisałem wyżej, można mieć wątpliwości czy ktokolwiek toto redagował. Wszystko spieprzono -- za przeproszeniem -- jedną decyzją redakcyjną, mianowicie ktoś wpadł na ,,genialny'' pomysł dziwacznego wyróżniania fragmentów programów w R. [A w tej książce circa połowa objętości to przykłady i fragmenty programów.] Otóż konstrukcje języka R w tekście są wyróżnione przez złożenie ich (wąskim) krojem bezszeryfowym w stopniu optycznie znacząco większym od kroju otaczającego tekstu. Tego się po prostu nie da czytać!

Żeby jeszcze było trudniej czytelnikowi, to stosując kursywę i pogrubienie, autorzy usiłują graficznie rozróżnić: zmienne, wartości zmiennych, nazwy poleceń, argumenty poleceń, nazwy pakietów. Pytanie po co? Co to daje oprócz mieniącej się w oczach pstrokacizny? Oczywiście taki skomplikowany schemat skutkuje multum błędów i niekonsekwencji (co raz jest pochyłe-grube, innym razem jest proste-chude). Kolejne kuriozum, to wyśrodkowane krótkich fragmentów składni umieszczonych w całości w oddzielnym akapicie (tzw. choinka, po co?). Jednocześnie w tych wyśrodkowanych akapitach nie ma już graficznego rozróżnienia co jest zmienną, nazwą polecenia itp... Pełen chaos... lepiej zaciemnić już nie można było... Nie prościej było wszystko złożyć Courierem, wyróżniając nieliteralne fragmenty np. kursywą? Cóś w stylu: merge(zbiórA, zbiórB, by.X=zmiennaA, by.Y=zmiennaB,all=TRUE).

Inne błędy to już pestka:

Wydawanie takich książek to wstyd, obciach i żenada. Przede wszystkim dotyczy to wydawnictwa (CeDeWu wydawnictwa fachowe), bo to ono odpowiada za redakcję merytoryczną i techniczną. Tak się kończy powielanie maszynopisu przygotowanego przez dyletantów i na dokładkę w programie, który do tego się średnio nadaje. Rekomendacja: Nie kupować, szkoda pieniędzy...

Katarzyna Kopczewska, Tomasz Kopczewski i Piotr Wójcik, Metody ilościowe w R, Aplikacje ekonomiczne i finansowe, CeDeWu Warszawa 2009, ISBN 978-83-7556-150-0.

url | Fri, 26/02/2010 08:11 | tagi: , , ,
Dyskretyzacja w R

Jeżeli plik data1.dat zawiera liczby rzeczywiste, to poniższy skrypt wykona następujacą transformację:

# zamień macierz liczb rzeczywistych z pliku data1.dat na 
# odpowiednie wartości na skali porządkowej...
categories <- 5;
df <- read.table("data1.dat", header=T);
dl <- lapply(df, function(x){ cut(x, breaks=categories, labels=FALSE)} );
write.table(data.frame(dl), "data1.txt", sep="\t", na="", row.names=F)

Dla każdej kolumny wyznaczy categories przedziałów o jednakowej szerokości, przyporządkuje każdą liczbę do określonego przedziału, przypisze tej liczbie numer tego przedziału. Polecenie write.table wypisze wyznaczone numery do pliku data1.txt.

url | Fri, 11/09/2009 23:04 | tagi: ,
Środowisko R i pakiet ESS: instalowanie i rozpoczęcie pracy

R to środowisko do obliczeń statystycznych i wchodzi w skład każdej praktycznie dystrybucji Linuksa. Zainstalować można go bez problemu używając yuma, jeżeli już wcześnie nie został zainstalowany domyślnie. Dokumentację w formacie html odnaleźć można w katalogu /usr/lib/R/html/.

Emacs ma wsparcie do R w postaci pakietu ESS. Instalowanie ESS jest proste: należy rozpakować i dodać do plików startowych Emacsa następujące dwa wiersze (katalog ~/.emacs-local/ess/lisp oczywiście należy dopasować do własnych ustawień):

(add-to-list 'load-path "~/.emacs-local/ess/lisp")
(require 'ess-site)

Uruchamianie ESS jest jakby nieco mniej oczywiste; być może nawet to co opisałem poniżej jest nieoptymalne. Startuję R z wnętrza Emacsa za pomocą M-x R Enter. Zostanie wyświetlone w minibuforze pytanie o katalog roboczy, np.:

ESS [S(R): R] starting data directory ...

Należy wybrać odpowiedni katalog. Po pewnej chwili Emacs przejdzie do bufora *R*, który umożliwia interaktywną pracę z R. W buforze *R* można działać w środowisku R z wnętrza Emacsa dzięki czemu pracuje się wygodniej: działa dopełnienie (Tab) oraz help (C-c C-v). Tyle, że w buforze *R* polecenia R i wyniki obliczeń są przemieszane i szybko można się pogubić. Lepiej pisać program (skrypt) R w osobnym buforze a wyniki oglądać w buforze *R* (ogólnie *R:numer-procesu*, jeżeli działamy z więcej niż jednym skryptem, tj. dla drugiego skryptu zostanie utworzony bufor *R:2*, dla trzeciego *R:3*, itd.). Aby to osiągnąć należy otworzyć (nowy) plik za pomocą standardowego polecenia C-x C-f. Plik powinien mieć rozszerzenie .r. Bufor przejdzie do trybu ESS co zostanie zasygnalizowane pojawieniem się napisu ESS w wierszu trybu (modeline).

ESS

W tym buforze także działa pomoc (C-c C-v) i dopełnianie (C-c C-Tab). Pojedynczy wiersz ze skryptu R można uruchamiać za pomocą ess-eval-line (C-c C-j; uwaga: polecenia podzielone na wiersze wymagają naciśnięcia C-c C-j dla każdego wiersza); cały blok poleceń zaś za pomocą ess-eval-region (C-c C-r). Drobna niedogodność to przechodzenie pomiędzy różnymi oknami: tematów pomocy, R oraz bufora ze skryptem R (ESS otwiera/zamyka okna mało ,,intuicyjnie''). Ponieważ skrypty R są krótkie dobrym pomysłem jest podział ekranu na pół (C-x 2) i wyświetlanie w drugim oknie bufora *R*.

ESS

Prosty przykład wykorzystania R do określenia związku między poziomem korupcji a sposobem głosowania w sprawie zaakceptowania przez ISO specyfikacji OOXML można znaleźć w Corrupt countries were more likely to support the OOXML document format (Kai Puolamäki). Rysunek obok pokazuje wykonanie skryptu R z ,,wnętrza'' Emacsa (jak widać nawet okno zawierające histogram też się ładnie wyświetliło).

url | Mon, 08/10/2007 11:59 | tagi: , , , , ,