Weblog Tomasza Przechlewskiego [Zdjęcie T. Przechlewskiego]


scrum
random image [Photo gallery]
Zestawienie tagów
1-wire | 18b20 | 1wire | 2140 | 3rz | alsamixer | amazon | anniversary | antypis | apache | api | applebaum | arm | armenia | astronomy | asus | atom.xml | awk | aws | bachotek | bakłażan | balcerowicz | balta | bash | berlin | bibtex | bieszczady | biznes | blogger | blogging | blosxom | bono | borne-sulinowo | breugel | bt747 | budapeszt | canon | cedewu | chello | chiller | chillerpl | chown | chujowetaśmy | ciasto | cmentarz | contour | cron | css | csv | curl | cycling | d54250wykh | dbi | 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 | food | fop | foto | france | francja | fripp | fuczki | fuji | fuse | gammu | garmin | gawk | gazwyb | gdańsk | gdynia | gender | geo | geocoding | georgia | gft | git | github | gmail | gmaps | gnokii | gnus | google | googlecl | googleearth | googlemaps | gotowanie | gphoto | gphoto2 | gps | gpsbabel | gpsphoto | gpx | gpx-viewer | greasemonkey | gruzja | grzyby | haldaemon | handbrake | 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 | kociewie kołem | kod | kolibki | komorowski | konwersja | krutynia | kuchnia | kurski | latex | latex2rtf | latex3 | lcd | legend | lenny | lesund | lewactwo | liberation | linksys | linux | lisp | lisrel | litwa | lizbona | logika | ltr | lubowla | lwp | lwów | m2wś | mapquest | mapsource | marvell | math | mathjax | mazury | mbank | mediolan | mencoder | mh17 | michalak | michlmayr | microsoft | monitor | mp4box | mplayer | ms | msc | mssql | msw | mswindows | mtkbabel | museum | muzyka | mymaps | mysql | nanopi | natbib | navin | nekrolog | neo | neopi | netbook | niemcy | niemieckie zbrodnie | nikon | nmea | nowazelandia | nuc | nxml | oauth | oauth2 | obituary | okular | olympus | ooffice | ooxml | opera | osm | otf | otftotfm | other | overclocking | panoramio | pdf | pdfpages | pdftex | pdftk | perl | photo | photography | picasa | picasaweb | pim | pine | pis | pit | plotly | pls | plugin | po | podróże | politics | polityka | polsat | portugalia | 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 | sernik | sheevaplug | sienkiewicz | signature | sks | skype | skytraq | smoleńsk | sqlite | srtm | sshfs | ssl | staszek wawrykiewicz | statistics | stats | statystyka | stix | stretch | suwałki | svg | svn | swanetia | swornegacie | szwajcaria | słowacja | tbilisi | terrorism | tex | texgyre | texlive | thunderbird | tomato | totalnaopozycja | tourism | tramp | trang | truetype | ttf | turystyka | tusk | tv | tv5monde | twitter | typetools | ubuntu | uchodźcy | udev | ue | ukraina | umap | unix | upc | updmap | ups | utf8 | varia | video | vienna | virb edit | vostro | wammu | wdc | wdfs | webcam | webdav | wh2080 | wiedeń | wikicommons | wilno | win10 | 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 | zdrowie | łeba | świdnica | żywność
Archiwum
O stronie
wykorzystywany jest blosxom plus następujące wtyczki: tagging, flatarchives, rss10, lastbuilddatexhtmlmime. Niektóre musiałem dopasować nieco do swoich potrzeb. Więcej o blosxom jest tutaj
Subskrypcja
RSS 1.0
Hojarskiej ekstremalnie skośny rozkład głosów

Uczono mnie, że współczynnik skośności aczkolwiek może przyjąć wartości większe od 3, to w praktyce taka sytuacja się nie zdarza. Dlatego z rezerwą podszedłem do obliczeń, z których wynikało, że dla pewnego zbioru danych wynosi on 14:

