Miesiąc: wrzesień 2016

Czasem diabeł tkwi w szczegółach.

Ostatnio trafił do mnie taki ciekawy problem.

Problem:
Mamy dwa zbiory monet (oznaczmy typ 1/typ 2) z końca XIII w., chcemy sprawdzić się czy zostały wyprodukowane tą samą techniką = czy rozkłady ciężarów monet są takie same.
Wiemy, że ich średnie wielkości wagowe są takie same ale mogą różnić się np. skośnością. Jeżeli uznamy, skośności są różne to można przypuszczać, że monety w różnych zbiorach wyprodukowano różnymi technikami.
Średnia wielkość monet jest taka sama, bo (o ile nie dochodziło do oszustwa) z określonej ilości surowca (jednej marki czy jednej grzywny) wiadomo było ile powinno powstać monet. Ale nawet jeżeli średnia jest taka sama to techniki produkcji monet mogą się różnić udziałem lżejszych i cięższych monet.
Monet jednego typu jest 40, drugiego 32 egz.

A więc jak porównać te dwie grupy?

Gdyby obudzić w nocy studenta statystyki i zapytać go o to jak porównywać rozkłady w dwóch grupach, to ten pewnie by tylko wymamrotał przez sen 'test Kołmogorowa Smirnowa’ i przewrócił się na drugi bok.
W programie R ten test dla dwóch prób można wykonać funkcją ks.test{stats} (wiki)

Wikipedia mówi ’The two-sample K–S test is one of the most useful and general nonparametric methods for comparing two samples, as it is sensitive to differences in both location and shape of the empirical cumulative distribution functions of the two samples.’, znaczy że ten test wykrywa różnice i w średnich i skośności i innych nierównościach pomiędzy próbami. Można więc mieć nadzieję, że nawet gdy nie ma różnicy w średnich to różnica w skośnościach będzie wykryta.
Ale…
Test K-S wykrywa różne odstępstwa od równości z różną czułością. Najbardziej czuły jest na przesunięcia ale już na różnice w skośności jest mniej czuły.
Gdyby szukać analogii, można by porównać ten test ze strategią poszukiwania okularów rano. Strategią w której przez 60% czasu szuka się okularów na parapecie, przez 30% na biurku a przez 10% pod łóżkiem. Ta strategia pozwala na znalezienie okularów w każdym z tych trzech miejsc o ile wystarczająco długo szukamy. Ale zgodzicie się, że jeżeli mamy dodatkową wiedzę, że okularów nie ma na parapecie to lepiej wykorzystać inną, bardziej specyficzną strategię.

Podobnie jest z testem K-S. Odpowiednikiem czasu szukania okularów jest tutaj liczba monet potrzebnych do identyfikacji różnicy. Jeżeli nie wiemy co może różnić porównywane grupy to test K-S jest być może nienajgorszym wyborem. Ale jeżeli wiemy gdzie szukać różnic (np. w skośności) to lepiej użyć innych narzędzi.

Jakich?

Niestety [?] teoria oparta o rozkłady gaussowskie się nie przyda, bo te nie są skośne. Ale zwolennicy rozwiązań parametrycznych mogą wykorzystać rozkłady z parametrem położenia, skali i kształtu (GAMLSS) aby testować skośność. Więcej o tych modelach można przeczytać np. w pracy magisterskiej mojej dyplomantki Anny Lis.

Jednak, jeżeli lubicie dźwięk rozpędzającego się wiatraczka chłodzącego procesor, to interesującą alternatywą będzie test oparty o metodę bootstrap lub test permutacyjny.

Technika bootstrap (wiki) pozwala na szacowanie przedziałów ufności dla skośności w obu próbach. Szacując dokładność oszacowania skośności możemy już prosto skonstruować test dla dwóch prób. Schemat jest następujący: 1. z próby pierwszej losujemy ze zwracaniem 40 monet, 2. na otrzymanej próbie wyznaczamy skośność, 3. powtarzamy 1-2 wiele razy, np. 10000 razy, 4. robimy 1-3 dla drugiej próby, 5. porównujemy oszacowania skośności.

Test permutacyjny przemiesza obie grupy (wiki) i sprawdzi, czy po resamplingu połączonej próby i po ponownym podzieleniu wymieszanych obserwacji na dwie podpróby, skośności w nowo utworzonych podpróbach różnią się bardziej niż różniły się w oryginalnych populacjach. Tutaj schemat jest następujący: 1. Łączymy monety w jeden wektor 72 elementowy, 2. permutujemy wektor, 3. dzielimy wektor na dwie podpróby o wielkości 40 i 32, 4. wyznaczamy skośności na podpróbach, 5. wyznaczamy moduł różnicy w skośnościach z punktu 4, 6. powtarzamy 1-5 wiele razy, np. 10000 razy, 7. sprawdzamy czy ta różnica w skośnościach którą widzimy w oryginalnych danych jest większa niż różnice w takich przepermutowanych próbach.

