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

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


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

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.

url | Tue, 29/09/2015 16:53 | 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: , , , , , , ,
Box-plot chart with plot.ly and Perl API

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/

url | Mon, 07/04/2014 14:01 | tagi: , , ,
GDP per capita in purchasing power standards of EU member states

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.

url | Mon, 17/12/2012 09:52 | tagi: , , , ,
Grim prospects of EU economy

Total production in EU 2012 Unemployment in EU 2012
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%.

url | Fri, 14/12/2012 09:34 | tagi: , , , ,
Where they came from?

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.

Raw data are available here: div#1 and div#2.

url | Thu, 09/02/2012 10:41 | tagi: , , , , ,
Google Fusion Tables map rendering issue?

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

url | Mon, 06/02/2012 10:54 | tagi: , , , , , ,
Trzecie ćwiczenie z Google Fusion Tables - sejm RP siódmej kadencji

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]+)\&amp;type=A@g ) {
    $id = $1; ## id posła
    $pos = "http://www.sejm.gov.pl//sejm7.nsf/posel.xsp?id=$id&amp;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\&#322;os\&oacute;w:</p>[ \n\t]*<p class="right">([^<>]+)</p>@) {
   $glosy = $1 ; }

# Data i miejsce urodzenia:</p><p class="right">29-06-1960<span>,&nbsp;</span>Gda&#324;sk</p>
if (m@Data i miejsce urodzenia:[ \n\t]*</p>[ \n\t]*<p class="right">([0-9\-]+)<span>,[ \n\t]*\&nbsp;[ \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\&oacute;d:</p>[ \nt]*<p class="right">([^<>]+)</p>@ ) { $zawod = recode($1) ; }

if (m@klub.xsp\?klub=([^"<>]+)@ ) { $klub = recode($1); }

# Klub poselski:
if (m@Okr\&#281;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/\&nbsp;/ /g; $s =~ s/\&#322;/ł/g; $s =~ s/\&oacute;/ó/g;
  $s =~ s/\&#324;/ń/g; $s =~ s/\&#321;/Ł/g; $s =~ s/\&#378;/ź/g;
  $s =~ s/\&#380;/ż/g; $s =~ s/\&#346;/Ś/g; $s =~ s/\&#379;/Ż/g;
  $s =~ s/\&#281;/ę/g; $s =~ s/\&#261;/ą/g; $s =~ s/\&#347;/ś/g;
  $s =~ s/\&#263;/ć/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.

url | Sun, 29/01/2012 21:31 | tagi: , , , , , ,
Drugie ćwiczenie z Google Fusion Tables - mapa zatopionych Ubootów

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

Literatura

Kemp Paul, U-Boats Destroyed: German Submarine Losses in World Wars, US Naval Institute Press, 1997 (isbn: 1557508593)

url | Wed, 11/01/2012 21:23 | 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: , , ,
Uruchomienie Smart PLS w Linuksie

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

url | Fri, 12/11/2010 09:39 | tagi: , , , ,
Dlaczego nie należy używać Excela

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

AlmironetAl2010
Marcelo G. Almiron, Bruno Lopes, Alyson L. C. Oliveira, Antonio C. Medeiros, and Alejandro C. Frery. On the numerical accuracy of spreadsheets. Journal of Statistical Software, 34(4):1--29, 4 2010.
Berger2007
Roger L. Berger. Nonstandard operator precedence in excel. Computational Statistics and Data Analysis, 51(6):2788--2791, 2007.
KeelingPavur2007
Kellie B. Keeling and Robert J. Pavur. A comparative study of the reliability of nine statistical software packages. Computational Statistics and Data Analysis, 51(8):3811--3831, 2007.
McCulloughHeiser2008
B.D. McCullough and David A. Heiser. On the accuracy of statistical procedures in microsoft excel 2007. Computational Statistics and Data Analysis, 52(10):4570--4578, 2008.
McCullough2008
B.D. McCullough. Microsoft excel's not the wichmann-hill random number generators. Computational Statistics and Data Analysis, 52(10):4587--4593, 2008.
Nash2008
John C. Nash. Teaching statistics with excel 2007 and other spreadsheets. Computational Statistics and Data Analysis, 52(10):4602--4606, 2008.
Su2008
Yu-Sung Su. It s easy to produce chartjunk using microsoft excel 2007 but hard to make good graphs. Computational Statistics and Data Analysis, 52(10):4594--4601, 2008.

Powyższy spis literatury w formacie BiBTeX.

url | Fri, 03/09/2010 08:41 | tagi: ,
Implementacja prostej procedury ustalania trafności różnicowej

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

BagozziDholakia2006
Bagozzi, R. P., Dholakia, U. M. Open Source Software User Communities: A Study Of Participation In Linux User Groups Management Science 7/52, 2006, p. 1099--1115.
url | Mon, 23/08/2010 11:55 | tagi: , , ,
Obliczanie AVE i CR z pliku out LISRELa

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

FornellLarcker81
Fornell, C. i Larcker, D. F. (1981). Evaluating structural equation models with unobservable variables and measurement error. Journal of Marketing Research, 18(1):39--50.
HairetAl98
Hair, J. F., Black, B., Anderson, R. E., i Tatham, R. L. (1998). Multivariate Data Analysis. Prentice Hall.
url | Sat, 24/07/2010 22:33 | tagi: , , ,
Metody ilościowe w R. Aplikacje ekonomiczne i finansowe

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

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

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

Inne błędy to już pestka:

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

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

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

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

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

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

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

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

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

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

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

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

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

ESS

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

ESS

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

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