library(moments)
# Wynik 1-ki z listy komitetu Z.Stonogi (D.Hojarska, wybory 2015)
s <- read.csv("hojarska.csv", sep = ';', header=T, na.string="NA");
skewness(s$glosy)
[1] 14.08602

nrow(s)
[1] 657

sum(s$glosy)
[1] 1671

Pierwsza hipoteza jest taka, że formuła liczenia współczynnika może być egzotyczna. Ustalmy zatem jak toto jest liczone:

?skewness
Description:  
This function computes skewness of given data

Wiele to się nie dowiedziałem. Po ściągnięciu źródła można ustalić, że to współczynnik klasyczny czyli iloraz trzeciego momentu centralnego przez odchylenie standardowe do trzeciej potęgi. Zatem jak najbardziej klasyczny a nie egzotyczny. Sprawdźmy co wyliczy Perl dla pewności:

#!/usr/bin/perl

print STDERR "** Użycie: $0 plik numer-kolumny (pierwszy wiersz jest pomijany)\n";
print STDERR "** Domyślnie: wyniki kandydata nr1 na liście komitetu Z.Stonoga (okr.25 / D.Hojarska)\n";

$file = "hojarska.csv"; ## default file
$column = 1; ## first column = 0

if ( $#ARGV >= 0 ) { $file=$ARGV[0]; } 
if ( $#ARGV >= 1 ) { $column=$ARGV[1]; } 

open (S, "$file") || die "Cannot open $file\n";

print "\nDane z pliku $file (kolumna: $column):\n"; 
$hdr = <S> ; # wczytaj i pomin naglowek ($n będzie prawidłowe)
while (<S>) { chomp();
  @tmp = split (/;/, $_); 
  $v = $tmp[$column];
   push(@L, $v) ; 
   $sum += $v; 
   $Counts{"$v"}++;
   $n++;
}

# Wyznaczenie średniej:
$mean = $sum /$n; 

## Wyznaczenie dominanty:
## przy okazji wydrukowanie danych pogrupowanych
print "+--------+---------+--------+--------+--------+\n";
print "| Głosy  | Obwody  |  cumGł | cumGł% |  cumN% |\n";
print "+--------+---------+--------+--------+--------+\n";
$maxc = -1;
for $c (sort {$a <=> $b } keys %Counts ) { 
  $sumCum += $Counts{"$c"} * $c; 
  $nCum += $Counts{"$c"}; 
  if ($maxc_pos < $Counts{$c} ) { 
    $maxc_pos = $Counts{$c} ; $maxc_val = $c
  }
  printf "| %6d |  %6d | %6d | %6.2f | %6.2f |\n", 
    $c,  $Counts{$c}, $sumCum, $sumCum/$sum *100, $nCum/$n *100;
}
print "+--------+---------+--------+--------+--------+\n\n";

$mode = $maxc_val; # dominanta

## Wyznaczenie mediany:
$half = int(($#L +1)/2 );
@L = sort ({ $a <=> $b } @L) ; #numerycznie
##print "$half of @L\n";
$median = $L[$half];

printf "Średnie: x̄ = %.3f (N=%d) Me =%.3f D = %.3f\n", $mean, $n, $median, $mode;

## Odchylenia od średniej:
for my $l (@L) {##
  $sd1 = ($l - $mean) ; 
  $sd2 = $sd1 * $sd1; # ^2 
  $sd3 = $sd2 * $sd1; # ^3

  $sum_sd1 += $sd1; $sum_sd2 += $sd2;
  $sum_sd3 += $sd3;
}

# odchylenie std./3-ci moment:
$sd = sqrt(($sum_sd2 / $n)); 
$u3 = $sum_sd3/$n;

printf "Sumy (x-x̄): %.2f(¹) %.2f(²) %.2f(³)\n", 
   $sum_sd1, $sum_sd2, $sum_sd3;
printf "Rozproszenie: σ = %.3f µ³ = %.3f\n", $sd, $u3;

printf "Skośność: (x̄-D)/σ = %.2f", ($mean -$mode)/$sd;
printf " 3(x̄-Me)/σ) = %.2f", 3 * ($mean - $median)/$sd;
printf " (µ³/σ³) = %.2f\n", $u3/($sd*$sd*$sd);
##//

Uruchomienie powyższego skryptu daje w wyniku to samo. Moja pewność, że wszystko jest OK jest teraz (prawie że) stuprocentowa.

Wracając do R:

library(ggplot2);
p <-  ggplot(data = s, aes(x = glosy)) + geom_histogram(binwidth = 4)
p
breaks <- ggplot_build(p)$data
breaks

[[1]]
     y count   x xmin xmax      density      ncount    ndensity PANEL group
1  482   482   0   -2    2 0.1834094368 1.000000000 2628.000000     1    -1
2  140   140   4    2    6 0.0532724505 0.290456432  763.319502     1    -1
...

Pierwszy przedział jest definiowany jako (-2;2], ale z uwagi na wartość danych de facto liczy głosy w komisjach, w których Hojarska dostała 0--2 głosów (drugi to (2;6], itd) i ni cholery nie da się tego zmienić na coś bardziej zbliżonego do prawdy. Bo tak jak jest to faktycznie pierwszy przedział jest dwa razy węższy niż każdy pozostały. W szczególności prawdziwa średnia, w tym przedziale wynosi 1,259 głosa, zaś liczona jako środek przedziału przez liczebność wynosi 1.0 (482 *1/482). Wg formuły ggplota środek przedziału to zero zatem średnia też wyjdzie 0, tj. wartość jest wyznaczona z dużym/nieokreślonym błędem.

Dzieje się tak ponieważ: the histogram is centered at 0, and the first bars xlimits are at 0.5*binwidth and -0.5*binwidth. Dopiero gdy dane są dodatnie, to ggplot zaczyna od zera. No ale nasz zbiór zawiera zera i klops. Zmiana tego jest nietrywialna.

Zamiast ggplota moża użyć po prostu polecenia hist:

h <- hist(s$glosy, breaks=seq(0,max(s$glosy), by=4) )
h

$breaks
 [1]   0   4   8  12  16  20  24  28  32  36  40  44  48  52  56  60  64  68  72
[20]  76  80  84  88  92  96 100 104 108 112 116 120 124 128 132 136 140 144 148
[39] 152

$counts
 [1] 592  42   6   4   2   2   2   1   2   1   2   0   0   0   0   0   0   0   0
[20]   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   1

$density
 [1] 0.2252663623 0.0159817352 0.0022831050 0.0015220700 0.0007610350
 [6] 0.0007610350 0.0007610350 0.0003805175 0.0007610350 0.0003805175
[11] 0.0007610350 0.0000000000 0.0000000000 0.0000000000 0.0000000000
[16] 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000
[21] 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000
[26] 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000
[31] 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000
[36] 0.0000000000 0.0000000000 0.0003805175

$mids
 [1]   2   6  10  14  18  22  26  30  34  38  42  46  50  54  58  62  66  70  74
[20]  78  82  86  90  94  98 102 106 110 114 118 122 126 130 134 138 142 146 150

No teraz liczy (prawie) jak trzeba pierwszy przedział [0;4], drugi (4;8] itd. Prawie bo na upartego pierwszy przedział jest szeroki na 5 wartości a kolejne na 4. Eh te detale i te upierdliwe zero:-)

Przy okazji się zadumałem, a cóż to jest to density. W pierwszym podejściu, że to są częstości, ale nie, bo ewidentnie nie sumują się do 1:

  ?hist
  density: values f^(x[i]), as estimated density values. If
          `all(diff(breaks) == 1)', they are the relative frequencies
          `counts/n' and in general satisfy sum[i; f^(x[i])
          (b[i+1]-b[i])] = 1, where b[i] = `breaks[i]

No dość kryptycznie, ale można się domyśleć, że nie $\sum_i d_i =1$, ale $\sum_i d_i * w_i$, gdzie $w_i$, to szerokość przedziału $i$ (a $d_i$ oznacza gęstość w przedziale $i$ oczywiście). Pole pod krzywą będzie zatem równe 1, jak tego należy się spodziewać po gęstości.

Wprawdzie już to padło mimochodem wyżej, ale analizowany zbiór danych to liczba oddanych głosów na numer jeden na liście komitetu wyborczego Zbigniewa Stonogi w okręgu wyborczym 25 w wyborach do Sejmu w 2015 roku. Tym numerem jeden była Danuta Hojarska, która faktycznie dostała 152 głosy w jednej komisji (we wsi, w której mieszka) oraz 0--4 głosów w 90% pozostałych komisjach (w tym 0 głosów w 191/657.0 = 29% komisjach).

Analiza przestrzenna

Ze starych zapasów odgrzebałem plik CSV ze współrzędnymi obwodowych komisji wyborczych i połączyłem go z wynikami uzyskanymi przez Hojarską z zamiarem wyświetlenia tego na mapie za pomocą Google Fusion Tables (GFT):

obwod;coordexact;glosy;icon
101884;54.2701205 18.6317438;3;measle_white
101885;54.260714 18.6282809;0;measle_white
101886;54.262187 18.632082;2;measle_white
101887;54.257501 18.63147;1;measle_white
101888;54.25786 18.6306574;7;measle_grey
...

Kolumna icon zawiera nazwę ikony, która jest definiowana wg schematu: 0--3 głosy biala (measle_white); 4--9 głosów szara (measle_grey) 9-19 głosów żółta (small_yellow) 20--99 głosów czerwona (small_red) 100 i więcej głosów purpurowa (small_purple). Zestawienie dostępnych w GFT ikon można znaleźć tutaj. Konfigurowanie GFT aby korzystała z kolumny z nazwami ikon sprowadza się do wyklikania stosownej pozycji (Configure map →Feature map →Change feature styles →Use icon specified in a column.)

Przy okazji odkurzyłem zasób skryptów/zasobów używanych do geokodowania i/lub przygotowywania danych pod potrzeby GFT, w szczególności: joincvs.pl (łączy dwa pliki csv w oparciu o wartości n-tej kolumny w pierwszym pliku oraz m-tej kolumny w drugim); geocodeCoder0.pl (uruchamia geokodera Google). Oba skrypty i jeszcze parę rzeczy można znaleźć: tutaj.

url | Sun, 25/02/2018 17:38 | tagi: ,
time plot tygodniowej liczby twitów

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

require(ggplot2)

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

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

##str(d)

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

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

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

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

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

ggsave(file=filePDF)  

url | Wed, 21/02/2018 13:56 | tagi: , ,
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: , , ,
Numerologia stosowana (lies, damned lies and statistics)

,,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ń.

url | Mon, 28/12/2015 20:03 | tagi: , ,
Afera madrycka: taka tam analiza wyjazdów posłów 7 kadencji

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:

Edytuj→Znajdź i zamień
usuwam odstępy, tj. zamieniam spację na pusty napis
Format→Komórki
wybieram numer z dwoma miejscami po przecinku.

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ą:

  1. współrzędne celu podróży (w tym celu zamieniam adres Państwo, miasto na współrzędne geograficzne korzystając z geokodera Google);

  2. 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);

  3. 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.

url | Tue, 09/12/2014 19:09 | 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: , , ,
Średnia ruchoma aka moving average

Poniżej nie do końca sprawdzona procedura wyznaczania średniej ruchomej zwykłej, zaimplementowana w Perlu:

#!/usr/bin/perl
my @qq= (1, 2, 6, 7, 10, 11, 11, 12, 13); # 9 elementów

print "Oryginalne:: @qq\n";

$ref_av = moving_avg(\@qq, 3);
print "Uśrednione:: @{$ref_av}\n";

$ref_av = moving_avg(\@qq, 7);
print "Uśrednione:: @{$ref_av}\n";

# ## ## ## ## ## ## ##
sub moving_avg {
  my $wlista = shift; # wskaźnik do listy, która będzie uśredniona
  my $lokr = shift; # liczba okresów

  my $i; my $total ;

  my @lista_av = @{$wlista};

  for ($i=0; $i < $lokr; $i++) {  $total += ${$wlista}[$i] ; }

  ## trzeba sprawdzić czy jest OK; nie jestem pewien
  for ($i = $lokr; $i <= $#{$wlista} + 1; $i++) {
    $lista_av[$i] = $total / $lokr ;
    $total += ${$wlista}[$i] - ${$wlista}[$i - $lokr];
  }

  return \@lista_av; ## zwróć wskaźnik do listy z uśredn. wartościami 
}

url | Thu, 10/02/2011 19:54 | 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: ,
Rysowanie diagramów ścieżkowych za pomocą MetaPosta

Diagram ścieżkowy (path diagram)-- popularna notacja służąca do graficznego definiowania równań strukturalnych -- to rodzaj grafu, w którym kółka i prostokąty oznaczają węzły a strzałki krawędzie. Do tego zarówno kółka/prostokąty jaki i strzałki mogą być oznaczone symbolami. No i tu pech, bo tradycyjnie używa się greckich liter, do tego z subskryptami.

Path Diagrams

Dia to fajny program, ale przy bardziej zaawansowanych diagramach wychodzą jego ograniczenia, jak przykładowo niemożność umieszczania formuł matematycznych. Wstawianie symboli greckich próbowałem obejść modyfikując LOCALE i klawiaturę (da się) ale subskrypty/superskrypty to już się okazały nie do przejścia. Kombinowałem zatem jakby tu łatwo i przyjemnie zrobić diagramy ścieżkowe w czymś innym. W szczególności książka Skrondala i Rabe-Hesketh (Generalized latent variable modeling: multilevel, longitudinal, and structural equation models, ISBN: 9781584880004) -- ewidentnie składania LaTeXem -- zawiera bardzo schludne diagramy ścieżkowe. Podobno też w R coś tam się samo składa... Nic wszakże konkretnego nie ustaliłem i stanęło jak zwykle na MetaPoście.

Tutaj umieszczone bardzo proste makra pozwalają na rysowania diagramów ścieżkowych. Przykładowo narysowanie jednoczynnikowego modelu pomiaru, zawierającego trzy miary refleksyjne będzie wyglądało jakoś tak:

z1  = (40mm,20mm); %% srodek koła
%% Rysowanie czynnika eta_1 z błędem delta_1
ErCircle(z1, "$\eta_1$", "$\delta_1$", "r");

%% miary (refleksyjne):
z15 = (10mm,0mm);  
%% miara x_1 z błędem \epsilon_1
ErRectangle(z15, "$x_1$", "$\epsilon_1$", "l"); 
CircleToRectangle(z15,z1, "$\lambda_1$");

z10 = (10mm,10mm); 
%% miara x_2 z błędem \epsilon_2
ErRectangle(z10, "$x_2$", "$\epsilon_2$", "l"); 
CircleToRectangle(z10,z1, "$\lambda_2$");

z11 = (10mm,20mm); 
%% miara x_3 z błędem \epsilon_3
ErRectangle(z11, "$x_3$", "$\epsilon_3$", "l"); 
CircleToRectangle(z11,z1, "$\lambda_3$");

Bardziej zaawansowany model jest obok na rysunku. Może komuś się też się przyda...

url | Sun, 19/07/2009 00:34 | tagi: , ,