Akurat dla tych konkretnych danych obie te techniki dają bardzo zbliżone wyniki, p-wartości w okolicach 0.04, co przy tak małych próbach można uznać nawet za wyraźny sygnał pochodzenia z różnych populacji.
K-S test nie widział żadnych różnic, raportował p-wartość powyżej 0.8.

Uczę, bo lubię?

screen-shot-2016-09-26-at-21-57-43

,,Uczę, bo lubię?’’ to tytuł wykładu inaugurującego nowy rok akademicki Uniwersytetu Otwartego Uniwersytetu Warszawskiego. Zazwyczaj unikam wykładów inaugurujących, ponieważ sam wykład trwa krócej niż poprzedzające go ceremonie uroczystego powitania połowy sali (tutaj doszła jeszcze ceremonia rozdania nagród). Ale znak zapytania w tytule skusił mnie na kampus główny UW i cieszę się, że na ten wykład trafiłem.

Punktem centralnym wykładu było pytanie, dlaczego robimy to co robimy, czyli dlaczego nauczyciele akademiccy uczą.

Jedni uczą, bo muszą (nie mają stanowisk naukowych, a na naukowo-dydaktycznych trzeba wyrobić pensum), inni uczą, bo chcą. Ale chcieć mogą z różnych powodów, jedni lubią robić show (a tutaj rzesza widzów), inni lubią poczucie władzy (wystawianie ocen), innym może podobać się stabilność zatrudnienia (choć z tą jest różnie).
Są też tacy, którzy uczą z rzemieślniczą starannością robienia porządnie tego co robią. Oraz tacy, którzy po prostu lubią mówić o tym co ich interesuje. Są wreszcie tacy, którzy lubią pracować z młodymi ludźmi, którzy lubią potyczki intelektualne ze studentami. Powodów jest pewnie bez liku, tyle co nauczycieli, niektóre bardziej inne mniej wzniosłe.

Gdy przyjrzeć się w jakiej roli nauczyciele występują w filmach, książkach czy serialach to wizerunek nauczyciela jest różny, od Siłaczki Żeromskiego, przez Alcybiadesa Niziurskiego po postawy z zachodnich filmów takich jak John Keating w Stowarzyszeniu Umarłych Poetów.

Te generyczne wizerunki zazwyczaj mało pasują do nauczycieli akademickich matematyki, którzy wydają się być osobnym plemieniem wśród nauczycieli. Jakim? Bardzo ciekawy wizerunek dydaktyka i naukowca przedstawiony jest w filmie Przestrzenie Banacha (o Stefanie Banachu), który polecił mi znajomy po ww wykładzie.
Warto zobaczyć przed nadchodzącym rokiem akademickim.

screen-shot-2016-09-26-at-23-07-40

Który telefon się wcześniej zepsuje?

Uczestniczyłem ostatnio w ciekawym (od strony metodologicznej) projekcie (patrz też komentarz na dole). A ponieważ przy okazji używałem kilku ciekawych narzędzi, więc poniżej krótko o nich opowiem.

Problem:

Przypuśćmy, że mamy taki problem. Mamy telefony dwóch marek, np A i B. Załóżmy też, że w tych telefonach może zepsuć się bateria lub ekran, innych uszkodzeń nie rozważamy. Jeżeli się cokolwiek zepsuje, telefon jest wymieniany na nowy, nie ma sensu go naprawiać.
Chcemy teraz sprawdzić czy tempo psucia się ekranu w telefonach marki A jest inne niż w telefonach marki B. Podobnie chcemy porównać działanie baterii.
Przypuśćmy, że w badaniu obserwujemy dużo (tysiące) telefonów każdej z marek. Badanie prowadzimy przez 3 lata, ale niektóre telefony uczestniczą w badaniu od roku, inne od 2 lub 3 lat.
Pytanie, jaką metoda porównać te marki?

Podejście brutalne 1:

Najbardziej brutalne podejście to policzenie jaki procent ekranów zepsuł się w marce A a jaki w marce B i porównanie tych procentów np. dokładnym testem Fishera (w R funkcja fisher.test {stats}, więcej o teście tutaj) lub regresją logistyczną.
Ale…
To nie jest dobre podejście, ponieważ telefony nie były obserwowane przez tak samo długi czas. Mogło się zdarzyć, że w jednej grupie pozyskano więcej telefonów na początku, 3 lata temu, więc ich czas obserwacji był dłuższy. A może w innej grupie pozyskano telefony niedawno, np. rok temu, czas obserwacji jest krótszy i dlatego mniej ich się zepsuło.

Podejście brutalne 2:

