![[Photo gallery]](/staff/tp/icons/ico_gallery.png){kind=icon}
Niewątpliwie piramida wieku to rysunek okropnie transfobiczny. Statistics Canada rozwiązał ten problem.
Kto uważa się za kobietę jest liczony jako W, kto za mężczyznę jest liczony do M. Prosta sprawa. Dla statystyka liczy się dokument lub oświadczenie; a czy pytany ma fiuta albo brodę go (statystyka) nie obchodzi. Problem jednak jest z tymi co nie mogą się zdecydować:
Given that the non-binary population is small, data aggregation is necessary to protect the confidentiality of responses provided. Most 2021 Census information is disseminated using a two-category gender variable. In these cases, people in the category ,,non-binary persons'' are distributed in the other two gender categories and are denoted by the + symbol.
Przyznam że mocno niejasne. Ludzie niebinarni są distributed in the other two gender categories? Liczy się ich dwa razy? Albo zlicza do jakiś other two categories? Ale jakich... Mocno to mętne i żadnych other two tutaj nie widać (cf https://www12.statcan.gc.ca/census-recensement/2021/dp-pd/dv-vd/pyramid/index-en.htm) Zakładając więc, że liczy się takiego jednak dwa razy (raz jako kobietę a raz jako mężczyznę) mamy razem 36,991,985. W innym zestawieniu jest liczba ludności łącznie (https://www150.statcan.gc.ca/t1/tbl1/en/tv.action?pid=9810000101) 36,991,981. Różnica wynosi 4 osoby. To te persony binarne?
Tyle się narobić dla czegoś co jest 4 rzędy mniejsze od typowego błędu statystyczneg? Tylko ideologiczne zaślepienie może to wytłumaczyć...
Reuters donosi, że While Russia has confirmed the world's fourth largest tally of coronavirus cases, it has a relatively low death toll from the associated disease, COVID-19[...] But data released by the Rosstat State Statistics Service on Sept. 4 show there were 57,800 excess deaths between May and July, the peak of the outbreak. The figure was calculated by comparing fatalities over those three months in 2020 with the average number of May-July deaths between 2015 and 2019. The excess total is more than three times greater than the official May-July COVID-19 death toll of 15,955. Co mnie zmotywowało do poszukania danych źródłowych.
Okazało się że dość prosto jest je pobrać ze strony
Urzędu Statystycznego Federacji Rosyjskiej (https://eng.gks.ru/
a konkretnie https://showdata.gks.ru/finder/)
Ponieważ okazało się to aż tak proste, to pobrałem
nie tylko zgony ale także urodzenia. W okresie 2015--2020 (miesięcznie).
Dane są w formacie XSLX. Kolumny to miesiące wiersze poszczególne regiony i ogółem Federacja. Eksportuję do CSV używając LibreOffice i wybieram wiersz z danymi dla całej federacji. W rezultacie mam trzy wiersze: rok-miesiąc, urodzenia i zgony. Obracam moim super skryptem do transpozycji plików CSV:
month;urodzenia;zgony 2015-01-01;149269.99;174722.99 2015-02-01;144591.99;156737 2015-03-01;160974.99;175809.99 ...
Trzy wykresy liniowe przedstawiające dynamikę zgonów, urodzin i przyrostu naturalnego (czyli różnicy między urodzinami a zgonami)
![[Zdjęcie T. Przechlewskiego]](/personal/rysunki/home-icon.jpg){kind=icon}
library("dplyr")
library("ggplot2")
library("scales")
spanV <- 0.5
## przyrost naturalny
d <- read.csv("UZRT.csv", sep = ';', header=T, na.string="NA");
d$diff <- d$urodzenia - d$zgony
pz <- ggplot(d, aes(x= as.Date(month), y=zgony)) +
geom_line(color="steelblue", size=.6, alpha=.5) +
geom_point(color="steelblue", size=.8) +
## Trend
geom_smooth(method="loess", se=F, span=spanV, color="red3") +
## Skala na osi X-ów
scale_x_date( labels = date_format("%y/%m"), breaks = "4 months") +
xlab(label="") +
ylab(label="deaths") +
geom_vline(xintercept = as.Date("2015-07-01"), alpha=.25, size=1) +
geom_vline(xintercept = as.Date("2016-07-01"), alpha=.25, size=1) +
geom_vline(xintercept = as.Date("2017-07-01"), alpha=.25, size=1)+
geom_vline(xintercept = as.Date("2018-07-01"), alpha=.25, size=1) +
geom_vline(xintercept = as.Date("2019-07-01"), alpha=.25, size=1) +
geom_vline(xintercept = as.Date("2020-07-01"), alpha=.25, size=1) +
labs(caption="https://showdata.gks.ru/finder/") +
ggtitle("Russian Federation: Deaths (ths)",
subtitle= sprintf("red line: loess with span=%.2f ; gray vertical: july", spanV))
pu <- ggplot(d, aes(x= as.Date(month), y=urodzenia)) +
geom_line(color="steelblue", size=.6, alpha=.5) +
geom_point(color="steelblue", size=.8) +
geom_smooth(method="loess", se=F, span=spanV, color="red3") +
scale_x_date( labels = date_format("%y/%m"), breaks = "4 months") +
xlab(label="") +
ylab(label="births") +
labs(caption="https://showdata.gks.ru/finder/") +
geom_vline(xintercept = as.Date("2015-08-01"), alpha=.25, size=1) +
geom_vline(xintercept = as.Date("2016-08-01"), alpha=.25, size=1) +
geom_vline(xintercept = as.Date("2017-08-01"), alpha=.25, size=1)+
geom_vline(xintercept = as.Date("2018-08-01"), alpha=.25, size=1) +
geom_vline(xintercept = as.Date("2019-08-01"), alpha=.25, size=1) +
geom_vline(xintercept = as.Date("2020-08-01"), alpha=.25, size=1) +
ggtitle("Russian Federation: Births (ths)",
subtitle= sprintf("red line: loess with span=%.2f ; gray vertical: august", spanV))
pp <- ggplot(d, aes(x= as.Date(month), y=diff)) +
geom_line(color="steelblue", size=.6, alpha=.5) +
geom_point(color="steelblue", size=.8) +
geom_smooth(method="loess", se=F, span=spanV, color="red3") +
scale_x_date( labels = date_format("%y/%m"), breaks = "4 months") +
geom_vline(xintercept = as.Date("2015-05-01"), alpha=.25, size=1) +
geom_vline(xintercept = as.Date("2016-05-01"), alpha=.25, size=1) +
geom_vline(xintercept = as.Date("2017-05-01"), alpha=.25, size=1)+
geom_vline(xintercept = as.Date("2018-05-01"), alpha=.25, size=1) +
geom_vline(xintercept = as.Date("2019-05-01"), alpha=.25, size=1) +
geom_vline(xintercept = as.Date("2020-05-01"), alpha=.25, size=1) +
xlab(label="") +
ylab(label="balance") +
labs(caption="https://showdata.gks.ru/finder/") +
ggtitle("Russian Federation: natural balance (ths)",
subtitle= sprintf("red line: loess with span=%.2f ; gray vertical: may", spanV))
Wykres liniowy dla miesięcy jednoimiennych (MJ):
library(ggplot2)
library(dplyr)
library(tidyr)
library(scales)
options(dplyr.print_max = 1e9)
size0 <- 1.2
size1 <- 0.6
size2 <- 0.6
tt <- read.csv("UZRT.csv", sep = ';', header=T, na.string="NA");
tt$rok <- substr(tt$month, 1, 4)
tt$yyyymmdd <- as.Date(sprintf ("%s-%s", substr(tt$month, 6, 7),
substr(tt$month, 9, 10) ), format="%m-%d")
## Obliczenie przyrostu naturalnego
tt$diff <- tt$urodzenia - tt$zgony
## Obliczenie
tt.yy2020 <- tt %>% filter(rok==2020) %>% as.data.frame
tt.yy2019 <- tt %>% filter(rok==2019) %>% as.data.frame
tt.yy2018 <- tt %>% filter(rok==2018) %>% as.data.frame
tt.yy2017 <- tt %>% filter(rok==2017) %>% as.data.frame
tt.yy2016 <- tt %>% filter(rok==2016) %>% as.data.frame
tt.yy2015 <- tt %>% filter(rok==2015) %>% as.data.frame
pf <- ggplot() +
ggtitle("Russian Federation: births 2015--2020") +
geom_line( data = tt.yy2020, mapping = aes(x=yyyymmdd, y = urodzenia, colour = '2020'), alpha=.5, size=size0) +
geom_point( data = tt.yy2020, mapping = aes(x=yyyymmdd, y = urodzenia, colour = '2020'), alpha=.5, size=size0) +
geom_line( data = tt.yy2019, mapping = aes(x=yyyymmdd, y = urodzenia, colour = '2019'), alpha=.5, size=size1) +
geom_point( data = tt.yy2019, mapping = aes(x=yyyymmdd, y = urodzenia, colour = '2019'), alpha=.5, size=size0) +
geom_line( data = tt.yy2018, mapping = aes(x=yyyymmdd, y = urodzenia, colour = '2018'), alpha=.5, size=size1) +
geom_point( data = tt.yy2018, mapping = aes(x=yyyymmdd, y = urodzenia, colour = '2018'), alpha=.5, size=size1) +
geom_line( data = tt.yy2017, mapping = aes(x=yyyymmdd, y = urodzenia, colour = '2017'), alpha=.5, size=size1) +
geom_point( data = tt.yy2017, mapping = aes(x=yyyymmdd, y = urodzenia, colour = '2017'), alpha=.5, size=size1) +
###
geom_line( data = tt.yy2016, mapping = aes(x=yyyymmdd, y = urodzenia, colour = '2016'), alpha=.5, size=size1) +
geom_point( data = tt.yy2016, mapping = aes(x=yyyymmdd, y = urodzenia, colour = '2016'), alpha=.5, size=size1) +
geom_line( data = tt.yy2015, mapping = aes(x=yyyymmdd, y = urodzenia, colour = '2015'), alpha=.5, size=size2) +
geom_point( data = tt.yy2015, mapping = aes(x=yyyymmdd, y = urodzenia, colour = '2015'), alpha=.5, size=size1) +
scale_x_date( labels = date_format("%y/%m"), breaks = "2 months") +
labs(caption="https://showdata.gks.ru/finder/") +
ylab(label="ths") +
xlab(label="") +
labs(colour = "") +
theme(legend.position="top") +
theme(legend.text=element_text(size=10));
qf <- ggplot() +
ggtitle("Russian Federation: deaths 2015-2020") +
geom_line( data = tt.yy2020, mapping = aes(x=yyyymmdd, y = zgony, colour = '2020'), alpha=.5, size=size0) +
geom_point( data = tt.yy2020, mapping = aes(x=yyyymmdd, y = zgony, colour = '2020'), alpha=.5, size=size0) +
geom_line( data = tt.yy2019, mapping = aes(x=yyyymmdd, y = zgony, colour = '2019'), alpha=.5, size=size1) +
geom_point( data = tt.yy2019, mapping = aes(x=yyyymmdd, y = zgony, colour = '2019'), alpha=.5, size=size0) +
geom_line( data = tt.yy2018, mapping = aes(x=yyyymmdd, y = zgony, colour = '2018'), alpha=.5, size=size1) +
geom_point( data = tt.yy2018, mapping = aes(x=yyyymmdd, y = zgony, colour = '2018'), alpha=.5, size=size1) +
geom_line( data = tt.yy2017, mapping = aes(x=yyyymmdd, y = zgony, colour = '2017'), alpha=.5, size=size1) +
geom_point( data = tt.yy2017, mapping = aes(x=yyyymmdd, y = zgony, colour = '2017'), alpha=.5, size=size1) +
###
geom_line( data = tt.yy2016, mapping = aes(x=yyyymmdd, y = zgony, colour = '2016'), alpha=.5, size=size1) +
geom_point( data = tt.yy2016, mapping = aes(x=yyyymmdd, y = zgony, colour = '2016'), alpha=.5, size=size1) +
geom_line( data = tt.yy2015, mapping = aes(x=yyyymmdd, y = zgony, colour = '2015'), alpha=.5, size=size2) +
geom_point( data = tt.yy2015, mapping = aes(x=yyyymmdd, y = zgony, colour = '2015'), alpha=.5, size=size2) +
scale_x_date( labels = date_format("%y/%m"), breaks = "2 months") +
labs(caption="https://showdata.gks.ru/finder/") +
ylab(label="ths") +
xlab(label="") +
labs(colour = "") +
theme(legend.position="top") +
theme(legend.text=element_text(size=10));
qf <- ggplot() +
ggtitle("Russian Federation: natural balance 2015-2020") +
geom_line( data = tt.yy2020, mapping = aes(x=yyyymmdd, y = diff, colour = '2020'), alpha=.5, size=size0) +
geom_point( data = tt.yy2020, mapping = aes(x=yyyymmdd, y = diff, colour = '2020'), alpha=.5, size=size0) +
geom_line( data = tt.yy2019, mapping = aes(x=yyyymmdd, y = diff, colour = '2019'), alpha=.5, size=size1) +
geom_point( data = tt.yy2019, mapping = aes(x=yyyymmdd, y = diff, colour = '2019'), alpha=.5, size=size0) +
geom_line( data = tt.yy2018, mapping = aes(x=yyyymmdd, y = diff, colour = '2018'), alpha=.5, size=size1) +
geom_point( data = tt.yy2018, mapping = aes(x=yyyymmdd, y = diff, colour = '2018'), alpha=.5, size=size1) +
geom_line( data = tt.yy2017, mapping = aes(x=yyyymmdd, y = diff, colour = '2017'), alpha=.5, size=size1) +
geom_point( data = tt.yy2017, mapping = aes(x=yyyymmdd, y = diff, colour = '2017'), alpha=.5, size=size1) +
###
geom_line( data = tt.yy2016, mapping = aes(x=yyyymmdd, y = diff, colour = '2016'), alpha=.5, size=size1) +
geom_point( data = tt.yy2016, mapping = aes(x=yyyymmdd, y = diff, colour = '2016'), alpha=.5, size=size1) +
geom_line( data = tt.yy2015, mapping = aes(x=yyyymmdd, y = diff, colour = '2015'), alpha=.5, size=size2) +
geom_point( data = tt.yy2015, mapping = aes(x=yyyymmdd, y = diff, colour = '2015'), alpha=.5, size=size2) +
scale_x_date( labels = date_format("%y/%m"), breaks = "2 months") +
labs(caption="https://showdata.gks.ru/finder/") +
ylab(label="ths") +
xlab(label="") +
labs(colour = "") +
theme(legend.position="top") +
theme(legend.text=element_text(size=10));
Wreszcie wykres kombinowany liniowo-słupkowy dla miesięcy jednoimiennych (MJ).
Linie to 1) średnia liczby zgonów w latach 2015--2020,
2) średnia liczba zgonów w latach 2015--2020 plus/minus odchylenie
standardowe liczby zgonów (dla MJ); 3) liczba zgonów w roku 2020. Wykres słupkowy to
różnica między liczbą zgonów w roku 2020
a maksymalną liczbą zgonów w latach 2015-2019 ORAZ liczba zgonów z powodu COVID19
(z bazy ECDC). Ten skrypt korzysta z danych z pliku covid_ru.csv, który
powstał przez zagregowanie danych dziennych dostępnych ze strony ECDC (dla Federacji Rosyjskiej)
oraz policzenie średnich/odchyleń
dla jednoimiennych miesięcy z danych w pliku UZRT.csv
library(ggplot2)
library(dplyr)
library(tidyr)
library(scales)
options(dplyr.print_max = 1e9)
#
tt <- read.csv("covid_ru.csv", sep = ';', header=T, na.string="NA");
tt$diff <- tt$urodzenia - tt$zgony
cols <- c("mean19"="navyblue","mean+sd19"="#3591d1", "mean-sd19"="#3591d1",
"deaths20"='brown4', "diff19"="red", "c19d"='cyan')
pf <- ggplot(data=tt) +
ggtitle("Russian Federation: deaths 2015--2020",
subtitle='mean19/sd19: mean/sd 2015--19; diff19: 2020 - max(2015-019); deaths20: deaths 2020; c19d: C19 deaths') +
geom_line(mapping = aes(x=yyyymmdd, y = mean, colour = 'mean19'), alpha=.5, size=.8) +
geom_line(mapping = aes(x=yyyymmdd, y = mean + sd, colour = 'mean+sd19'), alpha=.5, size=.8) +
geom_line(mapping = aes(x=yyyymmdd, y = mean - sd, colour = 'mean-sd19'), alpha=.5, size=.8) +
geom_text(mapping = aes(x=yyyymmdd, label=difflabel, y=diff), vjust=-1.0, size=3 ) +
##geom_line(mapping = aes(x=yyyymmdd, y = c19d, colour = 'c19d'), alpha=.5, size=1.2) +
labs(caption = "Source: https://showdata.gks.ru/finder/; https://www.ecdc.europa.eu/en/covid-19-pandemic") +
##
geom_bar(mapping = aes(x=yyyymmdd, y = diff, fill = 'diff19' ), color='brown', stat="identity", alpha=.25) +
geom_bar(mapping = aes(x=yyyymmdd, y = c19d, fill = 'c19d' ), color ='navyblue', stat="identity", alpha=.25) +
###
geom_line(mapping = aes(x=yyyymmdd, y = deaths, colour = 'deaths20'), alpha=.5, size=1.2) +
scale_x_date( labels = date_format("%y/%m"), breaks = "2 months") +
ylab(label="ths") +
xlab(label="") +
##labs(colour = "Legend") +
scale_colour_manual(name="Lines: ", values=cols) +
scale_fill_manual(name="Bars: ", values=cols) +
theme(legend.position="right") +
theme(legend.text=element_text(size=12));
Skrypt pomocniczy liczący średnie 2015--2019 dla miesięcy jednoimiennych itp rzeczy
library(dplyr)
library(tidyr)
## Agreguje zgodny wg miesięcy jednoimiennych (MJ) dla lat 2015-19
## drukuje średnią/sd/max dla MJ
## Drukuje z2020 - max(z2015--2019)
tt <- read.csv("UZRT.csv", sep = ';', header=T, na.string="NA");
tt$rok <- substr(tt$month, 1, 4)
tt$yyyymmdd <- as.Date(sprintf ("%s-%s", substr(tt$month, 6, 7), substr(tt$month, 9, 10) ), format="%m-%d")
tt$diff <- tt$urodzenia - tt$zgony
tt.yy2020 <- tt %>% filter(rok==2020) %>% as.data.frame
tt.yy2019 <- tt %>% filter(rok==2019) %>% as.data.frame
tt.yy2018 <- tt %>% filter(rok==2018) %>% as.data.frame
tt.yy2017 <- tt %>% filter(rok==2017) %>% as.data.frame
tt.yy2016 <- tt %>% filter(rok==2016) %>% as.data.frame
tt.yy2015 <- tt %>% filter(rok==2015) %>% as.data.frame
## max zgony
tt.yy2020$max.so.far <- pmax (tt.yy2019$zgony, tt.yy2018$zgony,
tt.yy2017$zgony, tt.yy2016$zgony, tt.yy2015$zgony)
## średnia
tt.yy2020$mean <- (tt.yy2019$zgony + tt.yy2018$zgony
+ tt.yy2017$zgony + tt.yy2016$zgony + tt.yy2015$zgony)/5
## odchylenie std:
tt.yy2020$sqmean <- (tt.yy2019$zgony^2 + tt.yy2018$zgony^2 +
tt.yy2017$zgony^2 + tt.yy2016$zgony^2 + tt.yy2015$zgony^2)/5
tt.yy2020$sq <- sqrt(tt.yy2020$sqmean - tt.yy2020$mean^2)
tt.yy2020$diff20 <- tt.yy2020$zgony - tt.yy2020$max.so.far
sprintf ("%s;%.2f;%.2f;%.2f;%.2f",
tt.yy2020$yyyymmdd, tt.yy2020$diff20, tt.yy2020$sqmean, tt.yy2020$sq, tt.yy2020$max.so.far )
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ć.
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.)
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
Ż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!
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
,,Po dziennik Fakt (Ringier Axel Springer Polska) w okresie od maja do października 2015 roku sięgało 10,33%. Drugie miejsce zajęła Gazeta Wyborcza (Agora), której poziom czytelnictwa wzrósł w ciągu roku o 0,39 pkt. proc. do 8,58%.'' (za Gazeta Wyborcza dogania...).
Na samym końcu tekstu: ,,Podane liczby to wskaźniki CCS -- Czytelnictwo Cyklu Sezonowego''.
Ki ch#j ten CCS? Bo liczby wydają się mooocno napompowane. Okazuje się, że CCS oznacza jaki procent w populacji stanowią osoby, które zetknęły się z tym tytułem choćby raz w okresie cyklu sezonowego... (za Wskaźniki PBC). Jeszcze żeby było wiadomo co oznacza ,,zetknęły się''? Zawinęły karpia w sklepie albo pocięły jako ersatz Toilettenpapier?
Skądinąd wiadomo (z innej strony na tym samym serwerze), że sprzedaż ww. jest licha (por. Spadła sprzedaż kioskowa wszystkich dziennikóww). Proste przemożenie 30136052 (Polacy w wieku 15--75 lat) x 10,33% i 8,58% odpowiednio daje w wyniku 3,1 i 2,6 mln czytelników przy 300 i 180 tys sprzedanych egz. co oznacza, że 1 egz. Faktu czyta 10 osób, a GW nawet 15... Kto chce niech wierzy.
Ja bym zamiast pytać ,,czy się zetknął'' pytał się o treść 1 dowolnego tekstu z gazety ,,z którą się miało zetknięcie''... Ciekawe ile by wtedy wyszło...
Dane z artykułów w formacie CSV są poniżej:
tytul;CCS2015;CCS2014;platne2015;druk2015;platne2014;druk2014;cz2015;cz2014;razem2015u;razem2014u;druk2015u;druk2014u Fakt;10.33;11.74;308616;305178;326748;323002;3113054;3537973;1008.7;1082.8;1020.1;1095.3 Gazeta Wyborcza;8.58;8.19;177447;106794;191127;113292;2585673;2468143;1457.2;1291.4;2421.2;2178.6 SuperExpress;3.58;4.55;145647;143954;153387;151735;1078871;1371190;740.7;893.9;749.5;903.7 Przegląd Sportowy;3.57;2.90;31484;30518;35546;34500;1075857;873946;3417.2;2458.6;3525.3;2533.2 Dziennik Gazeta Prawna;1.68;1.59;54089;9887;58027;11412;506286;479163;936.0;825.8;5120.7;4198.8 Rzeczpospolita;1.62;1.94;56560;7897;58056;9352;488204;584639;863.2;1007.0;6182.1;6251.5 PulsBiznesu;0.31;0.20;13314;1655;13996;1743;93422;60272;701.7;430.6;5644.8;3458.0 Parkiet Gazeta Giełdy;0.09;0.08;5332;1319;5015;1532;27122;24109;508.7;480.7;2056.3;1573.7
Gdzie: CCS2015 -- wskaźniki CCS dla 2015 r.; platne2015 -- rozpowszechnianie płatne razem w 2015;
druk2015 -- sprzedaż egz. wydań drukowanych; cz2015 = CCS2015 * 30136052;
razem2015u = platne2015 / druk2015 * 100%; druk2015u = druk2015 / druk2015 * 100%.
Dla roku 2014 analogicznie.
Ponadto:
Sprzedaż ogółem to suma sprzedaży egzemplarzowej wydań drukowanych, sprzedaży egzemplarzowej e-wydań oraz wszystkich form prenumeraty wydań drukowanych i prenumeraty e-wydań.
Rozpowszechnianie płatne razem to suma sprzedaży egzemplarzowej wydań drukowanych, sprzedaży egzemplarzowej e-wydań, prenumeraty wydań drukowanych i prenumeraty e-wydań oraz innych płatnych formy rozpowszechniania wydań drukowanych i innej płatnej dystrybucji e-wydań.
Scrapping various Web pages I managed to gather data on players participating in last 4 Rugby World Cups. Is there a trend in body mass of rugby players participating in RWC tournaments?
Using Plotly API via Perl script (described here Box-plot chart with plot.ly and Perl API) I can quickly plot series of boxplots:
# cpan install WebService::Plotly in case WebService::Plotly is not installed plotly_boxplot.pl -col=5 -by=0 -title='RWC players by weight' -sep=';' rwc1999-2015.csv
Resulting boxplots and data can be viewed here.
UWAGA: Ten tekst nie jest o polityce ale o [elementarnej] statystyce.
Media informowały, że posłowie PiS Adam Hofman, Mariusz A. Kamiński i Adam Rogacki wzięli na podróż do Madrytu na posiedzenie komisji Zgromadzenia Parlamentarnego Rady Europy po kilkanaście tysięcy złotych zaliczki, zgłaszając wyjazd samochodem; w rzeczywistości polecieli tanimi liniami lotniczymi. Ponieważ kontrola wydatków posłów jest iluzoryczna różnica pomiędzy kosztem podróży samochodem a samolotem [za dużo mniejsze pieniądze] miała stanowić dodatkowy przychód wyżej wymienionych. Według prokuratury, która wszczęła śledztwo, zachodzi podejrzenie popełnienia oszustwa.
Łapiąc wiatr w żagle [sprawa się upubliczniła tuż przed ostatnimi wyborami samorządowymi] koalicja rządząca w osobie Marszałka Sejmu RP Sikorskiego zarządziła audyt, którego efektem było udostępnienie m.in. dokumentu pn. Wyjazdy zagraniczne posłów VII kadencja (kopia jest tutaj).
Jak przystało na kraj, w którym od lat działa Ministerstwo cyfryzacji zestawienie jest w formacie PDF, zatem pierwszym ruchem była zamiana na coś przetwarzalnego. Wpisanie w google PDF+Excel+conversion skutkuje ogromną listą potencjalnych konwerterów. Bagatelizując skalę problemu spróbowałem dokonać konwersji narzędziami dostępnymi on-line, ale z marnym rezultatem (za duży dokument przykładowo; serwis za free jest tylko dla PDFów mniejszych niż 50 stron). W przypadku Convert PDF to EXCEL online & free coś tam skonwertował, nawet wyglądało toto na pierwszy rzut oka OK ale na drugi już nie: dokument niekompletny oraz nieprawidłowo zamienione niektóre liczby (przykładowo zamiast 837,50 zł w arkuszu jest 83750 -- 100 razy więcej!).
Ostatecznie skończyło się na ściągnięciu 30 dniowej wersji Adobe Acrobata Pro XI, który faktycznie sprawdził się w roli konwertera PDF→XLSX. Do konwersji wykorzystałem służbowego laptopa Elki wyposażonego w legalny Office 2010, na którym zainstalowałem ww. AA Pro XI. OOffice niby czyta XLSX, ale z koszmarnymi błędami, więc żeby dalej móc obrabiać arkusz w Linuksie wczytałem wynikowy XLSX do Excela 2010 po czym zapisałem go w (starszym) formacie XLS. Ten plik wyświetlił się w OO Calcu bez problemu.
Arkusz jest tak sformatowany, że 4 pierwsze komórki oraz są często wielowierszowe i scalone, zawierają bowiem liczbę porządkową, datę, miejsce i cel wyjazdu delegacji posłów. Po zamianie na plik CSV zawartość komórek scalonych pojawi się w pierwszym wierszu, a pozostałe będą puste. Prostym skryptem Perlowym mogę wypełnić puste komórki wg. algorytmu: jeżeli cztery pierwsze pola są puste, to skopiuj wartości ostatnich niepustych:
if ($tmp[0] eq '' && $tmp[1] eq '' && $tmp[2] eq '' && $tmp[3] eq '' ) { ... }
Pierwszy problem: wielowierszowe komórki z kolumn 1--4 nie zawsze są scalone. Czasem tekst jest podzielony na wiersze co psuje konwersję. Ręcznie scalam niescalone komórki (trochę to trwa). Przed scaleniem usuwam z kolumn 1--4 końce wiersza.
Drugi problem: część liczb nie jest liczbami z uwagi na użycie separatora tysięcy, który się zamienił w PDFie na odstęp (spację). Zatem zaznaczam kolumny zawierające różne pozycje kosztów po czym:
Po uporządkowaniu arkusza, zapisuję go w formacie CSV. Następnie prostym skryptem Perlowym zamieniam na taki plik CSV, w którym puste komórki są wypełniane zawartością z poprzednich wierszy. Kolumna Państwo - miasto jest kopiowana. Kopia jest zmieniana na jednoznaczne: Państwo, miasto (pierwszy-kraj, przecinek, pierwsze miasto z listy celów podróży -- żeby geokoderowi było łatwiej.)
Innym skryptem Perlowym dodaję do pliku CSV 3 kolumny, które zawierają:
współrzędne celu podróży (w tym celu zamieniam adres Państwo, miasto na współrzędne geograficzne korzystając z geokodera Google);
odległość w kilometrach pomiędzy punktem o współrzędnych
21.028075/52.225208 (W-wa, Wiejska 1) a celem podróży (obliczoną przy wykorzystaniu pakietu
GIS::Distance);
linię zdefiniowana w formacie KML o końcach 21.028075/52.225208--współrzędne-celu-podróży (do ewentualnego wykorzystania z Google Fusion Tables).
#!/usr/bin/perl
#
use Storable;
use Google::GeoCoder::Smart;
use GIS::Distance;
$geo = Google::GeoCoder::Smart->new();
my $gis = GIS::Distance->new();
my $GeoCodeCacheName = 'geocode.cache';
my $NewCoordinatesFetched=0; # global flag
my $SLEEP_TIME = 2 ;
my $coords_okr = "21.028075,52.225208"; # Warszawa = środek świata
my %GeoCodeCache = %{ retrieve("$GeoCodeCacheName") } if ( -f "$GeoCodeCacheName" ) ;
my ($wwa_lng, $wwa_lat) = split (",", $coords_okr);
my $linesNo = 0 ;
my $GCtotaluse = 1; # laczna liczba wywolan geocodera
while (<>) {
$linesNo++;
chomp(); $_ =~ s/[ \t]+;[ \t]+/;/g; ## usuń ew. niepotrzebne spacje
@line = split ";", $_; print STDERR "**$linesNo = $line[3] ... ";
# geokodowanie (uwaga na limit)
# Poprawki dla miejsc, których nie zna Google:
$line[3] =~ s/Erewań/Erywań/; ##
$line[3] =~ s/Sowayma/Madaba/; ## najbliższe miasto
$line[3] =~ s/Bołszowce/Iwano-Frankiwsk/; ## najbliższe miasto
my $coords = addr2coords( $line[3] );
($tmp_lat, $tmp_lng, $gcuse) = split " ", $coords;
if ($gcuse > 0) {$GCtotaluse++ ; }
$distance = $gis->distance($tmp_lat,$tmp_lng => $wwa_lat,$wwa_lng );
$distance_t = sprintf ("%.1f", $distance);
my $kml_line = "<LineString><coordinates>$tmp_lng,$tmp_lat $coords_okr</coordinates></LineString>";
print "$_;\"$coords\";$distance_t;\"$kml_line\"\n";
print STDERR "\n";
if ($GCtotaluse % 100 == 0 ) {# store every 100 geocoder calls
store(\%GeoCodeCache, "$GeoCodeCacheName");
print STDERR "\n... Cache stored. ***\n";
}
}
##
store(\%GeoCodeCache, "$GeoCodeCacheName");
## ## ## ####
sub addr2coords {
my $a = shift ;
my $r = shift || 'n';
my ($lat, $lng) ;
my $GCuse = 0;
##consult cache first
if (exists $GeoCodeCache{"$a"} ) {
print STDERR "Coordinates catched ... $a ";
($lat,$lng) = split (" ", $GeoCodeCache{"$a"} );
}
else {
print STDERR "Geocoding ... $a ";
my ($resultnum, $error, @results, $returncontent) = $geo->geocode("address" => "$a");
$GCuse = 1;
sleep $SLEEP_TIME; ## make short pause
$resultnum--;
$resultNo=$resultnum ;
if (resultNo > 0) { print STDERR "** Location $a occured more than once! **" }
if ($error eq 'OK') {
$NewCoordinatesFetched=1;
for $num(0 .. $resultnum) {
$lat = $results[$num]{geometry}{location}{lat};
$lng = $results[$num]{geometry}{location}{lng};
##print "*** LAT/LNG:$lat $lng ERROR: $error RES: $resultNo ***\n";
}
$GeoCodeCache{"$a"} = "$lat $lng"; ## store in cache
} else { print STDERR "** Location $a not found! due to $error **" }
}
if ($r eq 'r' ) { return "$lng,$lat,$GCuse"; } # w formacie KML
else { return "$lat $lng $GCuse"; }
}
Gotowy plik CSV zawierający zestawienie podróży jest dostępny tutaj.
Na podstawie zestawienia i z użyciem pakietu ggplot2 generują się takie oto śliczne wykresy.
Wszystkie podróże z zestawienie (N=1874; odpowiednio: koszt łączny, koszt transportu, długość w tys km):
Tylko podróże dla których koszt transportu był niezerowy (N=1423; odpowiednio: koszt łączny, koszt transportu, długość w tys km):
Poniższy skrypt R sumuje i drukuje wszystkie podróże każdego posła:
require(plyr)
d <- read.csv("W7RR_podroze_by_podroz1.csv", sep = ';', dec = ",", header=T, na.string="NA");
# Dodaj kolumnę której wartości to konkatenacja: "Poseł|Klub"
d[,"PosKlub"] <- do.call(paste, c(d[c("Posel", "Klub")], sep = "|"));
# Usuń wszystko za wyjątkiem tego co potrzeba:
d <- d[ c("PosKlub", "Klacznie", "Ktransp", "Dist") ];
# Sumowanie po PosKlub
PSums <- as.data.frame ( ddply(d, .(PosKlub), numcolwise(sum)) );
# Z powrotem rozdziel kolumnę "Poseł|Klub" na dwie
PSums <- as.data.frame ( within(PSums, PosKlub <-data.frame( do.call('rbind',
strsplit(as.character(PosKlub), '|', fixed=TRUE)))) )
# Drukuj
PSums;
Z pliku .Rout kopiuję zestawienie łącznych wydatków
posłów oraz łącznej pokonanej przez nich odległości:
PosKlub.X1 PosKlub.X2 KlacznieT KtranspT DistT
1 Adam Abramowicz PiS 4.02599 2.64595 1.3153
2 Adam Hofman PiS 119.55271 59.53315 26.1716
3 Adam Kępiński SLD 10.15754 7.93882 3.8069
4 Adam Kępiński TR 12.63098 8.02327 2.2107
...
Uwaga: kilkanaście nazwisk się powtarza ponieważ posłowie zmienili przynależność klubową w czasie trwania kadencji [Aby uwzględnić takich posłów sumowanie odbywało się po wartościach zmiennej zawierającej połączone napisy Poseł|Klub.]
Na podstawie takiego z kolei zestawienia i znowu z użyciem ggplot2 generują inne śliczne wykresy.
Uwaga: sumowane tylko podróże, dla których koszt transportu był niezerowy (N=1423; odpowiednio: koszt łączny, koszt transportu, długość w tys km):
Link do tabeli zawierającej zestawienie podróży w formacie Google Fusion Tables jest tutaj.
Dane + skrypty dostępne są także w: github.com/hrpunio/Data.
The following Perl script reads data from a CSV file and draws a series of Box-Plots. Usage:
perl plotly_boxplot.pl -col=number -by=number -title=TITLE
where: -col=number -- column number containig variable to plot,
-by=number -- column number containig grouping variable.
#!/usr/bin/perl
use WebService::Plotly;
use Getopt::Long;
# login to plotly script
require "$ENV{'HOME'}/bin/login2plotly.pl";
my $plotly = WebService::Plotly->new( un => $plotly_user, key => $plotly_key );
my $sep_sequence = ';';
my $col_number = -1;
my $by_col_number = -1;
my $chart_title='??Chart title??';
my $header='Y';
#my $boxpoints='outliers'; ## or all' | 'outliers' | False
my $USAGE="*** USAGE: -col=i -by=i -title=s -header=s -sep=s FILE *** \n";
# plot values from column 'col' grouped by column 'by'. If header is Y skip first row in data.
# Add title 'title'. Columns in csv data are separated by 'sep' (default ';')
GetOptions("col=i" => \$col_number, "by=i" => \$by_col_number, "title=s" => \$chart_title,
'header=s' => \$header, 'sep=s' => \$sep_sequence, );
##'boxpoints=s' => \$boxpoints ) ; ## this option not work!
if (($col_number == -1 ) || ($by_col_number == -1) ) { print $USAGE }
while (<>) { chomp ($_); $nr++;
if (($nr < 2) << ( $header eq 'Y' ) ) { next }
$_ =~ s/"//g;
my @fields = split(/$sep_sequence/, $_);
push @{$data{$fields[$by_col_number]}}, $fields[$col_number];
# http://stackoverflow.com/questions/3779213/how-do-i-push-a-value-onto-a-perl-hash-of-arrays
}
my @variants = sort keys %data;
print STDERR "*** No of rows scanned: $nr ***\n";
print STDERR "*** Groups found: @variants ($boxpoints) \n";
for $k (keys %data ) { print "$k"; push (@boxes, { y =>$data{$k}, type => 'box', #'boxpoints' => 'none',
name => "$k" } ) }
my $layout = { 'title' => $chart_title };
my $response = $plotly->plot(\@boxes, layout => $layout );
my $url = $response->{url};
my $filename = $response->{filename};
print STDERR "*** done: filename: '$filename' url: '$url' ***\n"
Example: Age of Nobel Prize winners by discipline (grouping wariable) plot.ly/~tomasz.przechlewski/28/
While GDP per capita is mainly an indicator reflecting the level of economic activity, Actual Individual Consumption (AIC) per capita is an alternative indicator better adapted to describe the material welfare situation of households.
GDP and AIC per capita in PPS, EU27 = 100
+-------------------------------------------------+
GDP per capita AIC per capita
2009 2010 2011 2009 2010 2011
+-------------------------------------------------+
EU27 100 100 100 100 100 100
Euro area (EA17) 109 108 108 107 107 107
+-------------------------------------------------+
Luxembourg 255 267 271 144 141 140
Netherlands 132 131 131 118 114 113
Ireland 130 129 129 103 103 101
Austria 125 127 129 116 118 119
Sweden 120 124 127 116 114 116
Denmark 123 128 125 116 116 113
Germany 115 119 121 115 117 120
Belgium 118 119 119 109 111 111
Finland 114 113 114 110 111 112
United Kingdom 111 111 109 121 120 118
France 109 108 108 113 113 113
Italy 104 101 100 103 102 101
Spain 103 99 98 96 95 94
Cyprus 100 97 94 102 99 98
Malta 83 85 85 85 83 84
Slovenia 87 84 84 81 80 81
Czech Republic 83 80 80 73 71 71
Greece 94 87 79 104 97 91
Portugal 80 80 77 83 84 81
Slovakia 73 73 73 72 71 70
Estonia 63 63 67 58 56 58
Lithuania* 55 57 66 63 61 70
Hungary 65 65 66 62 60 61
Poland 61 63 64 64 67 69
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Latvia 54 54 58 52 53 57
Romania 47 47 49 46 46 47
Bulgaria 44 44 46 43 43 45
+-------------------------------------------------+
Norway 177 181 186 134 136 135
Switzerland 150 154 157 128 129 130
Iceland 120 112 111 111 106 107
Croatia 62 59 61 58 57 59
Turkey 46 50 52 51 54 57
Montenegro 41 42 42 50 52 53
FYR Macedonia 36 36 35 41 40 40
Serbia 36 35 35 44 44 43
Albania** 28 27 30 32 30 34
BosniaHerzegovina 31 30 30 37 36 36
+-------------------------------------------------+
Source: Eurostat newsrelease 180/2012, December 2012 (http://epp.eurostat.ec.europa.eu)
* 2011 population figures adjusted on the basis of the 2011 Census. Therefore the per capita indices for 2011 are not entirely comparable with previous years due to this break in time series.
** Figures for all years based on Eurostat estimate of GDP
The euro area (EA17) consists of Belgium, Germany, Estonia, Ireland, Greece, Spain, France, Italy, Cyprus, Luxembourg, Malta, the Netherlands, Austria, Portugal, Slovenia, Slovakia and Finland.
Source: Eurostat
Click on image to enlarge
The collapse of production during April-2008--April-2009 period (left figure) corresponds to a significant rise in unemployment (right figure), however in April 2009, when production was at the lowest level, unemployment rate was approx 9% while currently is is approx 12%.
Google fusion tables another excercise.
Two data sets describe football players who plays in Polish t-Mobile ekstraklasa (1st division) and Pierwsza Liga (2nd division) in 2011/2012 (autumn).
To show from where the player came a straight line is drawn from a player's birthplace to club's stadium, the player plays for.
Figure 1. 1st division.
Figure 2. 2nd division.
Players from 2nd division seems to be born closer to the clubs they play for:-)
Warning: in considerable number of cases the geocoding as performed by Google maybe wrong due to poor data quality--have no time to check/correct.
Yet another excercise, which tests Google Fusion Tables.
The data set contains--among other things--birth place for 600+ rugby players who last year take part in Rugby World Cup in New Zealand. The raw data is available here and here.
On the following diagram (cf. Figure 1) the players are mapped by column containing birthplace coordinates
Figure 1. World's concentration of top Rugby Union players.
The map do not show which player plays for which country. To show that a straight line is drawn from each player's birthplace to the country's capital, the player plays for (cf. Figure 2).
Figure 2. The origin of top Rugby Union players by federation.
Some lines look strange and the problem is particularly evident around New-Zealand-Samoa-Tonga-Fiji. For example it seems that many Samoan players were born at high ocean (cf. Figure 3).
Figure 3. The origin of top Rugby Union Samoan players.
BTW mapping Samoans one-by-one is OK (try it). The problem is when all rows are mapped together by Visualize→Map function.
Accidentaly around 600 miles to the East of New Zealand lies antipodal meridian ie. a meridian which is diametrically opposite the Greenwich meridian (a pime one). The longitude of points lying on antipodal meridian can be referenced (at least in GoogleMaps) both as -180° or 180° (ie. -36,-180° and -36,+180° refers to the same place). Perhaps it is the cause of observed errors...
Wprawdzie Sejm RP szybkimi krokami zmierza w stronę wybranego w 1936 r. Reichstagu określanego jako ,,najlepiej opłacany męski chór ma świecie'' (The best paid male chorus in the world -- określenie użyte podobno w amerykańskim czasopiśmie The Literary Digest) i być może już w kadencji ósmej ten stan ideału zostanie osiągnięty. Ale zanim to nastąpi zbierzmy dane dotyczące posłów wybranych do Sejmu 7 kadencji.
Zaczynam od ściągnięcia listy stron bibliograficznych posłów:
## Najpierw zestawienie posłów:
wget -U XX http://www.sejm.gov.pl/sejm7.nsf/poslowie.xsp?type=A -O lista-poslow.html
## Potem dla każdego posła z zestawiena:
cat lista-poslow.html | perl -e '
undef $/;
$_ = <>;
s/\n/ /g;
while ($_ =~ m@href="/sejm7.nsf/posel.xsp\?id=([0-9]+)\&type=A@g ) {
$id = $1; ## id posła
$pos = "http://www.sejm.gov.pl//sejm7.nsf/posel.xsp?id=$id&type=A";
print STDERR $pos, "\n";
system ('wget', '-U', 'XXX', "-O", "7_$id.html", $pos);
sleep 3; ## --let's be a little polite--
}
Powyższe pobiera 460 stron, każda zawierająca biografię jakiegoś posła.
Z kolei poniższy skrypt wydłubuje co trzeba (data scraping) ze strony zawierających biografię pojedynczego posła (np. tego).
#!/usr/bin/perl
#
use Storable;
use Google::GeoCoder::Smart;
$geo = Google::GeoCoder::Smart->new();
# Domyślny kraj
my $GeoCodeCacheName = 'geocode.cache';
my $NewCoordinatesFetched=0; # global flag
my $kraj = 'Polska';
my $SLEEP_TIME = 3 ;
# Rok wyborów
my $baseyr = 2011;
my ($data_day, $data_mc, $data_yr);
undef $/;
$_ = <>;
s/\n/ /g;
# Retrieve geocode cash (if exists)
# http://curiousprogrammer.wordpress.com/2011/05/11/faking-image-based-programming-in-perl/
my %hash = %{ retrieve("$GeoCodeCacheName") } if ( -f "$GeoCodeCacheName" ) ;
if (m@<h2>([^<>]+)</h2>@) { $posel = recode($1); }
# zgadnij płeć patrząc na imię (jeżeli kończy się na a->kobieta):
my @tmp_posel = split " ", $posel;
if ( $tmp_posel[0] =~ m/a$/ ) { $sex='K' } else { $sex='M' }
# id posła:
if (m@http://www.sejm.gov.pl/sejm7.nsf/posel.xsp\?id=([0-9]+)@) {$id = $1; }
# liczba oddanych na posła głosów:
if (m@<p class="left">Liczba g\łos\ów:</p>[ \n\t]*<p class="right">([^<>]+)</p>@) {
$glosy = $1 ; }
# Data i miejsce urodzenia:</p><p class="right">29-06-1960<span>, </span>Gdańsk</p>
if (m@Data i miejsce urodzenia:[ \n\t]*</p>[ \n\t]*<p class="right">([0-9\-]+)<span>,[ \n\t]*\ [ \n\t]*</span>([^<>]+)</p>@ ) {
$data = recode ($1);
($data_day, $data_mc, $data_yr) = split "-", $data;
##print STDERR "R/M/D: $data_day, $data_mc, $data_yr\n";
$mce = recode ($2);
}
# Zawód:
if (m@<p class="left">Zaw\ód:</p>[ \nt]*<p class="right">([^<>]+)</p>@ ) { $zawod = recode($1) ; }
if (m@klub.xsp\?klub=([^"<>]+)@ ) { $klub = recode($1); }
# Klub poselski:
if (m@Okr\ęg wyborczy:</p>[ \t\n]*<p class="right">([^<>]+)</p>@) {
$okr = recode($1);
## Pierwszy wyraz to numer okręgu, reszta nazwa (może być wielowyrazowa)
($okrnr, $okrnz) = $okr =~ m/([^ \n\t]+)[ \n\t]+(.+)/; ## -- sprawdzić czy działa --
}
# Poprawienie błędów
if ($mce =~ /Ostrowiec Świetokrzyski/) {$mce = 'Ostrowiec Świętokrzyski' }
elsif ($mce =~ /Stargard Szaczeciński/) {$mce = 'Stargard Szczeciński' ; }
elsif ($mce =~ /Białograd/ ) {$mce = 'Białogard'; }
elsif ($mce=~ /Szwajcaria/) {$mce = 'Szwajcaria,Suwałki' }
elsif ($mce=~ /Kocierz Rydzwałdzki/) { $mce= "Kocierz Rychwałdzki" }
elsif ($mce=~ /Stąporów/) { $mce= "Stąporków" }
elsif ($mce=~ /Szmotuły/) { $mce= "Szamotuły" }
# Wyjątki:
if ($posel=~ /Arkady Fiedler/ ) {$kraj = 'Wielka Brytania' }
elsif ($posel=~ /Vincent-Rostowski/) {$kraj = 'Wielka Brytania' }
elsif ($mce=~ /Umuahia/) {$kraj = 'Nigeria' }
elsif ($posel=~ /Munyama/) {$kraj ='Zambia'; }
## Imię kończy się na `A' ale to facet:
elsif ($posel=~ /Kosma Złotowski/ ) {$sex = "M"}
# Sprawdź czy data jest zawsze w formacie dd-mm-yyyy:
if ($data_yr < 1900) { print STDERR "*** $data_yr: $data for $posel ($id)\n"; }
# geokodowanie (uwaga na limit)
my $coords = addr2coords($mce);
my $wiek = $baseyr - $data_yr;
my $coords_okr = addr2coords($okrnz, 'r'); ($tmp_lat, $tmp_lng) = split " ", $coords;
my $kml_line = "<LineString><coordinates>$tmp_lng,$tmp_lat $coords_okr</coordinates></LineString>";
print STDERR "$id : $sex : $posel : $wiek\n";
print "$id;$sex;$posel;$data;$wiek;$mce,$kraj;$coords;$klub;$glosy;$zawod;$okrnr;$okrnz;$kml_line\n";
## Retrieve geocode cash
if ($NewCoordinatesFetched) { store(\%hash, "$GeoCodeCacheName"); }
## ## ## ## ## ##
sub recode {
my $s = shift;
$s =~ s/\ / /g; $s =~ s/\ł/ł/g; $s =~ s/\ó/ó/g;
$s =~ s/\ń/ń/g; $s =~ s/\Ł/Ł/g; $s =~ s/\ź/ź/g;
$s =~ s/\ż/ż/g; $s =~ s/\Ś/Ś/g; $s =~ s/\Ż/Ż/g;
$s =~ s/\ę/ę/g; $s =~ s/\ą/ą/g; $s =~ s/\ś/ś/g;
$s =~ s/\ć/ć/g; $s =~ s/[ \t\n]+/ /g;
return $s;
}
## ## ## ## ## ##
sub addr2coords {
my $a = shift ;
my $r = shift || 'n';
my ($lat, $lng) ;
##consult cache first
if (exists $GeoCodeCache{"$a"} ) {
($lat,$lng) = split (" ", $GeoCodeCache{"$a"} );
}
else {
my ($resultnum, $error, @results, $returncontent) = $geo->geocode("address" => "$a");
$resultnum--;
$resultNo=$resultnum ;
if (resultNo > 0) { print STDERR "** Location $a occured more than once! **" }
if ($error eq 'OK') {
$NewCoordinatesFetched=1;
for $num(0 .. $resultnum) {
$lat = $results[$num]{geometry}{location}{lat};
$lng = $results[$num]{geometry}{location}{lng};
##print "*** LAT/LNG:$lat $lng ERROR: $error RES: $resultNo ***\n";
}
} else { print STDERR "** Location $a not found! due to $error **" }
}
$GeoCodeCache{"$a"} = "$lat $lng"; ## store in cache
sleep $SLEEP_TIME;
if ($r eq 'r' ) { return "$lng,$lat"; } # w formacie KML
else { return "$lat $lng"; }
}
Skrypt wydłubuje id posła (id), imię i nazwisko (imnz), datę urodzenia (data_ur),
miejsce urodzenia (mce_ur), skrót nazwy partii do której należy (partia),
liczbę oddanych na niego głosów (log), zawód (zawod),
numer okręgu wyborczego (nrokregu) i nazwę okręgu (nzokregu). Do tego
zgadywana jest płeć posłanki/posła oraz obliczany jest wiek w chwili wyboru jako $baseyr - $data_yr, gdzie
$data_yr to rok urodzenia.
Miejsce urodzenia oraz nazwa okręgu są geokodowane. Współrzędne miejsca urodzenia są zapisywane jako zmienna wsp. Jako zmienna odleglosc zapisywana jest linia definiowana (w formacie KML) jako: współrzędne miejsce urodzenia--współrzędne okręgu wyborczego.
<LineString><coordinates>lng1,lat1 lng2,lat2</coordinates></LineString>
Zamiana kupy śmiecia (tj. wszystkich 460 plików HTML) na plik CSV zawierający wyżej opisane dane sprowadza się teraz do wykonania:
for i in 7_*html ; do perl extract.pl $i >> lista_poslow_7_kadencji.csv ; done
Jako ciekawostkę można zauważyć, że strony sejmowe zawierają 6 oczywistych błędów: Ostrowiec Świetokrzyski, Stargard Szaczeciński, Białograd, Kocierz Rydzwałdzki, Stąporów oraz Szmotuły. Ponadto wieś Szwajcaria k. Suwałk wymagała doprecyzowania. Czterech posłów urodziło się za granicą a pan Kosma jest mężczyzną i prosta reguła: jeżeli pierwsze imię kończy się na ,,a'' to poseł jest kobietą zawiodła.
Ostatecznie wynik konwersji wygląda jakoś tak (cały plik jest tutaj):
id;plec;imnz;dat_ur;wiek;mce_ur;wsp;partia;log;zawod;nrokregu;nzokregu;odleglosc 001;M;Adam Abramowicz;10-03-1961;50;Biała Podlaska,Polska;52.0324265 23.1164689;PiS;\ 12708;przedsiębiorca;7;Chełm;<LineString><coordinates>23.1164689,52.0324265 23.4711986,51.1431232</coordinates></LineString> ...
Teraz importuję plik jako arkusz kalkulacyjny do GoogleDocs, a następnie, na podstawie tego arkusza tworzę dokument typu FusionTable. Po wizualizacji na mapie widać parę problemów -- ktoś się urodził niespodziewanie w USA a ktoś inny za Moskwą. Na szczęście takich omyłek jest tylko kilka...
Geocoder nie poradził sobie z miejscowościami: Model/prem. Pawlak, Jarosław/posłowie: Kulesza/Golba/Kasprzak, Orla/Eugeniusz Czykwin, Koło/Roman Kotliński, Lipnica/J. Borkowski, Tarnów/A. Grad. We wszystkich przypadkach odnajdywane były miejsca poza Polską -- przykładowo Orly/k Paryża zamiast Orli. Ponadto pani poseł Józefa Hrynkiewicz została prawidłowo odnaleziona w Daniuszewie ale się okazało, że ta miejscowość leży na Białorusi a nie w Polsce. Też ręcznie poprawiłem...
Ewidentnie moduł Google::GeoCoder::Smart ustawia
geocoder w trybie partial match, w którym to
trybie Google
intelligently handles incomplete input aka stara się być mądrzejszym od pytającego.
Może w trybie full match byłoby lepiej, ale to by wymagało doczytania.
Na dziś, po prostu poprawiłem ręcznie te kilka ewidentnie błędnych przypadków.
Po poprawkach wszystko -- przynajmniej na pierwszy rzut -- oka wygląda OK
Link do tabeli jest zaś tutaj.
Nawiasem mówiąc kodowanie
dokumentów na stronach www.sejm.gov.pl jest co najmniej dziwaczne.
Klasyczny zbiór danych pn. niemieckie łodzie podwodne z 2WWŚ wg. typu, liczby patroli, zatopionych statków i okrętów oraz -- dodane dziś -- miejsca zatopienia (źródło: www.uboat.net). W zbiorze są 1152 okręty (pomięto 14, które były eksploatowane przez Kriegsmarine ale nie były konstrukcjami niemieckimi), z tego 778 ma określoną pozycję, na której zostały zatopione (z czego kilkadziesiąt w ramach operacji DeadLight). Można zatem powiedzieć, że dane dotyczące miejsca zatopienia/zniszczenia są w miarę kompletne, bo np. Kemp (1997) podaje liczbę 784 U-bootów zniszczonych w czasie 2WWŚ (dane są różne w zależności od źródła.)
Jeżeli U-boot zaginął/przepadł bez wieści/został zniszczony, ale nie wiadomo dokładnie gdzie, to nie występuje na mapie. Operacja DeadLight to ta wielka kupa kropek ,,nad'' Irlandią....
Link do danych/mapy na serwerze Google jest tutaj.
Dane były tylko zgrubie weryfikowane, więc mogą zawierać błędy (ale ,,na oko'' jest ich niewiele)...
Kemp Paul, U-Boats Destroyed: German Submarine Losses in World Wars, US Naval Institute Press, 1997 (isbn: 1557508593)
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.
Program smartpls, który służy do modelowania ścieżkowego cząstkową metodą najmniejszych kwadratów (PLS path modelling), jest typowym przykładem aplikacji napisanej przez amatora:-). Program jest za free, ale autor nie udostępnia kodu źródłowego (pewnie się wstydzi).
W systemie FC8 działał a w FC11 przestał działać, zgłaszając na starcie błąd:
lookup error: ~/spls/configuration/org.eclipse.osgi/bundles/58/1/.cp/libswt-mozilla-gtk-3232.so: \ undefined symbol: _ZN4nsID5ParseEPKc
Pomogło doinstalowanie starego Firefoxa (w wersji 2). Przykładowo jeżeli
stary FF znajduje się w katalogu ~/bin/ff2, to program działa uruchomiony
w następujący sposób:
MOZILLA_FIVE_HOME=~/bin/ff2 spls
Napisałem o ww. problemie na forum spls -- może autor się przejmie i skompiluje spls w wersji korzystającej z FF3; albo --co też jest bardzo prawdopodobne -- zmiany umożliwiające kompilację z FF3 są tak duże, że mu się nie będzie chciało...
Arkusza kalkulacyjnego (a zwłaszcza Excela f-my Microsoft) nie należy używać do zadań bardziej skomplikowanych od podliczania budżetu domowego [McCulloughHeiser2008,Berger2007,McCullough2008] zobacz także [AlmironetAl2010,KeelingPavur2007,Nash2008,Su2008].
Literatura
Powyższy spis literatury w formacie BiBTeX.
Bagozzi i Dholakia [BagozziDholakia2006] stosują następującą procedurę ustalania trafności różnicowej (discriminant validity) skali wieloczynnikowej: 1) oszacowanie bazowego modelu CFA (swobodnie korelujące ze sobą czynniki); 2) oszacowanie ograniczonego modelu CFA, tj. modelu w którym korelacja pomiędzy dwoma czynnikami skali jest ustalona jako równa 1 (co oznacza, że czynniki te de facto stanowią jeden czynnik). Oszacowanie modelu ograniczonego dla każdej pary czynników skali; 3) ustalenie czy wartość różnicy statystyk χ2 jest istotna statystycznie.
Istotność statystyki χ2 świadczy iż jakość dopasowania modelu ograniczonego (krok 2) jest istotnie gorsza od modelu bazowego (krok 1), a zatem korelacja pomiędzy parą czynników jest mniejsza od 1, czyli czynniki te się różnią... Trafność dyskryminacyjna skali jest potwierdzona jeżeli wszystkie (a przynajmniej zdecydowana większość) wartości różnic są istotne statystycznie.
Jeżeli skala składa się z n czynników, to należy oszacować (n (n-1))/2 modeli ograniczonych (np. dla n=6, jest to 15 modeli), co jest pracochłonne. Skrypt LISRELa dla modelu CFA oraz modelu ograniczonego różni się zaś od modelu bazowego tylko jednym wierszem:
VA 1 PHI(i,j)
gdzie i oraz j są numerami odpowiednich czynników. Można to wszystko zautomatyzować, w sposób następujący: 1) oblicz model bazowy; z pliku OUT Lisrela pobierz wartość statystyki χ2; 2) w pętli dla każdej pary i,j wykonaj skrypt wyznaczający model ograniczony, pobierz wartość statystyki χ2; 3) oblicz różnicę i wydrukuj...
Skrypt Perlowy wykonujący powyższe (z przykładem wykorzystania) jest tutaj.
Literatura
Zaproponowane przez Fornella i Larckera [HairetAl98,FornellLarcker81] rzetelność łączna (composite reliability, CR) oraz przeciętna wariancja wyodrębniona (average variance extracted, AVE) stały się często wykorzystywanymi miarami rzetelności wewnętrznej oraz trafności zbieżnej. CR obliczana jest według następującej formuły [por. także zencaroline.blogspot.com/2007/06/composite-reliability]:
CRη = (∑i ληi )2 / ( (∑i ληi )2 + ∑i var(εi))
gdzie: var(εi) = 1 - (ληi2); ληi to wektor (zestandaryzowanych) ładunków czynnikowych dla zmiennej ukrytej η.
Minimalną akceptowaną wartością CR jest 0,7
Podobny do CR jest wskaźnik przeciętnej wariancji wyodrębnionej (average variance extracted, AVE), służący do oceny trafności zbieżnej (convergent validity). Jest on obliczany według następującej formuły:
AVEη = (∑i (ληi )2 ) / ( ∑i (ληi)2 + ∑i var(εi) )
Znaczenie symboli jest identyczne, jak we wzorze określającym rzetelność łączną. Minimalną akceptowaną wartością AVE jest 0,5.
Znając zestandaryzowane wartości ładunków czynnikowych ληi
policzenie CR oraz AVE jest proste.
Dodanie opcji SC (completely standardized solutions)
do polecenia OUTPUT spowoduje wydrukowanie zestandaryzowanych wartości ładunków czynnikowych. Poniższy
skrypt obliczy na tej podstawie CR/AVE dla każdej zmiennej ukrytej:
#!/usr/bin/perl
# Computes/prints Composite Reliability (CR) and Average Variance Extracted (AVE)
# from LISREL OUT file. Measures can load on several factors
#
# (c) 2010; t.przechlewski http://pinkaccordions.homelinux.org/staff/tp/
# GPL license
#
# The formulas for CR and AVE are as follows [cf. Hair, Anderson, Tatham and Black, Multivariate Data Analysis, 5th Ed., Pearson Education, p. 637]:
# CR = (\sum standardized loading)^2 / ( (\sum standardized loading)^2 + \sum indicator measurement error ),
# AVE = (\sum (standardized loading)^2 ) / ( \sum (standardized loading)^2 + \sum indicator measurement error )
# where: indicator measurement error = 1 - loading^2
# see also: http://zencaroline.blogspot.com/2007/06/composite-reliability.html
#
# usage: perl print_cr_and_ave lisrel-output-file
#
my $scan = 0; ## flag to figure out where we are
my $initial_latent_var_no = 0;
while (<>) {
chomp;
# We are looking for the line with `Completely Standardized Solution' (CSS), which
# starts the block containing relevant data
if (/Completely Standardized Solution/) { $scan = 1 ;
print STDERR "*** Found *** $_ ***\n"; next ; }
# After CSS line we look the line with LAMBDA-* (there are up to two such lines)
if ( $scan > 0 && /LAMBDA-[XY]/ ) { $scan++;
print STDERR "*** Found *** $_ (initial: $initial_latent_var_no) ***\n";
$_ = <> ; $_ = <> ; $_ = <> ; ## eat exactly next three lines
## ok we are about to scan LAMBDA-X/Y matrix
while (<>) { chomp;
if (/^[ \t]*$/) { # exmpty line ends LAMBDA-X/Y matrix
## before reading next block store the number of latent vars from the 1st block
## first latent var number in the second block = number of vars in the previous block +1
$initial_latent_var_no += $#loadings;
print STDERR "*** Initial var number = $initial_latent_var_no ***\n";
last ; ## OK, all rows in matrix was read...
}
s/- -/xxx/g; # change `- -' to `xxx'
@loadings = split ' ', $_;
# column number = latent var number ; column `0' contains measurement variable name
for ($l=1; $l <= $#loadings; $l++) {
if ($loadings[$l] !~ /xxx/) { ## store in hash
$Loadings{$l + $initial_latent_var_no }{$loadings[0]} = $loadings[$l];
}
}
}
}
}
print STDERR "*** Latent variables = $initial_latent_var_no ***\n";
print "=======================================================\n";
### Compute/print CR i AVE ### ### ### ###
for $l (sort keys %Loadings ) {
$loadings = $sqloadings = $errors = 0;
print STDERR "*** Xi/Eta: $l ***\n";
for $m ( sort keys %{ $Loadings{ $l }} ) {
$load = $Loadings{$l}{$m};
print STDERR "$m ($l) = $load\n";
$loadings += $load ;
$sqloadings += $load * $load ;
$errors += (1 - $load * $load);
}
$cr = ($loadings * $loadings) / ( ($loadings * $loadings) + $errors ) ;
$ave = $sqloadings / ($sqloadings + $errors ) ;
printf "Xi/Eta_%2d -> CR = %6.3f AVE = %6.3f\n", $l, $cr, $ave;
}
print "=======================================================\n";
Skrypt można także pobrać tutaj.
Literatura
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.