Skoro czas obserwacji jest różny, to porównajmy liczbę dni do czasu zepsucia się ekranu lub baterii. Możemy to zrobić np. testem Wilcoxona (w R funkcja wilcox.test {stats} więcej o teście np. tutaj).
Ale…
To nie jest dobre podejście, ponieważ nie wszystkie telefony się zepsuły, nie dla wszystkich znam czas do zepsucia. Być może w jednej grupie średni czas do uszkodzenia jest niższy, ale też mniej telefonów się zepsuło. Nie chcemy też odrzucać informacji o tych wszystkich telefonach, które się nie zepsuły.

Mniej brutalne podejście 3:

Skoro mamy różne czasy obserwacji i nie dla wszystkich obserwacji wiemy jaki był czas uszkodzenia ekranu lub baterii, to zastosujmy metody analizy przeżycia. Przedstawmy czas działania ekranu/baterii za pomocą krzywych przeżycia (np. Kaplan-Meier) i sprawdźmy czy marki różnią się używając np. testu log rank (w R funkcja survdiff {survival}, więcej o teście).
Ale…
Obserwujemy dwa zdarzenia, uszkodzenie ekranu lub uszkodzenie baterii. Jeżeli wystąpi jedno uszkodzenie to nie widzimy drugiego (telefon oddajemy i nie wiemy kiedy nastąpiłaby druga awaria). Ale analizując uszkodzenia baterii, nie możemy traktować uszkodzeń ekranu jako obserwacje cenzorowane ponieważ nie musi być to cenzorowanie nieinformatywne (brak spełnionych założeń). Być może uszkodzony ekran wpływa na sposób pracy baterii, a być może wadliwe ekrany i baterie produkowane są w jednej fabryce, itp.

Podejście najmniej brutalne 4:

Skoro różne telefony obserwujemy przez różny czas, nie dla wszystkich znamy całkowity czas działania, czyli mamy cenzorowanie, ale obserwujemy dwa (lub więcej) elementów, które mogą ulec uszkodzeniu, to możemy modelować je łącznie używając modelu konkurujących ryzyk.

W R można to zrobić używając pakietu mstate (artykuł z JSS) czy też od jakiegoś czasu można to zrobić pakietem survival (teraz funkcja Surv{survival} może przyjąć argument type=mstate).
Do porównania krzywych można wykorzystać test Gray’s log-rank test. A w pakiecie ggplotit jest nawet przeciążona funkcja ggplotit() rysujące skumulowane funkcje zdarzeń dla obiektów klasy survfitms lub Cuminc.

btw:
Projekt na którym pracowałem dotyczył biologii molekularnej i analizy przeżycia pacjentów, ale dla przykładu i wnioskowania to akurat bez znaczenia. Podobna metodologię można zastosować do analizy usterek w samochodach czy kliknięć internautów na stronach www.

Dzień Popularyzacji Matematyki – w czwartek 15 IX 2016 – MiNI PW

screen-shot-2016-09-10-at-08-18-10

Już w czwartek (15 IX) na MiNI PW w godzinach 11-17 odbędzie się druga edycja Dnia Popularyzacji Matematyki.

Miejsca na warsztatach pozostały już nieliczne, ale na wielu wykładach jeszcze miejsca są.

Screen Shot 2016-06-08 at 18.42.39

Przejrzeć listę wykładów i warsztatów oraz zapisać można się przez stronę http://dpm.mini.pw.edu.pl/?q=og/1.
Jest wiele ciekawych propozycji.

A na 14:50-15:50 planowana jest MiNI gra terenowa ,,Łamacze Wykresów” w rozszyfrowywanie wykresowych łamigłówek.

Warsaw R Ladies – Warsztaty z R tylko dla kobiet

cut
Kolejne warsztaty w ramach Spotkań Entuzjastów R zaplanowane są na 20 października. To spotkanie będzie wyjątkowe, na ten dzień zaplanowane są dwa warsztaty tylko dla kobiet.
Spotkania R Ladies działają w Londynie, San Francisco i innych miastach. Na celu mają zwiększenie liczby entuzjastek R.

Poinformujcie znajome, które mogłyby być takimi warsztatami zainteresowane.

Planowane są dwa warsztaty (R dla początkujących / ggplot2 w wizualizacji danych), każdy po 20 osób. Miejsca mogą się skończyć szybko, więc warto je szybko zarezerwować.

Więcej informacji o wydarzeniu oraz formularz zgłoszeniowy znajduje się tutaj.

MinechaRts #1 (Minecraft + R + Edgar Anderson’s Iris Data)

How to use R to draw 3D scatterplots in Minecraft? Let’s see.

Minecraft is a game about placing blocks and going on adventures (source). Blocks are usually placed by players but there are add-ons that allow to add/modify/remove blocks through external API.
And this feature is being used in educational materials that show how to use Minecraft to learn Python (or how to use Python to modify Minecraft worlds, see this book for example). You need to master loops to build a pyramid or a cube. And you need to touch some math to build an igloo or a fractal. Find a lot of cool examples by googling ’minecraft python programming’.

So, Python+Minecraft are great, but how to do similar things in R?
You need to do just three things:

  1. Install the Spigot Minecraft Server along with all required dependencies. The detailed instruction how to do this is here.
  2. Create a socket connection to the Minecraft Server port 4711. In R it’s just a single line 
    conn <- socketConnection(host="localhost", port = 4711, blocking=T, server=F, open="r+")
  3. Send building instructions through this connection. For example 
    writeLines("world.setBlocks(0,70,0,10,80,10, 46)", conn)

    will create a cube 11x11x11 made of TNT blocks (id=46 is for TNT, see the full list here) placed between coordinates (0,70,0) (10,80,10). You can add and remove blocks, move players, spawn entities and so on. See a short overview of the server API.

The R code below creates a connection to the minecraft server, builds a flat grassland around the spawning point and plots 3d scatterplot with 150 blocks (surprise surprise, blocks coordinates correspond to Sepal.Length, Sepal.Width, Petal.Length variables from the iris dataset).

#
# A useful function
addBlock <- function(x, y, z, b, conn) {
  writeLines(paste0("world.setBlock(",round(x),",",round(y),",",round(z),",",round(b),")"), conn)
}

#
# Connect to the server (install and run https://www.nostarch.com/download/LTPWM_ch01_update_online.pdf)
conn <- socketConnection(host="localhost", port = 4711, blocking=TRUE, server=FALSE, open="r+")
baseline <- 70

#
# Add two layers of grass and wipe out everything above
writeLines(paste0("world.setBlocks(-80,",baseline-2,",-80,180,",baseline,",180,2)"), conn)
writeLines(paste0("world.setBlocks(-50,",baseline+1,",-50,150,",baseline+50,",150,0)"), conn)

#
# And now add blocks based on iris data
for (i in 1:nrow(iris)) {
  addBlock(10*iris[i,"Sepal.Length"],
           baseline + 2 + 10*(iris[i,"Sepal.Width"] - min(iris[,"Sepal.Width"])),
           10*iris[i,"Petal.Length"],
           c(41,42,45)[as.numeric(iris[i,"Species"])], 
           conn)
}

If you do not like scatterplots try barcharts 😉

Mecze Matematyczne

Podczas wakacji dowiedziałem się o bardzo ciekawej inicjatywie, tj. o Dolnośląskich Meczach Matematycznych. Temat w sam raz na 1. września.

Dolnośląskie Mecze od kilkunastu lat organizuje Małgorzata Mikołajczyk z UWr. Od kilku lat rozgrywane są też mecze w ligach na Pomorzu i Wielkopolsce. W rozgrywkach udział biorą 10 osobowe zespoły, które wystawiają zainteresowane szkoły. Rozgrywki są prowadzone w ligach dla podstawówek, gimnazjów i liceów. W poprzednim roku można było się zgłaszać do końca października. Zespoły spotykają się parami aby rozegrać mecz. Zwycięzcy przechodzą dalej i tak aż do finałów.

Dokładny regulamin meczu znajduje się tutaj. Poniżej opiszę kilka cech, które najbardziej przypadły mi do gustu.

1) Zespoły mają po 10 osób. Zespoły mają godzinę na opracowanie rozwiązań, ale rozwiązania prezentowane są indywidualnie. Jedna osoba może przedstawić tylko jedno rozwiązanie. Wprowadza to do rozgrywki ciekawy element optymalizacji, gdzie zespół musi ustalić kto będzie prezentował które rozwiązanie.

2) W ocenie rozwiązania liczy się nie tylko końcowy wynik, ale (a nawet przede wszystkim) ścieżka wnioskowania i uzasadnienie wyniku. Jeżeli zaprezentowane wnioskowanie ma słabe strony (lub znacznie prostsze rozwiązanie), to drużyna przeciwna może je wytknąć przejmując część punktów. Konieczność argumentacji rozwiązania oraz analiza tej argumentacji ma niesamowitą wartość edukacyjną.

3) W tej rozgrywce matematycznej planowanie i praca zespołowa jest bardzo istotna aby wygrać mecz. To jest najbardziej niesamowity element całej rozgrywki. W zawodach ,,zespołowych” w których brałem kiedyś udział, wkład drużyny nie był taki ważny. Zespół mógł mieć jednego silnego zawodnika by mieć dobre wyniki. W DMM każde zadanie przedstawia inna osoba, więc nawet jeżeli część drużyny jest silniejsza to musi wyjaśnić reszcie cześć rozwiązań.

Przyznacie, ciekawa inicjatywa. Czy znacie coś podobnego w Warszawie?