hakank.blogg: Statistik/data-analys Archives

oktober 28, 2007

Svar på Hakkes frågor angående monovokala ord

Den alltid nyfikne hakke skickade följande mail till Jonas Söderström och mig häromdagen. (Mailet är trivialt redigerat.)

Det går ju inte att kommentera den här gamla godingen :( Ooo, så långa ord
Så jag får mejla kommentaren istället :)

/Håkan (hakke)

Vissa saker man läst sitter liksom kvar i hjärnan som en lös liten
skruv, som ibland vickar till och gör sig påmind. Dit hör det här
inlägget, bland annat för att det var här jag upptäckte bloggarna men
också på grund av dess språkglädje.

Nyss stötte jag på ordet "kulturutbud" i en text. Det är inte särskilt
lättläst, en egenskap det troligen delar med många andra monovokala
ord.

Håkan Kjellerstrand hade ju vänligheten att publicera följande lilla
lista över förekomsten av monovokala ord för var och en av svenskans
vokaler:
Vokal: antal ord ---------------- a: 3137 e: 791 i: 1299 o: 960 u: 625 y: 192 ä: 590 å: 280 ö: 272
Några saker jag blir nyfiken på är:

Efter kom sedan tre (3) frågor som av besvaras (eller åtminstone bekommenteras) var och en i det nedanstående. Notera att ordlistan som användes för att skapa ovanstående förekomstfördelning är äldre i relation till den ordlista som används i dessa svar.

Hakkefråga 1. Vilka är de längsta fem monovokala orden för varje vokal?

Svar fråga 1.
Här är de inte bara de fem utan även de sex längsta ord för respektive monovokal . Ordlängden visas efter ordet. Not: Det kan finnas flera ord med samma ordlängd som den minst långa ordlängden för respektive vokal. Programmet visar då endast att det blir exakt sex stycken ord (och slumpmässighetens underbara men samtidigt starkt underskattade men starkt påverkande hand styrde exakt vilka som visas).

Vokal a: brandalarmapparat: 17 andrahandsmarknad: 17 partssammansatt: 15 branschanpassad: 15 brandhandgranat: 15 tandspacklarnas: 15



Vokal e:

referensfrekvens: 16

frekvensmeterns: 15

telexreglemente: 15

referenselement: 15

meddelelsemedel: 15

pendelfrekvens: 14

Vokal i:

lindningsriktning: 17

tillfriskningstid: 17

lindningsstigning: 17

drivningsriktning: 17

visningsspridning: 17

stigningsriktning: 17

Vokal o:

motorfordonskontroll: 20

kontrollprotokoll: 17

domstolsprotokoll: 17

torvjordskompost: 16

fordonskontroll: 15

kontrollmottolk: 15

Vokal u:

ursprungspunkt: 14

djupbrunnspump: 14

sunhultsbrunns: 14

grundstruktur: 13

sunhultsbrunn: 13

kugghjulspump: 13

Vokal y:

skyddshytt: 10

styckfryst: 10

krymptryck: 10

tryckstyrd: 10

plymprydd: 9

frysskydd: 9

Vokal å:

stålspåntlås: 12

språngstråk: 11

tvångsvård: 10

stånggång: 9

nålsprång: 9

ståltråds: 9

Vokal ä:

rättshjälpsnämnd: 16

kärrsnäppsägg: 13

vändskärsfräs: 13

rännhärdsjärn: 13

skräntärnsägg: 13

ändskärsfräs: 12

Vokal ö: mörkrödglöd: 11 bröstsköld: 10 förströtts: 10 bröstmjölk: 10 slöjdbjörk: 10 bröstböld: 9

Hakkefråga 2. De korta orden intresserar mig inte särskilt mycket. Det skulle vara intressant att se motsvarande sammanställning begränsad till de ord där det finns minst 3, 4 respektive 5 vokaler. Jag gissar att fördelningen mellan ordrikedomen per vokal då också kan komma att förändras något. Kanske blir ledningen för a och i ännu tydligare?

Svar fråga 2

Först kommer den totala fördelningen av antal vokaler per ord som har
minst 2 vokaler för att få en känsla för vad som kommer:

Fördelning av antal vokaler per ord: 2: 9868 3: 3215 4: 539 5: 60 6: 6

Sedan med hakkes föreslagna begränsningar om minst v vokaler.

Vanligaste bokstaven (minst 2 monovokaler): a: 7377 e: 2080 i: 1612 o: 1275 u: 436 ä: 410 ö: 253 å: 196 y: 49



Vanligaste bokstaven (minst 3 monovokaler):

a: 2447

e: 717

i: 350

o: 229

u: 64

ä: 9

å: 2

ö: 2

Vanligaste bokstaven (minst 4 monovokaler):

a: 354

e: 179

i: 38

o: 33

u: 1

Vanligaste bokstaven (minst 5 monovokaler):

a: 42

e: 16

i: 4

o: 4

Vanligaste bokstaven (minst 6 monovokaler): e: 4 a: 1 o: 1

Det finns inga ord i ordlistan med 7 eller fler monovokaler.

För fullständighetens skulle visas här även fördelningen av ordlängden (för ord med minst 2 monovokaler):

3: 22 4: 300 5: 1226 6: 2486 7: 2964 8: 2744 9: 1886 10: 995 11: 517 12: 299 13: 151 14: 62 15: 16 16: 6 17: 13 20: 1

Hakkefråga 3. Undrar om fördelningen ändras över tiden? Sedan listan skapades har det ju kommit en ny version av saol.

Svar (eller snarare kommentar till) fråga 3

Hakke har troligen en poäng att ovanstående beskrivna fördelningar förändras över tiden. Det är dock utanför mitt experimenterande eftersom jag inte använder SAOL utan Den stora svenska ordlistan (eller snarare ett derivat av den ordlista man kan ladda ner här och om vilket kommenteras något mer här nedan).

Svar på anticiperad följdfråga: Nej, jag har inte sparat olika DSSO-versioner för denna typ av jämförelse.

Några vidarekommentarer
För ett antal (cirka 2) månader sedan förnyades monovokaldiskussionen på Blind Höna, i Ooo, så många o:n! Monovokal toppnotering tangerad (där mina findings bygger på samma ordlista som ovanstående analyser). Se även Söndagspyssel där den ursprungliga monovokaldiskussionen fortsatte att diskuteras.

DSSO-listan är samma ordlista som man hittar på http://sv.speling.org/files/ (det görs en omdirigering till DSSO-sajten). Denna ordlista har även används i andra språk-/ordprojekt, t.ex.
* Visa ordklasser (presenteras i Svenska ordklasser samt gissning med hjälp av ordsuffix)
* Consonants Away (presentation i Consonants Away)
* samt ett gäng andra s.k. useless-projekt.

Posted by hakank at 09:50 FM Posted to Språk | Statistik/data-analys | Comments (6)

september 26, 2007

Videoföreläsning om R (del 1, en kort introduktion)

Decision Science News (en av favoritbloggarna kring beslutsteori, dataanalys etc) har nu börjat med videoföreläsningar i dataanalysverktyget R. Se vidare R video tutorial number 1. Det är en kort introduktion men kan vara början på något trevligt.

Som presentatören säger så rekommenderas en bok (eller två) i R (eller det mycet snarlika språket S). En förhoppningsvis fullständig lista över relaterare böcker finns här. Några av mina favoriter är:
* William N. Venables, Brian D. Ripley. S Programming (ISBN: 9780387989662). En introduktion för den som vill gå in på djupet i själva R-språket (eller S).
* William N. Venables, Brian D. Ripley: Modern Applied Statistics with S, 2002 (ISBN: 9780387954578). Det är en relativt avancerad bok men innehåller mycket matnyttig information, många tips och funktioner.

Se även Några videoföreläsningar om statistik, sannolikhet och data mining.

Posted by hakank at 10:31 EM Posted to Statistik/data-analys | Video podcasts | Comments (1)

augusti 21, 2007

Upptäcka kraftiga förändringar i tidsserier (changepoint detection)

I Recension av Keith Devlin, Gary Lorden: The Numbers Behind Numb3rs: Solving Crime with Mathematics nämndes att kapitel 4 When does the Writing First Appear on the Wall? - Changepoint Detection var det som gav mig mest ny information. Så naturligtvis har jag lekt lite med detta.

I Numb3rs-avsnittet Hardball (tredje säsongen) nämner Charlie Epps (matematikern i Numb3rs) Shiryayev-Roberts metod i samband med baseballresultat för att upptäcka förändringar i tidsserier, men den är "slightly more technical to describe" än den metod som E.S. Page skapade varpå man beskriver Pages. (Kapitlet är troligen skrivet av Gary Lorden eftersom han tidigare skrivit om changepoint detection. Och troligen är det E.S. Page paper "A test for a change in a parameter occurring at an unknown point" som avses. Men ingen av detta står uttryckligen i boken. Gnäll^2.)

Först presenteras kortfattat Pages metod, följt av min R-implementation. Hur väl den följer Page har jag faktiskt inte lyckats röna ut. Här används emellanåt "förändringspunkt" i stället för "changepoint"; "breakpoint" används också nedan.

Introduktion och beskrivning av Pages metod
En specifik händelse inträffar eller inte under en dag. Om den inträffar markeras detta i en tidsserie med 1, annars 0.

Bokens exempel är en händelse som inträffat c:a 1 gång per månad, dvs sannolikheten är 1/30. Sannolikheten för att händelsen inte inträffar är således 1-1/30 = 29/30. Detta är normalfallet. För att citera bokens lättsamma exempel på några sådana händelser som kan tänkas inträffa (i genomsnitt) en gång i månaden:

Examples abound - a New Yorker finds a parking place on the street in front of her apartment, a husband actually offers to take out the garbage, a local TV news show doesn't lead off with a natural disaster or violent crime, and so on.

Mer seriöst (och det beskrivs mer sådant i kapitlet) skulle man kunna tänka sig någon form av olycka eller sjukdomsfall som rapporteras varje dag. Notera att vi här endast arbetar med om händelsen inträffar eller inte; det är inte frågan om antalet olyckor etc. (Se i slutet av denna anteckning för hur man arbetar med numeriska värden.)

Vi vill nu upptäcka ("detektera") om/när det blir en drastisk förändring i hur ofta dessa händelser inträffar, t.ex. att de plötsligt inträffar i genomsnitt en gång i veckan i stället, dvs med sannolikheten 1/7. Låt oss kalla detta för extremfallet. Sannolikheten för att en extremhändelse inte inträffar en gång i veckan är således 1 - 1/7 = 6/7.

Falsklarm: Problemet med händelser som inträffar i genomsnitt en gång i månaden är att två händelser som i genomsnitt händer en gång i månaden faktiskt kan inträffa två intilliggande dagar eller klustra sig på andra sätt. Vi vill undvika falsklarm när slumpen gör det den ska (nämligen att slumpa).

Ett metod att minska antalet sådana falsklarm är alltså E.S. Pages metod som boken går igenom med ett exempel. Här har jag tagit beskrivningen från boken rakt av.

Det finns ett numeriskt index S som "spårar" aktiviteten dag för dag.
* Initialt är S = 1.
* För varje dag uppdateras värdet enligt följande:

om händelsen inträffar under dagen: multiplicera S med sannolikheten för att ~~normalfallet~~ extremfallet inträffar delat med sannolikheten för ~~extremfallet~~ normalfallet inträffar, här (1/7) / (1/30) = 30/7 = 4.286

om händelsen inte inträffar under dagen: multiplicera S med sannolikheten för att ~~normalfallet~~ extremfallet inte inträffar delat med sannolikheten för att ~~extremfallet~~ normalfallet inte inträffar, här (6/7) / (29/30) = 0.8867
om S efter en sådan uppdatering blir mindre än 1: återställ S till 1

* Till slut finns ett gränsvärde för när vi ska varna för att en s.k. changepoint har inträffat. Om S överstiger detta gränsvärde anses händelserna nu komma en gång i veckan (extremfallet). I boken har man satt 50 som detta värde. (Det förs en liten diskussion hur olika gränsvärden påverkar hur väl man kan upptäcka förändringarna, men det går jag inte in på här.)

Bokens exempel
Boken använder följande händelsesekvens som exempel, dvs de första två dagarna inträffade händelsen inte, på tredje dagen inträffade den osv.
0 0 1 0 0 0 0 0 1

Från denna serie får man motsvarande S-värden (tänk på att S återställs till 1 om S < 1 och att S är 1 från början). De första två dagarna blir 1 enligt:

S * 0.8867 = 1 * 0.8867 = 0.8867 < 1, varpå S återställs till 1

Här är en liten tabell för att visa utvecklingen av S-värdet.

H S-värde 0 1 0 1 1 4.285714 0 3.800141 0 3.369583 0 2.987808 0 2.649287 0 2.349122 0 2.082965 1 8.926994

Inget av S-värdena är större än gränsvärdet (50), så ingen varning utfärdas för dessa första 10 dagar. Om händelsen inträffar även nästa dag (dag 11) skulle S-värdet vara 38.25854 (8.926994*4.285714) vilket inte heller ger någon varning. Om den även inträffar påföljande dag (dag 12) skulle det en varning utfärdas eftersom 8.926994*4.285714*4.285714 = 163.9652 >= 50.

Poängen med Pages S-serie är att den stiger kraftigt om händelsen inträffar, men avtar långsamt om händelsen inte inträffar. (Eller snarare: i detta exempel är det "kraftigt", det exakta beteendet beror på de specifika sannolikheterna.)

R-funktionen changepoint
Pages metod har implementerat i statistikpaketet R med följande funktion:

changepoint <- function(x,normal=1/30,non.normal=1/7,detect.value=50,plot=T){ occurs <- non.normal/normal # händelsen inträffar not.occurs <- (1-non.normal)/(1-normal) # händelsen inträffar inte S <- 1 vec = NULL # vektor med S-värdena len <- length(x) alerted <- FALSE; cp <- NULL; # changepoint S.cp <- 0; # värdet för S i changepoint for (i in 1:len) { obs <- x[i] # observationen # justera S enligt metoden: S <- ifelse(obs, S*occurs, S*not.occurs) S <- ifelse(S<1,1,S) # reset om mindre än 1



      # varna om S-värdet överstiger gränsvärdet

      if (plot && !alerted && S >= detect.value) {

        cat("Alert! observation #", i, ": ", S, " >= ", detect.value,"\n")

        cp <- i; S.cp <- S; alerted <- TRUE

      }

      vec <- c(vec,S)

   }

   csum <- cumsum(x) # kummulativa summan av händelseserien

   cumsum.val <- 0 # antal observationer

   if (!is.null(cp) && cp > 0) { cumsum.val <- csum[cp] }

   m<-max(csum) # för grafen

   vec.len <- length(vec)

   # plotta tre grafer

   if(plot) {

     op <- par(mfrow=c(3,1)) # plotta grafer

     plot(x,type="l",

         main=paste("Changepoint value: ", cp, " S-value:", round(S.cp,2) ))

     if (!is.null(cp) && cp>0) { lines( rep( cp, 2), 0:1, col="red") }

     plot(csum,type="l", main=paste("Cumulative sum. Num obs.:",cumsum.val) )

     if (!is.null(cp) && cp>0) { lines( rep( cp, m+1), 0:m, col="red")}

     plot(vec,type="l",main=paste("S values. Detect value: ", detect.value) )

     abline( rep(detect.value, vec.len), 0:vec.len,col="red")

     par(op)

   }

# returvärde list(vec=vec,cp=cp,cumsum.val=cumsum.val); }

Funktionen har följande argument, där defaultvärdena är bokens exempel.

x: en vektor av respektive dagars värden om händelsen inträffade eller inte (1 eller 0)
normal:. Sannolikheten för att normalfallet inträffar, default 1/30 för att följa boken
non.normal: Sannolikheten för att extremfallet inträffar, default 1/7
detect.value: gränsvärde för changepoint-varning (default 50)
plot: huruvida man ska plotta respektive skriva ut en varning (när man gör flera simuleringar vill man inte att det plottas och håller på).

Följande värden returneras :

vec: vektor med S-värdena
cp: changepoint, dvs den observation där S-värdet överstiger varningsvärdet
cumsum.val: vilket är summan av antalet händelser som inträffat vid changepoint (egentligen endast för att jag ville kolla lite)

När man kör programmet visas tre grafer om man inte angivit plot=F i funktionen (visas nedan).

1) Först händelsegrafen, med 0 för ej inträffad och 1 för inträffad händelse. Den lodräta röda linjen visar changepointen om sådan finns. "S-value" är S-värdet för changepoint.
2) En graf för det kumulativa antalet inträffade händelser, samt en röd lodrät linje för changepoint. "Num. obs" är det antalet inträffade händelser för changepoint.
3) Den tredje grafen visar S-värdena. "Detect value" är gränsen för S-värdet (dvs då man ska slå larm) och visas som en röd vågrät linke. Notera att det kan vara svårt att se linjen om största S-värde är stort.

Körning av R-koden
Här är bokens exempel:
> x = c(0,0,1,0,0,0,0,0,0,1) > changepoint(x, 1/30, 1/7, 50, T) $vec [1] 1.000000 1.000000 4.285714 3.800141 3.369583 2.987808 2.649287 2.349122 [9] 2.082965 8.926994



$cp

NULL

$cumsum.val [1] 0

Som vi såg tidigare blev det alltså ingen changepoint eftersom S-värdet inte översteg 50. cp är därför NULL och antalet händelser vid changepoint är därför 0.

Graferna ser ut på följande sätt:

Ett mer intressant exempel (slumpat)
Här är ett mer intressant exempel. Det är en längre sekvens som jag helt enkelt slumpar fram med utgångspunkt från exemplets förutsättningar (normalfall resp. extremfall man vill detektera). I ärlighetens namn körde jag funktionen tills det blev en tidsserie med en changepoint.

För de första 30 dagarna inträffar händelsen (1) med sannolikheten 1/30, sedan kommer 30 dagar med sannolikheten 1/7 för 1. Här är slumpfunktionen:

c(sample(0:1,30,prob=c(29/30,1/30),rep=T),sample(0:1,30,prob=c(6/7,1/7),rep=T));

Så kör vi!

> x <-c(0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,1,1,0,0, 0,0,0,0,0,0,0,0,0,0,1,0,1,0,0,1,0,0,0,1,0,0,0,0,0,0) > changepoint(x) Alert! observation # 47 : 168.0480 >= 50 $vec [1] 1.000000 1.000000 1.000000 1.000000 1.000000 1.000000 [7] 1.000000 1.000000 1.000000 1.000000 1.000000 4.285714 [13] 3.800141 3.369583 2.987808 2.649287 2.349122 2.082965 [19] 1.846964 1.637702 1.452150 1.287621 1.141732 4.893139 [25] 4.338744 3.847162 3.411277 3.024778 2.682069 2.378189 [31] 10.192239 43.681023 38.731942 34.343594 30.452447 27.002170 [37] 23.942811 21.230079 18.824700 16.691853 14.800657 13.123736 [43] 11.636810 10.318354 44.221515 39.211196 168.047983 149.008064 [49] 132.125377 566.251614 502.095027 445.207413 394.765194 1691.850832 [55] 1500.163299 1330.194059 1179.482416 1045.846478 927.351557 822.282168



$cp

[1] 47

$cumsum.val [1] 6

Kommentar: Här ser vi alltså att det blev en varning vid observation #47 (dag 47) och då hade händelsen inträffat totalt 6 gånger. Det var på den 16 dagen (31 + 16 = 47) efter att "extremfallet" började gälla (dvs med 1/7 som sannolikhet). Man kan f.ö. notera att just i denna sekvens blev det 2 stycken händelser efter varandra precis i början på "extremfallet".

Grafen för detta exempel:

Not: Om man kör programmet på ovanstående sätt, dvs batchar hela serien på en gång, så förfelar det sitt syfte. Man ska naturligtvis köra programmet varje dag efter man fyllt i dagens "händelsevärde" och kontrollera om det kommer en varning. Programmet fortsätter alltså att köra efter varningen. I stället borde kanske lampor i lustiga färger börja blinka, sirener börja tjuta och den stora valvdörren börjar stängas så att man måste skynda sig in i valvet innan dörren stängs helt med ett tungt och djupt "swoonk". Minsann kommer där en jeep i högsta fart genom stängslet och fyra personer rusar fram till dörren. Kommer de att klara sig? Ojoj, vad spännade det är...

Andra changepoint-metoder
Pages metod (enligt bokens förklaring och exempel i alla fall) hanterar endast binära tidserier (0/1). Det finns även metoder som detekterar changepoints i tidsserier med numeriska värden. Här nämns några R-funktioner.

Ett standardexempel är datamängden Nile: (via ?Nile i R) Measurements of the annual flow of the river Nile at Ashwan 1871-1970. Det intressanta med denna tidsserie är att år 1898 minskade flödet drastiskt eftersom man byggde första Ashwan-dammen. Kan denna förändring uppptäckas med changepoint detection-metoderna (även kallad "change point analysis", "change detection" etc). Svaret är naturligtvis jakande.

R-paketet struccchange
Med functionen breakpoints() i paketet strucchange visas flera olika changepoints (benämns breakpoints och breakdates för respektive observationsnummer och årtal). Brytpunkterna sorteras på hur stora (viktiga) de är. Nedan ser man att observation 28 (år 1898) är den viktigaste förändringspunkten; det hände också något vid observation 83 (år 1953) etc.

> library(strucchange) > data(Nile) > summary(breakpoints(Nile~1))



         Optimal (m+1)-segment partition:

Call:

breakpoints.formula(formula = Nile ~ 1)

Breakpoints at observation number:

m = 1      28

m = 2      28       83

m = 3      28    68 83

m = 4      28 45 68 83

m = 5   15 30 45 68 83

Corresponding to breakdates:

m = 1        1898

m = 2        1898           1953

m = 3        1898      1938 1953

m = 4        1898 1915 1938 1953

m = 5   1885 1900 1915 1938 1953

Fit:

m 0 1 2 3 4 5 RSS 2835156.750 1597457.194 1552923.616 1538096.513 1507888.476 1659993.500 BIC 1318.242 1270.084 1276.467 1284.718 1291.944 1310.765

R paketet bcp
Ennan metod är funktionen bcp() från paketet bcp (Bayesian Change Point) som visar sannolikheterna för att de olika observationerna är en förändringspunkter. Den mest sannolika punkten har markerats med ett gäng utropstecken.

> library(bcp) > Nile.bcp <- bcp(Nile) > summary(Nile.bcp) > Nile.bcp



Probability of a change in mean, Posterior Mean,

and SD for each position:

Probability Mean SD 1 0.016 1085.8 14.15 2 0.012 1085.6 14.06 3 0.014 1085.1 13.61 4 0.018 1085.5 13.80 5 0.018 1085.2 14.49 6 0.026 1084.8 15.67 7 0.046 1079.2 37.90 8 0.038 1089.4 29.73 9 0.062 1093.0 40.96 10 0.058 1082.8 26.11 11 0.028 1075.1 25.60 12 0.016 1072.0 27.79 13 0.016 1071.6 27.32 14 0.010 1071.0 28.02 15 0.022 1070.6 28.80 16 0.026 1069.5 30.92 17 0.022 1070.0 30.39 18 0.044 1067.7 38.25 19 0.088 1071.5 35.46 20 0.062 1084.1 29.90 21 0.054 1091.7 29.46 22 0.038 1098.5 31.98 23 0.026 1102.0 32.26 24 0.024 1104.9 34.27 25 0.028 1105.3 36.10 26 0.136 1102.9 36.36 27 0.184 1071.1 77.60 28 0.618 1029.3 108.61 !!!! 29 0.078 869.8 66.26 30 0.034 856.8 40.52 31 0.042 852.5 28.16 32 0.020 849.5 14.42 33 0.022 850.4 14.29 34 0.008 849.9 14.03 35 0.018 849.6 14.19 36 0.018 850.2 14.20 ...

Den 28:e observationen har alltså högst sannolikhet. Den sista observationen har sannolikhet 1 och som jag ignorerar. En graf plottas också där man tydligt ser förändringen.



> which.max(Nile.bcp$posterior.prob[-c(100)])

[1] 28

> plot(Nile.bcp)

Grafen:

(För övrigt är detta blogganteckning #1200. Yeay!)

Uppdatering
Håkan Karlsson (aka hakke) uppmärksammade uppmärksamt att "normalfallet" och "extremfallet" hade förväxlats i pratbeskrivningen av hur S uppdateras. Tack hakke.

Posted by hakank at 07:39 EM Posted to Statistik/data-analys | Comments (2)

juli 03, 2007

Några videoföreläsningar om statistik, sannolikhet och data mining

Semester och det blir naturligtvis lite "hängmatte", t.ex. i form av videoföreläsningar.

Så här är några föreläsningar att förgylla sommaren med:

Davis Mease: Statistical Aspects of Data Mining
David Mease om "Statistical Aspects of Data Mining" med R (www.r-project.org) och Excel. Mease har detta även som en "riktig" (IRL inte bara URL) kurs med kurssajten www.stats202.com. Det är lugnt tempo som går igenom grunderna i både teori och verktyg.

Statistical Aspects of Data Mining (Stats 202) Day 1
Statistical Aspects of Data Mining (Stats 202) Day 2

I skrivande stund har endast publicerats ovanstående två föreläsningar. Efterkommande borde dyka upp via denna sökning på video.google.com .

Kursbok är Tan, Steinbach, Kumar: Introduction To Data Mining (ej läst).

Peter Donnelly: How juries get fooled by statistics
Peter Donnelly How juries get fooled by statistics. Visar flera av våra ointuitioner som finns inom området, ofta underhållare på ett traditionellt brittiskt-Grants sätt. Finns även på TEDTalks här.
Beskrivning:

Oxford mathematician Peter Donnelly explores the common mistakes we make in interpreting statistics, and the devastating impact these errors can have on the outcome of criminal trials. Statistical uncertainty and randomness, he says, confound many of our assumptions about the world. He shares the case of a British woman wrongly convicted of murdering her two infants -- a verdict reached, in part, by the misuse of statistics.

[För den som vill läsa mer t.ex. om det där mynt-experimentet finns en relativt djupgående redogörelse i Anirban DasGupta: Sequences, Patterns and Coincidences (PDF).]

Brian Brushwood: Scams, Sasquatch, and the Supernatural
Brian Brushwood Scams, Sasquatch, and the Supernatural. Är mest om debunkning av olika typer av pseudo-vetenskap, men innehåller en akt om sannolikheter (sammanträffanden). Underhållande och med högt tempo.

Ever wonder how those guys on TV seem to talk to the dead? What about ESP and psychic surgery? How do street scams and cons work? Want to ... all » know how YOU can trick your friends into believing you have psychic powers?

As a magician, Brian’s wise to all the tricks used by frauds, tricksters, and con artists …and now he’s ready to take YOU to scam school. This is no ordinary lecture: we’re talking hands-on experiments, a live performance of psychic surgery, free giveaways of cash and prizes, and all the secrets TV psychics DON’T want you to know.

Topics covered include: scams, cons, ESP, UFOs, skeptic, skepticism, dowsing, astrology, memory, alternative medicine, psychic surgery, pseudoscience, coincidence, and crop circles.

Relaterat

Persi Diaconis On Coincidence (som jag iofs skrev om för några år sedan men rekommenderar gärna igen).

Talks Hans Rosling: New insights on poverty and life around the world
Talks Hans Rosling: Debunking third-world myths with the best stats you've ever seen

En kul semesterövning är att analysera lottodata från det sydafrikanska lotteriet. Se vidare bloggen Freakonomics: Is the South African Lottery Rigged? A Hands-On Exercise for Bored Blog Readers.

Uppdatering
Äh, jag glömde ju Scott McClouds grafiska novella The Right Number:

The Right Number is a projected three part online graphic novella about math, sex, obsession and phone numbers presented in an unusual zooming format. Click above to read Parts One and Two. (Part Three will hopefully be completed and available before too long.)

(Not: Serien är skriven 2003 så frågan är om det kommer en tredje del. De två första delarna är dock ganska självständiga. Via EconLog.)

(Avslutningsvis - men egentligen helt orelaterat till ovanstående - kan noteras att det för några veckor sedan var 4-årsdagen av denna bloggs födelse, som inträffade helt utan några virtuella tårterier. Referens till tidigare bemärkelseskriverier finns i förra årets 3-årsdag samt en drapa om varför bloggar bör ses som mer än en dagbok.)

Posted by hakank at 10:56 FM Posted to Sammanträffanden | Statistik/data-analys | Video podcasts | Comments (0)

februari 11, 2007

Några länkar om dataanalys, data mining etc 20070211

Statistical Modeling, Causal Inference, and Social Science
The fallacy of the one-sided bet (for example, risk, God, torture, and lottery tickets)
Animated Social Network Visualization

Data Mining Research
Data mining application: terrorism, som även tipsar om en ny bok Data Mining and Predictive Analysis: Intelligence Gathering and Crime Analysis (ISBN: 075067796) skriven av Colleen McCue .

Geeking with Greg
Excellent data mining lecture notes tipsade för länge sedan om kursen CS345, Autumn 2006: Data Mining. Där finns en mängd presentationer av metoder.

Cognitive Daily
Is 17 the "most random" number?
Randomness wrap-up

Nyupptäckta bloggar i ämnet
Data Mining Research tipsade även om en ny blogg Crime Analysis and Data Mining: A meeting place for those interesting in analyzing or solving crimes using Predictive Analytics!, se presentationen.

Data Mining in MATLAB

Posted by hakank at 10:37 FM Posted to Machine learning/data mining | Statistik/data-analys | Comments (0)

februari 04, 2007

Ultimate Research Assistant (Web Edition)

Ultimate Research Assistant (Web Edition) är en trevlig (men tyvärr långsam) applikation för att sammanställa sökresultat på ett mer intelligent sätt än de sökmotorerna. Systemet listar ut fraser som anses vara signifikanta för sökordet och visar representativa sajter

Exempel: sökning på "blogg" ger följande nyckelord.

blogg
att det
jag har
att jag
blogg
Climate Change
det som
för att
har jag
jag har
om att
om det
Om man
som jag

Av ovanstående fraser anses denna blogg vara representativ för flera: "att det", "jag har", "att jag", "om det", "det som", "om man", "som jag". Vilket kan få en att undra...

Vad gäller urvalet av fraser kan man möjligen anta att stackarn blivit förvirrad av innehållet av de 50 första sökresultaten av "blogg" på Yahoo!, eller så är det helt enkelt att - som det heter på julafton - att "Tony förstår inte språket så bra". Intressant nog finns "Climat Change" med, en het potatis i både inom och utanför bloggvärlden.

Sökningar på mer stringenta fackfraser såsom "text mining" (som råkar vara den teknik som Ultimate Research Assistant använder) och "Diaconis" (en husgud som ofta används för testning av sökverktyg) ger betydligt bättre resultat och känns användbart.

Man bör dock notera den brasklapp som står på sajten: It is an experimental proof-of-concept prototype, and should not be used for any official purposes.

Se vidare
Andy Hoskinson:
Creating the Ultimate Research Assistant där tekniken bakom verktyget förklaras.

Samme Hoskinson har även skapat verktyget Keyword Analysis Tool - Advanced Keyword and Keyphrase Extraction Technology for Content Analysis and Search Engine Optimization (SEO).

Wikipedia: Text mining

Tyvärr koms det här även att tänkas på 200 dagar som bl.a. visar förekomsterna av ord på hakank.blogg där ordet "jag" kom på fjärde plats, samt "jag" i bloggen där vidare språkanalyser genomfördes.

(Verktyget funnet via webbserverloggen.)

Posted by hakank at 05:46 EM Posted to Språk | Statistik/data-analys | Comments (4)

januari 26, 2007

The Blog Authorship Corpus

Den som har lust att göra dataanalyser på engelskspråkiga bloggtexter bör kika på The Blog Authorship Corpus:

The Blog Authorship Corpus consists of the collected posts of 19,320 bloggers gathered from blogger.com in August 2004. The corpus incorporates a total of 681,288 posts and over 140 million words - or approximately 35 posts and 7250 words per person.

...

The corpus may be freely used for non-commercial research purposes.

Zipfilen (nedladdningsbar via ovanstående sida) är 305Mb. Datan uppackad är 840 Mb, bestående av 19320 XML-filer med enkelparsrade taggar såsom <Blog>, <date>, <post>.

Kodningen av författarna görs i filnamnet, t.ex.
3802222.female.13.Student.Gemini.xml, dvs löpnummer, kön, ålder, sysselsättning samt stjärntecken (!).

Det görs en analys i J. Schler, M. Koppel, S. Argamon and J. Pennebaker Effects of Age and Gender on Blogging .

Via Data Mining: Text Mining, Visualization and Social Media

(Det vore trevlig med motsvarande corpus av svenska bloggtexter. Och hellre YYYYMMDD-födelsedata än stjärntecken.)

Posted by hakank at 07:32 FM Posted to Blogging | Statistik/data-analys | Comments (3)

januari 15, 2007

Hur mycket kan man sänka världsrekorden? (extremvärdesanalys)

Daily Time-artikeln World record in 100 metres can be lowered, study predicts diskuterar hur mycket nuvarande världsrekord bör kunna sänkas i olika sportgrenar. Analysen är gjord med den statistiska metoden extremvärdesteori (extreme values).

The world record in the 100-metre sprint can be lowered by another half-second before man reaches his limits, according to an expert in the field of extreme values. The men will hit a wall when Asafa Powell’s current record of 9.77 seconds mark is reduced to 9.29, a study by Prof John Einmahl of Germany’s Tilberg university predicted on Thursday.
...

Den undersökning som diskuteras är John Einmahl, Jan R Magnus: Records in athletics through extreme-value theory (PDF). Från inledningen:

How fast can humans or cheetahs run? What is the total (insured) loss of hurricane Katrina? Is a specific runway at JFK airport long enough for a safe landing? How high should the Dutch dikes be in order to protect us against high water levels? How long can we live? These questions all relate to extremes. In this paper we shall be interested in two questions on extremes relating to world records in athletics. The first question is: what is the ultimate world record in a specific athletics event
(such as the 100m for men or the high jump for women), given today's state of the art? Our second question is: how `good' is a current athletics world record? An answer to the second question will enable us to compare the quality of world records in diffrent athletics events.

På sidan 12 finns en fin tabell Ultimate world records, their precisions, and the current world
records återspeglad här nedan. Dvs det nuvarande världsrekordet på 100m (män) är 9.77. Enligt analysen ska man kunna pressa ner detta till 9.29.

	Men			Women
Event	endpoint	precision	world record	endpoint	precision	world record
100m	9.29	0.39	9.77	10.11	0.40	10.49
110/100m h.	12.38	0.35	12.88	11.98	0.19	12.21
200m	18.63	0.88	19.32	20.75	0.57	21.34
400m	-	-	43.18	45.79	1.83	47.60
800m	1:39.65	1.44	1:41.11	1:52.28	1.39	1:53.28
1500m	3:22.63	3.31	3:26.00	3:48.33	2.78	3:50.46
10,000m	-	-	26:17.53	-	-	29:31.78
marathon	2:04:06	57	2:04:55	2:06:35	10:05	2:15:25
shot put	24.80	1.25	23.12	23.70	0.86	22.63
javelin throw	106.50	10.30	98.48	72.50	2.99	71.70
discus throw	77.00	2.85	74.08	85.00	8.10	76.80
long jump	-	-	8.95	-	-	7.52
high jump	2.50	0.05	2.45	2.15	0.05	2.09

Mer om detta
Econometrists Tilburg University Calculate the Ultimate Athletics Records

Den data som används i undersökningen finns tillgänglig från:
* Track Field Athlectics (Hans-Erik Pettersson)
* International Association of Athletics Federations (IAAF)

Mer om den extremvärdesanalys samt introducerande litteraturreferenser finns i Extremvärdesanalys av webbesök.

Posted by hakank at 05:39 EM Posted to Sport, idrott, hälsa | Statistik/data-analys | Comments (1)

oktober 29, 2006

Dagslänkar 2006-10-29 - tema lotterier

Maths theory bags lotto jackpot, via Minding the Planet

CBC News: Lottery retailers enjoying luck of the draw: Fifth Estate probe, via Sabermetric Research

Och varför inte de något relaterade
Lite om konsekutiva talpar i Lotto (och dess varianter).
Littlewood's Law of Miracles -The law of truly large numbers

Posted by hakank at 06:47 FM Posted to Statistik/data-analys | Comments (0)

rcompletion: Tab-komplettering i R

I statistik- och dataanalyspaketet R har det funnits tab-kompletteringar (för kommandoradsvarianten) men endast för filnamn. Länge har jag saknat att kunna tab-komplettera funktionsnamn à la Emacs eller Unix-skal.

I morse - i loggen över icke automatinstallerade paket - hittades paketet rcompletion som gör precis detta. Det är inte bara variabelnamn och funktionsnamn som kompletteras, utan man kan även få fram parametrarna i en funktion. Det sistnämnda känns mycket användbart.

En första not
rcompletion funkar endast för kommandoradsversionen av R och endast på system som har libreadline. Se nedan hur man installerar paketet, jag hade nämligen lite problem med detta.

Några exempel

Vi börjar med att ladda in paketet tseries för att testa lite
> library(tseries)

Så börjar vi vår tab-session, där <TAB> alltså betyder tab-tangentens nedtryckande.

> ari<TAB>

som ger följande förslag:

arima( arima0( arima0.diag( arima.errors( arima.sim(

Varpå vi sedan tabbar inom en funktion, dvs inklusive första parentesen och får listan över funktionens parametrar.
> arima(<TAB> fixed = kappa = n.cond = seasonal = xreg = include.mean = method = optim.control = transform.pars = init = model = order = x =

Nifty!

Mixtrande med installationen av rcompletion (version 0.8) hos mig
Tyvärr så var jag tvungen att mixtra lite för att kunna installera paketet. Anledningen till att det inte gick in vid automat-installationen (via R:s normala install.packages()) var att konfigurationen inte tyckte att min dator har libreadline.so, vilket ju är väldigt felaktigt.

Varpå följande gjordes:

1) ladda ner paketfilen (via paketets sida sida)

2) manuellt packa upp paketfilen i lämplig katalog (tar zxcv paketnamnet)

3) i katalogen som skapades finns filen configure. På ungefär rad 2287 (för version 0.8) finns dessa rader som utan pardon kommenterades bort enligt följande:
# if test $ac_cv_lib_readline_rl_completion_matches = yes; then



             PKG_CFLAGS="${PKG_CFLAGS} -DHAVE_LIBREADLINE"

             PKG_LIBS="${PKG_LIBS} -lreadline"

#else # # { { echo "$as_me:$LINENO: error: readline is required but was n\ ot found" >&5 #echo "$as_me: error: readline is required but was not found" >&2;} # { (exit 1); exit 1; }; } # exit 1 # #fi

dvs raderna med PKG_CFLAGS och PKG_LIBS ska vara kvar.

4) efter detta installerar man paketet manuellt (eventuellt som root om R är så installerat).
R CMD INSTALL paketnamn

5) I R laddar man paketet som vanligt med library(rcompletion)

Paketförfattaren är medveten om att det finns problem och skriver mer här.

Uppdatering torsdag 11 nov 2006
Installationsproblemet ser ut att ha löst sig. Paketet installeras nu helt sömlöst.

Posted by hakank at 06:22 FM Posted to Statistik/data-analys

juni 08, 2006

Vilken boll är rundast? Slumpmässighet i sporter

I P1 pratades precis om att fotboll är den mest slumpmässiga sporten. Matematikern som intervjuades, Torbjörn Lundh, har skrivit Which ball is the roundest? – a suggested tournament stability index (PDF).

Abstract:

All sports have components of randomness that team not to win every game. According to many spectators of the charm when following a competition or a match. less of this unpredictability? We suggest here a general stability index, which could measure this randomness factor and different sports. As an illustration we use exemplify squash, and soccer. Results will also be given on tournaments football, ice-hockey, and handball. Furthermore, we will optimization questions that turned up on the way.

(Funnen via sidan Vetenskapsfestival 2006)

Uppdatering
Intervjun kan (förhoppningsvis) ålyssnas via följande djuplänk.

Malin på Vetenskapsnytt skrev i vintras (onsdag, januari 04, 2006) om en liknande undersökning: Fotboll är den mest spännande sporten, säger forskare, där följande paper omnämns:
E. Ben-Naim, F. Vazquez, S. Redner What is the most competitive sport? (arXiv.org)

We present an extensive statistical analysis of the results of all sports competitions in five major sports leagues in England and the United States. We characterize the parity among teams by the variance in the winning fraction from season-end standings data and quantify the predictability of games by the frequency of upsets from game results data. We introduce a mathematical model in which the underdog team wins with a fixed upset probability. This model quantitatively relates the parity among teams with the predictability of the games, and it can be used to estimate the upset frequency from standings data. We propose the likelihood of upsets as a measure of competitiveness.

Andra bloggar om: matematik, sport.

Posted by hakank at 07:42 FM Posted to Matematik | Sport, idrott, hälsa | Statistik/data-analys | Comments (3)

januari 23, 2006

Något om prefixträd sorterade på lite olika sätt samt komprimering

Av en privat (och personlig) orsak behövdes ett program för att skapa unika prefix av ord i en lista, och sedan att presentera dessa prefix som ett träd. Anledningen var ungefär att korta ner en lista av ord till kortare ord (dess prefix).

Ett exempel på ett sådant prefixträd är för de svenska månadsnamnen:
a: ap: april au: augusti d: december f: februari j: ja: januari ju: jul: juli jun: juni m: ma: maj: maj (!!!) mar: mars n: november o: oktober s: september

Här ser vi t.ex. att månadsnamnen december, februari och november kan komprimeras till en enda bokstav, medan april och augusti båda börjar på bokstaven a så vi måste här ta till en extra bokstav för att skilja dem åt. Man kan också notera att maj inte blev något prefix alls, utan här behövs hela ordet eftersom mars, den rackaren, ligger och stör. Sådana fall markeras med !!!!

Program Prefix Trees
Det har naturligtvis skapats ett program för att kunna leka vidare med sådana här ord-listor, och som - möjligen något missledande - kallas för Prefix trees.

Programkörningen för ovanstående svenska månadsnamen finns ungefär här.

Komprimering och annat skoj
Även om det primära syftet med programmet var att skapa unika prefix och prefixträd så var det svårt att låta bli att studera vissa saker i mer detalj.

Om man kör programmet för de svenska månadsnamnet kommer det, förutom det fina prefixträdet, även lite statistik (längst ner på sidan):
... Original total length of words: 74 Total length of prefixes: 23 Mean prefix length: 1.92 Compression factor (original length/prefix length): 3.217

Man kan bl.a. se att det medellängden av de skapade prefixen är 1.92, vilket ska jämföras med medellängden av de ursprungliga orden på 6.17. Komprimeringsfaktorn är 3.217, vilket räknas ut med formeln
(sammanlagda längden av original orden) / (sammanlagda längden av prefixerna)

Dvs: genom att använda de unika prefixen istället för de hela orden sparar man en hel del, t.ex. träd (sic!) om man skulle använda papper för att skriva dem. Det kanske blir svårt att läsa det, men man skulle kunna spara både tid och pengar på detta. Nedan kommer vi att studera hur man kan spara något mer tid/papper/pengar men i gengäld blir det i stort sett oläsbart...

Frekvenssortering av bokstäver
En finess med programmet är att man kan studera lite olika varianter av representation av orden innan det skapas prefix. T.ex.

bokstäverna i ordet sorteras innan man gör prefixen, valet Sort (plain)
bokstäverna i ordet sorteras i omvänd ordning, valet Reverse. Not: om man både sorterar och reverse, blir det i omvänd sorteringsordning, dvs "ö" kommer före "ä" etc.

Sedan den mest intressanta varianten: bokstavsfrekvenssortering (valet Sort by letter frequency, som innebär att innan man skapar prefixen sorterar man ordets bokstäver i ordning av hur vanliga/sällsynta bokstäverna är bland samtliga ords bokstäver. T.ex. så är det troligt att bokstaven "z" är sällsynt vilket gör att ett ord med "z" kommer att prefixas som z + eventuellt någon annan bokstav.

Om vi fortsätter med månadsnamnsexemplet får man med bokstavsfrekvenssortering följande prefixträd:
d: december f: februari g: augusti j: maj k: oktober l: juli n: nj: nju: njui: juni njuia: januari p: pl: april ps: september s: mars v: november .... Original total length of words: 74 Total length of prefixes: 21 Mean original length: 6.17 Mean prefix length: 1.75 Compression factor (original length/prefix length): 3.524

Om man jämför med det "vanliga" prefixträdet så är komprimeringsfaktorn något lägre för bokstavsfrekvenssorteringsvarianten (3.524 jämfört med 3.217).

Denna skillnad tenderar att vara beständigt: Detta innebär att om man använder bokstavsfrekvenser istället blir det något bättre komprimering. Nackdelen är naturligtvis att "förkortningen" (prefixet) är fullständigt obegriplig. Se även en mer systematisk genomgång nedan.

Några noter kring detta
Möjligen skulle man här kunna göra en analogi med Huffman-kodning (ref till en.wikiedia) som skapar binära koder efter bokstävernas frekvens, där den kortaste koden tilldeleas den vanligaste bokstaven etc.

En mer avancerad variant av prefixträd skulle vara att analysera prefixträdet för att komprimera det ännu mer.

Det går naturligtvis också att göra det väldigt enkelt för sig och koda varje ord till ett eller flera fullständigt slumpmässiga tecken alltefter ordens/bokstävernas frekvens i texten. Men det är en övning som lämnas åt vidare öden (eller övning åt läsaren/skribenten).

Den ursprungliga poängen med prefixträden var att de skulle vara lätta att skapa och (relativt) lätt att uttyda när man ser dem. När man använder sorteringar av olika slag förloras detta syfte bort i skymningen varvid endast nästa dags solstrålar är glada att se det.

Liten undersökning
Det gjordes också en undersökning hur antalet olika ord i listan påverkar komprimeringsfaktorn respektive medelprefixlängden. Detta gjordes för ett olika antal slupmässiga svenska ord (ur en ordlista på nästan 400000 ord).

Det visar sig - inte speciellt förvånande - att ju fler slumpmässiga ord listan innehåller, desto sämre blir komprimeringsfaktorn och desto längre blir medellängden på prefixen. Som variansen visar är det relativt stabila värden förutom för den första körningen (10 ord i listan). Av tidsskäl kördes endast 10 gånger med samma liststorlek.

Prefix utan någon bokstavssortering

Antal ord	Komprimeringsgrad	Komprimeringsgrad varians	Prefix medellängd	Prefix medellängd, varians
10	7.80730	4.67526	1.47000	0.12456
100	4.14080	0.02874	2.63700	0.00842
200	3.54860	0.03530	3.11500	0.01812
500	3.04260	0.00201	3.60800	0.00775
1000	2.70590	0.00199	4.10200	0.00368
2000	2.39940	0.00078	4.60200	0.00717
3000	2.23620	0.00063	4.96000	0.00207
5000	2.06810	0.00030	5.36400	0.00158
10000	1.84500	0.00009	5.99800	0.00057
20000	1.65400	0.00002	6.68800	0.00040
50000	1.44040	0.00001	7.68000	0.00031
100000	1.31520	0.00000	8.41300	0.00020

Prefix med frekvensbokstavssortering

Antal ord	Komprimeringsgrad	Komprimeringsgrad varians	Prefix medellängd	Prefix medellängd, varians
10	9.68830	4.15306	1.15000	0.02500
100	4.24470	0.02675	2.56800	0.01280
200	3.81910	0.02128	2.94800	0.00544
500	3.26520	0.00193	3.40900	0.00283
1000	2.89340	0.00351	3.83500	0.00372
2000	2.61810	0.00032	4.24200	0.00044
3000	2.47130	0.00082	4.48900	0.00159
5000	2.28780	0.00018	4.82100	0.00054
10000	2.06480	0.00009	5.36200	0.00042
20000	1.87310	0.00003	5.91000	0.00018
50000	1.65010	0.00001	6.70800	0.00008
100000	1.50710	0.00000	7.34400	0.00014

Andra exempel
Jag har även kört programmet på några andra exempel, t.ex. 1000 vanligaste svenska orden, svenska förnamnen samt svenska efternamnen (plus "Kjellerstrand") från Språkbanken:

1000 vanligaste svenska orden (inklusive siffror, årtal och delar av förkortningar)
1000 vanligaste svenska orden, sorterade enligt bokstavsfrekvens
1000 vanligaste svenska förnamnen
1000 vanligaste förnamnen, sorterade enligt bokstavsfrekvens
1000 vanligaste svenska efternamnen (plus "kjellerstrand")
1000 vanligaste svenska efternamnen, sorterade enligt bokstavsfrekvens (plus "kjellerstrand")
1000 slumpmässiga svenska ord.
500 slumpade svenska ord, utan sortering samt med bokstavsfrekvenssortering.

Se även
En möjlig variant att göra en komprimering eller på annat sätt förkorta saker och ting är att använda reguljära uttryck, t.ex. som i programmet MakeRegex,.

Posted by hakank at 10:51 EM Posted to Program | Språk | Statistik/data-analys | Comments (2)

oktober 29, 2005

Kritik kring power law-forskning

Några korta länkningar till en intressant kritik av power law-forskningen och speciellt av Albert-László Barabasi (som jag skrivit om en del tidigare).

Joao Gama Oliveira och Albert-László Barabási skrev nyligen Darwin and Einstein correspondence patterns (PDF, Complementary Materials). New Scientist skrev sedan om detta i Email and letter writing share fundamental pattern. Jämför med Barabasis The origin of bursts and heavy tails in human dynamics (PDF, dess Supplement) som undersöker hur mail skickas och besvaras mellan en flera personer.)

Här är några aktuella och tidigare kritiska röster kring power law-forskningen. "Huvudpersonerna" är två statistiska fördelningar:
Log normal
Power law
"Normal" som Shalizi skriver om i citatet nedan är Normal distribution.

Cosma Shalizi: Gauss Is Not Mocked:

[T]he apparent power law is merely an artifact of a bad analysis of the data, which which is immensely better described by a log-normal distribution.
...

Log-normals are very common, for the same reasons that normals are. Unlike normals, they are very easy to mistake for power law distributions, especially if your knowledge of statistics is as limited as most theoretical physicists'. (The distribution of links to weblogs, for instance, is much better fit by a log-normal than a power law, as we've seen.)

Cosma Shalizi: Speaking Truth to Power About Weblogs, or, How Not to Draw a Straight Line.

Daniel B. Stouffer, R. Dean Malmgren, Luis A. N. Amaral: Comment on Barabasi, Nature 435, 207 (2005)

Michel L. Goldstein, Steven A. Morris, Gary G. Yen: Problems with Fitting to the Power-Law Distribution

Aaron Clauset: Links, links, links.

Pharyngula: λ >> µ, or why I haven't answered you yet

Geomblog: Darwin's and Einstein's (e)mail correspondence rates, or a rumination on power laws.

Michael Mitzenmacher har skrivit en översikt över power-law-forskningen: A Brief History of Generative Models for Power Law and Lognormal Distributions (PDF, publicerad i Internet Mathematics, volume 1); en tidigare version PS), och berättar i en videoföreläsning sina egna åsikter om denna forskning: New Directions for Power Law Research, där själva videoströmningen (QuickTime) finns här. (Från Models of Real-World Random Networks med en massa andra intressanta videos.)

Se även B. Conrad and M. Mitzenmacher Power Laws for Monkeys Typing Randomly: The Case of Unequal Probabilities (PDF).

Posted by hakank at 10:03 FM Posted to Social Network Analysis/Complex Networks | Statistik/data-analys | Comments (0)

augusti 29, 2005

Kortsamlarproblemet (coupon collector problem) i mindre och isolerade nätverk med fullständigt utbyte av kort

Introduktion

När jag var liten samlade jag en hel del på samlarkort, t.ex. flaggor, ishockey-/fotbollsspelare etc. Det var alltså sådana kort där man köpte osett, ett eller flera kort i en påse (med eventuellt medföljande tuggummi eller annat godis) och sedan antingen köpte man nya eller bytte man till sig de kort man inte hade. Den stora poängen var naturligtvis att få hela serien komplett. (Det fanns säkert en fin social poäng i att träffas och byta sådana kort, men sådana finesser uppfattade man nog inte i denna ålder, typ kring 10+ år.)

Speciellt för flaggorna kommer jag ihåg att det fanns vissa kort som var väldigt svåra att få tag på (var det Mexiko med nummer 104, eller kanske var det Kenya?). Om jag nu kommer ihåg detta korrekt var det några få lyckostar som fått detta kort och vi andra dräglade djupt över dessa rara kort.

Kortsamlarproblemet (coupon collector problem)
Det finns ett klassiskt problem inom statistiken/sannolikhetsläran som handlar om detta: kortsamlarproblemet (coupon collector problem, dvs kupongsamlarproblemet) där man studerar hur många kort man behöver köpa i genomsnitt (eller för en viss säkerhetsgrad) på att få samtliga kort i en sådan samlarserie. Det visar sig att det är ointuitivt många: För att få tag i samtliga kort i en serie av 100 kort behöver man i medeltal köpa 519 kort. (Enkel formel för detta : 100*sum(1/(1:100)) ~ 519, där 100 i formeln är just antal unika kort i samlarserien.)

Men detta antal fluktuerar naturligtvis: för att vara säker till 95% att få tag i samtliga kort krävs 791 inköpta kort. Se även min Simuleringssida för lite mer om detta; sök på "samlarkortproblemet ".

Man kan notera att ovanstående beräkningar endast gäller för en person som ensam inhandlar kort. Om man tillåter byte av kort är det en annan sak, och det är precis det som det följande kommer att handla om.

Samlarkortproblemet i ett isolerat nätverk med fullständigt utbyte
För ett tag sedan hade jag en diskussion med en kompis om detta. Denne - som vi kan kalla för E - hävdade att företagen troligen inte tillverkar lika många av varje kort, just för att få oss att köpa fler kort för att få tag i samtliga.

Denna konspirationsteori är möjligen korrekt, men jag började undra om denna effekt kunde uppstå även om företaget tillverkade lika många exemplar av alla kort, och att det vi råkat ut för var en "klustereffekt", dvs att de kort som vi i vårt lilla nätverk var en slumpmässigt urval med de traditionella egenskaperna hos sådana urval.

E antog att företagen som skapar korten alltså tillverkar vissa kort i färre upplagor, säg att de färsta (sic!) är en tiondel av de som tillverkas mest. Min utgångspunkt var att denna effekt skulle kunna vara subjektiv enligt "regeln" att om man slumpmässigt drar ett antal olikfärgade kulor från en påse med lika fördelning så är det mycket stor sannolikhet att något antal är färre än de övriga och att det finns en färg som dras fler gånger.

Detta är något man kan testa med en simulering. Vilket har gjorts.

Sammanfattning
Som vi kommer att se uppstår faktiskt denna typ av effekt (E-effekt) om det totalt köps in ett mindre antal samlarkort. För 700 inköpta kort blir förhållandet mellan minst antal kort och flest antal kort rätt exakt 1/10.

Här nedan kommer även andra saker att visas, t.ex. denna effekt hos Lotto-dragna nummer, där samma förhållande är 1/2, vilket även uppstår i en simulering. Kanske inte så remarkabelt men intressant.

Lokalt nätverk med fullständigt utbyte av samlarkort
Först en sak om nätverket. Den modell som vi här arbetar med innebär att det finns en mindre grupp, säg ett 10-tal personer, som har fritt utbyte av samlarkort inom gruppen men inte har något som helst kort-utbyte med några externa personer eller andra nätverk. Denna fullständiga isolering av inköpen gör modellen enkel att simulera.

(Man skulle även kunna studera komplexa nätverks-effekter kring detta, dvs att någon i Husieskolans krets känner någon i Stockholm och byter till sig de i Husie-kretsens ovanliga kort. Ju fler sådana svaga länkar mellan olika (inte helt isolerade) nätverk det finns, desto större blir naturligtvis den totala kretsen vilket leder till att denna lokala nätverkseffekt vi talar om här blir mindre.)

Modell och simulering
Tänk nu följande som översättning av samlarkortproblemet i ett lokalt nätverk till slumpmässig dragning av heltal i ett intervall.

En grundläggande förutsättning är att det faktiskt skapas exakt lika många kort, och att man globalt har lika stor chans att köpa alla kort, men vi studerar endast en mindre krets (nätverk) av personer som får tillgång till ett urval av dessa, t.ex. där tobakaffärens och den där kioskens urval.

Låt oss för enkelhetens skull anta att samlarkortserien har 100 olika kort. Dessa 100 olika korten motsvaras i simuleringen då av de första 100 positiva heltalen (1..100 alltså). Ett inköp av ett kort motsvarar att man slumpmässigt drar ett tal inom detta intervall, där alla tal har lika stor chans att dras. Detta görs s.a.s. med återläggning, så att ett tal har möjlighet att dragas igen.

Man drar slumpmässigt n stycken sådana tal, vilket alltså motsvarar inköp av n stycken samlarkort.

Efter det beslutade antalet "kort" (tal) är dragna, kontrolleras fördelningen bland dessa kort. Här jämför man antalet kort som kommit upp med minst antal och de som kommit med flest antal. Om ss är listan över de dragna talens fördelning är det sökta värdet alltså: min(ss)/max(ss), det som nedan kallas för min/max-förhållandet.

Vi behöver inte förutsätta något om hur stor gruppen är eller hur många kort en specifik deltagare införskaffar. Vi kommer senare se hur många fullständiga serier ett visst antal inköpta kort tenderar att generera, och man kan alltså anpassa antalet kort som (i medeltal) behöver köpas in för att få fullständiga serier.

Frågeställningen och hypotesen
Låt oss alltså se hur fördelningen mellan minimi-antalet kort och maximi-antalet kort utvecklas med det totala antal kort som inköpes.

Frågeställningen: Vad är förhållandet mellan min/max om det totalt köps in num.bought stycken kort av num.cards möjliga?

Hypotes: Om man slumpmässigt kommer fram till sådana min/max-förhållande på säg 1/10 eller något liknande har vi inget statistiskt stöd för E:s konspirationsteori.

Här följer R-koden som använts för simuleringen, som helt enkelt drar num.bought stycken tal ur intervallet 1..num.cards och returnerar fördelningen av de dragna talen (dvs hur många 1:or som dragits, hur många 2:or etc). Man kan här också notera hur enkelt denna simulering är. Det beror på att vi gjort det enkelt för oss och förutsatt att det fullständigt och fritt utbyte i detta isolerade nätverk.

collector.group.sim <- function(num.cards=100, num.bought=10000) { ss <- rep(0,num.cards); for (i in 1:num.bought) { s <- sample(1:num.cards,1); ss[s] <- ss[s] + 1; } ss }

Exempel:
Om vi nu antar att det är 100 olika kort och 10 personer som samlar, samt att dessa vardera köper respektive 200 kort, inalles 10*200 = 2000 kort. En sådan simulering ser ut så här (sorterad lista eftersom vi inte behöver bry oss om de specifika talen/korten som dragits):

> ss<-collector.group.sim(100, 10*200);sort(ss);min(ss);max(ss);min(ss)/max(ss) [1] 10 11 12 12 12 13 14 14 14 15 15 15 16 16 16 16 16 16 17 17 17 17 17 17 17 [26] 17 18 18 18 18 18 18 18 18 19 19 19 19 19 19 19 19 19 19 19 19 20 20 20 20 [51] 20 20 20 20 20 20 20 20 20 20 20 20 21 21 21 21 21 21 21 21 21 22 22 22 22 [76] 22 22 23 23 23 23 24 24 25 25 25 25 25 25 26 26 26 27 27 27 27 28 28 33 33 [1] 10 [1] 33 [1] 0.3030303 > 10*200 [1] 2000

Här kan vi alltså se att det minsta antalet dragna tal/kort är 10, dvs min(ss), och max är 33. Förhållandet mellan dessa min(ss)/max(ss) är 10/33 ~ 0.3. Redan enna enkla simulering (detta nätverk) ger en känsla av att vissa kort är ovanligare än varandra, och förhållandet mellan de ovanligaste och de vanligaste är 10 på 33.

Personligen tycker jag detta är emot intuitionen: efter en slumpmässig dragning har vi ett min/max-förhållande på så stort som 0.3.

Om vi istället säger att man tillsammans endast köper 519 kort (dvs så många som det i genomsnitt krävs för att en person ska få en komplett serie) får vi några intressanta resultat. Notera att även här är det endast en enda körning, för att få en känsla för datan.

> ss<-collector.group.sim(100, 519);sort(ss);min(ss);max(ss);min(ss)/max(ss) [1] 1 1 2 2 2 2 2 2 2 2 3 3 3 3 3 3 3 3 3 3 3 3 4 4 4 [26] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 5 5 5 5 5 5 5 5 5 [51] 5 5 5 5 5 5 5 5 5 5 5 5 5 6 6 6 6 6 6 6 6 6 6 6 6 [76] 7 7 7 7 7 7 7 7 7 7 7 8 8 8 8 8 8 8 9 9 10 10 11 12 13 [1] 1 [1] 13 [1] 0.07692308

Här kan till exempel följande noteras:

* min/max-förhållandet är värre än E:s föreslags 1/10, nämligen 1/13 (0.077)

* det minsta talet är 1 vilket innebär att det finns kort för en fullständig serie. (Fullständiga serier kommer att studeras mer här nedan.)

Vi gör nu en mer statistisk korrekt simulering med 100 sådana här simuleringar och tar medelvärdet för min/max-förhållandet:
> mean(replicate(1000,{ss<-collector.group.sim(100, 519);min(ss)/max(ss)})) [1] 0.05141527

min/max-talet 0.05 (~ 1/20) stödjer ovanstående känsla att det är värre min/max-förhålllande än E:s 1/10.

Men det är kanske inte så realistiskt att en grupp tillsammans köper endast 519 kort? Låt oss nu titta på hur min/max-förhållandet påverkas av antalet inköpta kort.

Tabell över (medelvärde) av min/max-förhållandet för olika antal kortinköp
Om vi nu studerar medelvärdet av min/max för olika antal kortinköp får man följande tabell.

Från 100 inköp till 3000
> sapply(1:30, function(n) rbind(n*100, mean(replicate(100,{ss<-collector.group.sim(100, n*100);min(ss)/max(ss)})) )) [,1] [,2] [,3] [,4] [,5] [,6] [,7] [1,] 100 200 300 400.00000000 500.00000000 600.00000000 700.0000000 [2,] 0 0 0 0.01154701 0.04447036 0.07848971 0.1029102



            [,8]        [,9]        [,10]        [,11]        [,12]

[1,] 800.0000000 900.0000000 1000.0000000 1100.0000000 1200.0000000

[2,]   0.1298539   0.1475774    0.1695608    0.1864951    0.1996583

            [,13]        [,14]        [,15]        [,16]        [,17]

[1,] 1300.0000000 1400.0000000 1500.0000000 1600.0000000 1700.0000000

[2,]    0.2132364    0.2434604    0.2501553    0.2677531    0.2798811

            [,18]        [,19]        [,20]        [,21]        [,22]

[1,] 1800.0000000 1900.0000000 2000.0000000 2100.0000000 2200.0000000

[2,]    0.2819174    0.2909844    0.3101849    0.3119899    0.3260513

            [,23]        [,24]        [,25]        [,26]        [,27]

[1,] 2300.0000000 2400.0000000 2500.0000000 2600.0000000 2700.0000000

[2,]    0.3373065    0.3520848    0.3504074    0.3622678    0.3721507

[,28] [,29] [,30] [1,] 2800.0000000 2900.0000000 3000.00000 [2,] 0.3841707 0.3844286 0.38901

Här kan man se att gränsen för 1/10 i förhållande mellan min och max är cirka 700 inköp. Detta motsvarar en liten grupp på 7 personer som köper vardera 100 kort, eller 4 personer som köper 175 kort osv. Jag tycker att det låter rätt lite, men å andra sidan var E:s 1/10 taget ur luften.

Om vi tar lite större steg: mellan 1000 till 30000 inköp:
> sapply(1:30, function(n) rbind(n*1000, mean(replicate(100,{ss<-collector.group.sim(100, n*1000);min(ss)/max(ss)})) )) [,1] [,2] [,3] [,4] [,5] [,6] [1,] 1000.0000000 2000.0000000 3000.0000000 4000.0000000 5000.0000000 6000.0000000 [2,] 0.1684101 0.3033053 0.3869192 0.4466518 0.4883167 0.5158891 [,7] [,8] [,9] [,10] [,11] [,12] [1,] 7000.0000000 8000.0000000 9000.0000000 10000.0000000 11000.0000000 12000.0000000 [2,] 0.5387887 0.5611546 0.5857873 0.6023973 0.6199602 0.6338816 [,13] [,14] [,15] [,16] [,17] [1,] 13000.0000000 14000.0000000 15000.0000000 16000.000000 17000.0000000 [2,] 0.6371665 0.6508695 0.6613642 0.665185 0.6797338 [,18] [,19] [,20] [,21] [,22] [1,] 18000.0000000 19000.0000000 20000.0000000 21000.0000000 22000.0000000 [2,] 0.6845181 0.6965812 0.6978221 0.7095483 0.7138271 [,23] [,24] [,25] [,26] [,27] [1,] 23000.0000000 24000.0000000 25000.0000000 26000.000000 27000.0000000 [2,] 0.7187094 0.7274481 0.7267149 0.733775 0.7309658 [,28] [,29] [,30] [1,] 28000.0000000 29000.0000000 30000.0000000 [2,] 0.7380007 0.7450008 0.7422353

Exempel: I vårt lilla samlarkortnätverk var vi väl 10 stycken och köpte kanske 150
kort per person, dvs total 10*150 = 1500 inhandlade kort. I tabellen ser man att
förhållandet mellan antalet färst kort och antalet flest kort är 0.25 (1/4), vilket
troligen skulle kunna uppfattas som den E-effekt vi talade om i början.

Ännu större inköp som tar riktigt lång tid att simulera: för 100000 inköpta kort blir medelvärdet av min/max cirka 0.85 (0.8506242). För en miljon inköp är det 0.9512054.

Vi ser alltså att min/max-förhållandet tenderar att utjämnas då antalet inköp ökar, vilket är vad man kan förvänta sig enligt de stora talens lag.

Fullständiga serier i nätverket
Ser man på minimivärderna för dessa får man antalet (genomsnittligt uppkomna) fullständiga serier i nätverket. Detta förutsätter alltså att alla i nätverket byter rakt av och inte är snikna eller har orealistiska byt-krav.

> sapply(1:30, function(n) rbind(n*1000, mean(replicate(100,{ss<-collector.group.sim(100, n*1000);min(ss)})) )) [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [1,] 1000.00 2000.00 3000.00 4000.00 5000.0 6000.0 7000.00 8000.00 9000.00 [2,] 3.08 9.71 17.08 24.98 33.2 41.2 49.87 59.03 66.67



        [,10]    [,11] [,12]    [,13]    [,14]    [,15]    [,16]    [,17]

[1,] 10000.00 11000.00 12000 13000.00 14000.00 15000.00 16000.00 17000.00

[2,]    76.01    84.79    94   102.37   111.52   120.41   128.45   138.62

        [,18]    [,19]    [,20]    [,21]   [,22]    [,23]    [,24]    [,25]

[1,] 18000.00 19000.00 20000.00 21000.00 22000.0 23000.00 24000.00 25000.00

[2,]   147.44   156.05   165.46   173.09   184.5   192.92   202.72   212.33

[,26] [,27] [,28] [,29] [,30] [1,] 26000.00 27000.00 28000.00 29000.00 30000.00 [2,] 220.74 229.74 238.11 248.47 257.23

Dvs om 3 (3.08 i tabellen ovan) personer tillsammans köper 1000 kort får de (i medelvärde) alla får en fullständig serie. På samma sätt om 9 (9.71) personer köper 2000 kort, 17 (17.08) köper 3000 kort så bör de i genomsnitt få fullständiga serier. Osv.

Fokuserar här in på 100 .. 3000 kortinköp:
> sapply(1:30, function(n) rbind(n*100, mean(replicate(100,{ss<-collector.group.sim(100, n*100);min(ss)})) )) [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [1,] 100 200 300 400.00 500.00 600.00 700.00 800.00 900.0 1000.00 1100.00 [2,] 0 0 0 0.13 0.57 1.05 1.48 1.97 2.5 3.17 3.69



       [,12]   [,13]   [,14]  [,15]  [,16]   [,17]  [,18] [,19]   [,20]   [,21]

[1,] 1200.00 1300.00 1400.00 1500.0 1600.0 1700.00 1800.0  1900 2000.00 2100.00

[2,]    4.44    4.99    5.63    6.1    7.3    7.97    8.3     9    9.91   10.41

       [,22]   [,23]   [,24]   [,25]   [,26]   [,27]   [,28]   [,29]   [,30]

[1,] 2200.00 2300.00 2400.00 2500.00 2600.00 2700.00 2800.00 2900.00 3000.00

[2,]   11.17   11.95   12.79   13.53   13.93   15.09   15.85   16.16   17.21

Här ser vi också en-person-problemet: för 600 inköp är det drygt en person som får en full serie, det "exakta" väntevärdet är alltså 519.

Några noteringar
Lite av det vi set ovan:

* Om det är få inköps (säg < 700) får man alltså den 1/10-effekt som både E och jag upplevde i våra respektive barndomar. Det behövs alltså inte någon konspirationsteori för att förklara detta. (Naturligtvis kan detta kombineras med att kort-företaget faktiskt gjorde denna 1/10-tillverkning, men dessa experiment stödjer alltså inte detta.)

* Nätverkets storlek, dvs antalet personer som ingår i det, är inte viktig i denna modell, utan det som är avgörande är hur många kort man tillsammans köper. Det kan vara en person som köper 10000 och de andra 100 stycken, och tillsammans har de alltså 10000+ 10*(N-1) kort vilka sedan kan bytas hur som helst.

Samt en sista notering:

För 519 inköp i ett utbytesnätverk blir det igenomsnitt så här många fullständiga serier:
> table(replicate(1000,{ss<-collector.group.sim(100, 519);min(ss)}))/1000

0 1 2 0.409 0.563 0.028

Vilket är intressant eftersom det i 2.8% av fallen faktiskt skapas 2 stycken fullständiga serier (trist nog befann man sig aldrig i dessa nätverk :). 519 inköpta kort är, som sagt, det magiska medeltalet för att en person skulle få en fullständig serie.

Jämförelse: min/max-analys av svenska Lottodragningar
Här görs en motsvarande min/max-analys av svenska Lottodragningar. I går inhämtades fördelningen av senaste halvårets lottodragningar. (Värdena hämtas från Svenska Spel: klicka på den gröna Lotto-bilden/knappen, och i den uppkomna popuppen klicka på den gröna "Statistik"-knappen.)

Fördelningen av de 35 Lotto-talen är att det minst dragna talet har endast dragits 22 gånger (faktiskt två tal), och det mest dragna talet har dragits 44 gånger. En sorterad fördelning av talen ser ut så här:
> lotto2<-read.table("svensk_lotto2.dat",head=T, sep=",") > sort(lotto2$dragna) [1] 22 22 23 24 27 28 28 29 30 30 30 31 31 31 31 32 32 32 32 33 34 35 35 35 35 [26] 36 37 37 38 38 39 40 40 43 44



> sum(lotto2$dragna)

[1] 1144

> min(lotto2$dragna);max(lotto2$dragna)

[1] 22

[1] 44

> min(lotto2$dragna)/max(lotto2$dragna) [1] 0.5

Det ger alltså ett min/max-förhållande på 22/44 = 0.5. Är det något lurt här? Låt oss ser hur en simulering beter sig. Det är alltså totalt 1144 dragna tal de senaste halvåret och det finns 35 tal, vilket - för att prata i samlarkortproblemets terminologi - motsvarar 35 "kort" och 1144 "inköpta kort"

En exempelkörning:
> sort(ss<-collector.group.sim(35, 1144)) [1] 22 22 22 24 24 26 27 27 28 29 29 29 31 31 32 32 32 33 33 33 33 34 34 35 36 [26] 37 39 39 39 40 41 42 43 43 43

> min(ss) [1] 22 > max(ss) [1] 43 > min(ss)/max(ss) [1] 0.5116279

Medelvärdet för 100 sådana min/max-simuleringar är cirka 0.47, vilket är tillräckligt nära Lotto:s min/max-värde på 0.5 för att vi ska känna oss trygga att Lotto-dragningar inte är påtagligt konstiga.

> mean(replicate(100, {ss<-collector.group.sim(35, 1144);min(ss)/max(ss)})) [1] 0.466839

Slutkomihåg
I en samling slumpmässiga dragningar som Lotto eller samlarkort är det alltid någon eller några bollar/kort/tal som slumpmässigt dras fler gånger än de övriga, liksom det är alltid någon eller några som dras färre gånger än de övriga.

Man kan förledas att tro att det är något signifikant (speciellt, konstigt, magiskt) med dessa olika extremvärdena, speciellt om det är ett mindre antal dragningar som studerats. Sådana felslut kallas för "The extreme value fallacy" (se nedan).

Se även
Extremvärdesanalys av webbesök, speciellt länken till Number Watch: The extreme value fallacy.

Födelsedagsparadoxen/födelsedagsproblemet: Se t.ex. Sammanträffanden - anteckningar vid läsning av Diaconis och Mosteller 'Methods for Studying Coincidences' (teknisk, men länkar i slutet till andra sidor) samt .

Födelsedagsproblemet kan kanske ses som ett komplement till samlarkortproblemet: I födeledagsproblemet studerar man hur många personer det minst måste vara för att två (eller fler) ska ha samma födelsedag. I samlarkortproblemet är det frågan om hur många personer det måste vara för att alla dagar på ett år ska ha åtminstone en person som har födelsedag på denna person: Enligt det klassiska samlarkortproblemet (för en person som själv samlar kort) är det 365*sum(1/(1:365)) ~ 2365.

En Java-simulering: Coupon Collector Problem, från den trevliga sajten Probability by Surprise.

Posted by hakank at 07:27 EM Posted to Matematik | Sammanträffanden | Statistik/data-analys | Comments (2)

augusti 06, 2005

Lite om konsekutiva talpar i Lotto (och dess varianter)

I Konstantinos Drakakis paper Distances between the winning numbers in Lottery (arXiv, PDF) undersöks sannolikheten för en förekomst av konsekutivt talpar i Lotto-spel. Man ger en formel för denna typ av sannolikhet, och kommer t.ex. fram till att det för ett 6/49-Lotto (dvs där 6 nummer dras av 49) är det ungefär 49.5% chans att en Lotto-rad innehåller åtminstone ett konsekutivt talpar. Vilket är förvånansvärt högt.

Här nedan är några kommentarer och utvikningar, bl.a. med hjälp av riktig data för den Lotto-form vi har i Sverige, 7/35-Lottoformen (dock ej svensk data) samt simulering av sådan Lotto-typ.

Not: Sorterad dragningslista
En not kan vara på sin plats så här i början. Det är frågan om (konsekutiva) talpar i den sorterade listan av de dragna Lotto-numren. Det gäller alltså inte talpar som uppstår i själva dragningen. Exempel: Om en dragning görs av följande "bollar" (1, 14, 34, 2, 19, 15, 20) så tittar vi här på den sorterade listan (1, 2, 14, 15, 19, 20, 34) som ha tre stycken konsekutiva talpar: (1,2), (14,15) samt (19,20).

Notera också att varken Drakakis eller jag här nedan har brytt oss om tilläggstal och andra finesser.

Matematisk formel
Problemet formuleras på följande sätt i papret:

What is the probability that, out of m > 0 numbers drawn uniformly randomly from the range 1,..., n > m, at least two are consecutive?

Dvs vad är sannolikheten att det finns minst ett konsekutivt talpar (dvs ett eller flera) om man drar m tal från intervallet 1..n.

Den slutna matematiska formeln för denna sannolikhet uträknades till
p(n,m) = 1 - p(n,m)= 1 - ( C(n-m+1, m) / C(n,m) )

där C motsvarar (binomialkoefficienten ("n över m") och som i sin tur är definierad som n!/(k!*(n-k)!). Papret visar med matematisk stringens hur man kommit fram till detta resultat, men det vågar jag inte ha någon åsikt om och kommer inte heller att beröra detta här nedan.

Konkret
Låt oss i stället bli konkreta. För den ovan nämnda Lotto-typen 6/49 studerade värdena (n=49, m=6) får man sannolikheten för att det ska finnas minst ett konsekutivt talpar om man drar 6 tal från intervallet 1..47.

I det symboliska matematisk paketet MuPad (fritt att ladda ner och fritt att använda) uttrycks det på följande sätt.
>> n:=49:m:=6:1- binomial(n-m+1,m) / binomial(n,m);float(%);



                                22483/45402

0.4951984494

dvs det är cirka 49.5% chans att det finns åtminstone ett talpar. Förvånande mycket, eller hur?

Simulering
Min vana trogen har jag även simulerat detta, och lika troget är det gjort i statistikpaketet R. I simuleringen gjordes 100000 "dragningar" - med sample(1:49, 6) - och sedan räknar man hur många tal som ligger intill varandra. (Som vanligt använder jag den kompakta one-liner-formen för denna typ av simuleringar. En liknande simulering har förklarats lite mer ingående här.)

> sum(sapply(1:100000, function(x) {ddd<-diff(xxx<-sort(sample(1:49,6)));sss<-sum(ddd==1);sum(sss>=1)})/100000) [1] 0.49559

dvs 49.5%, vilket stämmer hyfsat bra med formeln.

Den svenska Lottoforme(l)n: 7/35
OK, denna typ av Lotto kanske inte är så intressant för oss. Vad med den svenska formen 7/35, dvs där man drar 7 nummer av 35? Först tar vi den teoretiska beräkningen (i MuPad):

>> n:=35:m:=7:1- binomial(n-m+1,m) / binomial(n,m);float(%);



                                8903/11594

0.7678971882

Dvs det är nästan 76.8% chans att det finns åtminstone ett konsekutivt talpar i de svenska lottoraderna. En simulering av 100000 dragningar ger (nästan) samma resultat 76.6%:
> sum(sapply(1:100000, function(x) {ddd<-diff(xxx<-sort(sample(1:35,7)));sss<-sum(ddd==1);sum(sss>=1)})/100000) [1] 0.76618

Utvidningar av detta
Kommen så här långt funderade jag på några andra saker kring Lottodragningar:

- finns det data för den svenska Lotto-formen?
- hur många konsekutiva talpar finns det i genomsnitt?
- hur är fördelningen för andra typer av talpars-differenser?

Data för Lotto på formen 7/35 (dock ej svensk)
Det vore ju roligt att se hur 7/35 ser ut i praktiken. Finns det då någon lista över en massa dragningar på det detta "svenska" Lotto-format 7/35? Jepp, det finns det. Mattias Hästbacka har samlat resultaten av Lottoraderna på sin Lotto-sida. Den fullständiga listan finns på Lotto-rader från v40/1980 till v17/2005, dvs nästan 25 år.

Observera att detta inte är Svenska Spels Lottodragningar, utan motsvarande "Finska Spel". Men det spelar egentligen ingen roll för det jag vill komma åt.

Om någon har motsvarande dataanhopning för Svenska Spels Lotto så vore jag mycket intresserad av få tillgång till den. (Maila mig gärna i så fall på hakank@bonetmail.com".)

Jag har bearbetat datan och i filen lotto_all_just_numbers.txt finns "rena" rader utan någon annan information, förutom en tillbörlig källhänvisning till Mattias Hästbackas sajt.

Analys av data
OK, tillbaka till undersökningen. Det finns i Mattias Hästbackas datamängd 957 rader med åtminstone ett konsekutivt talpar av 1283 dragningar, vilket är 957/1283 = 0.746, dvs 74.6%.

Det teoretiska värdet är som vi såg ovan 76.8%, vilket skulle motsvara 1293 * 0.7678971882, dvs ungefär 985 rader. Skillnanden mellan det praktiska värdet (957) och det teoretiska (985) är 985 - 957 = 28, vilket är ett fel på 28/1283, dvs ungefär 2.2% , vilket väl får anses vara rätt hyfsat.

Hur många par blir det i genomsnitt?
Den andra frågan är något annorlunda än Drakakis ursprungliga. Hans utgångspunkt var hur stor sannolikheten är för att det ska bli ett eller flera konsekutiva talpar (dvs åtminstone ett) i en dragning. En annan fråga man kan ställa sig är: hur många konsekutiva talpar det i genomsnitt blir per dragning.

Här är den första dragningen som Mattias Hästbacka sparade (vecka 40 1980):
7, 8, 9, 16, 17, 18, 30

Det finns alltså fyra konsekutiva talpar (7,8), (8, 9), (16, 17) samt (17, 18). Så, frågan är då hur många sådana talpar det blir i genomsnitti per dragning.

En körning på Hästbackas data ger totalt 1429 konsekutiva talpar och på 1283 dragningar blir det 1429/1283 = 1.113796. Dvs för varje dragning blev det igenomsnitt 1.1 konsekutivt talbar.

Äpplen och päron
Men, kan man nu fråga sig, hur rimmar detta 1.2 talpar per dragning med att det finns 76.6% chans för att det blir åtminstone ett konsekutivt talpar per dragning? Borde sannolikheten inte bli 100% istället? Nej, det är två olika saker (olika äpplen och päron).

Konstantinos formel tar inte hänsyn hur många konsekutiva talpar det är utan endast om det är ett eller fler. Dvs det spelade ingen roll för honom om det är 1 sådant par eller 4. Denna senare uträkning tar däremot hänsyn till detta antal. Den andra skillnaden (päronet alltså) är att vi inte räknar med sannolikheter längre utan med det genomsnittliga antalet per dragning.

Fördelning av konsekutiva talpar
En fördelning av antalet konsekutiva talpar i Hästbacka-datan är
0: 325 (0.25) 1: 565 (0.44) 2: 318 (0.25) 3: 69 (0.05) 4: 4 (0.00) 5: 1 (0.00)

Dvs det finns 325 Lotto-rader utan något konsekutivt talpar (cirka 25% av alla dragningar), 565 rader med exakt 1 konsekutivt talpar (~44%) och en rad med 5 par, nämligen den 581:e dragningen (vecka 47 år 1991) som gav följande tal 15 32 33 34 35 36 37.

Låt oss nu för skoj skull simulera detta. Koden som används ser ut så här:

num.talpar = function(num.sims=1283, n=35, m=7) { # initiera tabellen som kommer att innehåla antalet # konsekutiva talpar. Det kan vara max m-1 (här 7-1=6) sådana talpar, # men vi måste även ha med 0, vilket innebär att det blir m stycken konsek.talpar = rep(0, m) for (i in 1:num.sims) { ss = sort(sample(1:n, m)) dd = diff(ss) # differensen mellan talen tt = sum(dd==1) # hur många konsekutiva talpar finns det # öka med 1. Not: "tt+1" är för att även få med dragninar som # inte har något konsekutivt talpar konsek.talpar[tt+1] = konsek.talpar[tt+1] + 1 } return(konsek.talpar) }

Här är en sådan körning som har värden som inte är helt olika de faktiska (dvs Hästbacka-datan). Först antal och sedan fördelningen.
> (ss = num.talpar(1283)) [1] 314 519 353 88 8 0

> round(ss/1283,2) [1] 0.24 0.40 0.28 0.07 0.01 0.00

För just denna simulering blev det alltså cirka 314 dragningar (cirka 24%, av 1283 rader) utan något konsekutivt talpar alls, 519 dragningar med ett sådant talbar, och hela 8 dragningar med 5 k.t.p (dvs där det finns 6 konsekutiva tal).

Om vi nu brassar på med 100 sådana här simuleringar om vardera 1283 dragningar, dvs 1283*100 = 128300 dragningar (och delar de värden vi får med 100) så får vi något mer exakta värden att jämföra med:

> (ss=num.talpar(1283*100)/100) [1] 299.31 540.63 339.48 92.53 10.60 0.45 0.00

> round(ss/1283,2) [1] 0.23 0.42 0.26 0.07 0.01 0.00 0.00

Vi ser alltså att det finns cirka 23% chans att det blir 0 stycken k.t.p., cirka 42% att det blir exakt ett k.t.p. osv.

Den intressanta frågan blir då vad är det (simulerade) genomsnittet konkekutiva par per dragning. I faktisk data blev det cirka 1.1 konsekutivt par per dragning.

sum.talpar = function(num.sims=1283, n=35, m=7) { antal.talpar = 0 for (i in 1:num.sims) { ss = sort(sample(1:n, m)) dd = diff(ss) # differensen mellan talen tt = sum(dd==1) # hur många konsekutiva talpar finns det antal.talpar = antal.talpar + tt } return(antal.talpar) } sum.talpar()

Först en enkel körning för att få en känsla för informationen.
> ss = sum.talpar(1283) [1] 1580

> ss/1283 [1] 1.231489

Och sedan en lite tyngre körning

> sum.talpar(1283*100)/(1283*100) [1] 1.200514

Dvs det blir cirka 1.2 konsekutiva talpar per dragning. Det är alltså mer än ett par i genomsnitt.

Hur fördelas differenserna av talen?
Ovan har vi endast studerat konsekutiva talpar, dvs tal där differensen är 1. Nen det är även intressant att studera övriga differenser som uppkommit i dragningarna. Det vi ovan kallat för "konsekutiva talpar" är talpar med differensen 1, vilket vi i det följande kallar "1-differens".

I Hästbacka-datan fördelas sig differenserna enligt följande (på formatet differens: antal förekomster). Notera att en dragning nästan alltid har flera olika differenser. T.ex. ovanstående 5-talpars-dragning (15 32 33 34 35 36 37) har fem stycken 1-differenser samt en 17-differens (dvs 32 - 15 = 17).
1: 1429 2: 1184 3: 939 4: 790 5: 683 6: 585 7: 447 8: 345 9: 283 10: 231 11: 207 12: 126 13: 108 14: 99 15: 66 16: 48 17: 44 18: 17 19: 22 20: 11 21: 16 22: 2 23: 5 24: 2 25: 3

Man kan här se att det är flest 1-differenser, sedan blir det ständigt färre och färre.

OK, för fullständighets skull så simulerar vi även detta. Koden är rätt lik den föregående.
diff.table.sim = function(num.sims=1283, n=35, m=7) { diff.table = rep(0, n) for (i in 1:num.sims) { ss = sort(sample(1:n, m)) dd = diff(ss) # differensen mellan talen # lägg in respektive differens i tabellen for (d in dd) { diff.table[d] = diff.table[d] + 1 } } return(diff.table) }

Det vi söker är alltså det genomsnittliga antalet differens per dragning, så vi simulerar 100000 dragningar.
> cbind(diff.table.sim(100000)/100000) [,1] [1,] 1.19969 [2,] 0.98629 [3,] 0.80547 [4,] 0.65678 [5,] 0.52912 [6,] 0.42424 [7,] 0.33684 [8,] 0.26439 [9,] 0.20978 [10,] 0.15981 [11,] 0.12235 [12,] 0.09001 [13,] 0.06614 [14,] 0.04662 [15,] 0.03512 [16,] 0.02358 [17,] 0.01583 [18,] 0.01097 [19,] 0.00734 [20,] 0.00406 [21,] 0.00281 [22,] 0.00149 [23,] 0.00073 [24,] 0.00033 [25,] 0.00013 [26,] 0.00006 [27,] 0.00002 [28,] 0.00000 [29,] 0.00000 [30,] 0.00000 [31,] 0.00000 [32,] 0.00000 [33,] 0.00000 [34,] 0.00000 [35,] 0.00000

Varvid vi kan se att det är i genomsnitt (nästan) 1.2 stycken 1-differenser (dvs det vi ovan kallade konsekutiva talpar) per dragning, i genomsnitt 0.99 2-differens per rad osv.

Genomsnittlig antal differenser
Det genomsnittliga antalet differenser är cirka 4.53 enligt följande simulering. Här delas summan av differensvärdena med det totala antalet uppkomna differenser. (Man kan trixa med föregående värden, men jag vill ha lite kontroll över saker och ting.)
diffs = function(num.sims=1283, n=35, m=7) { diffs = 0 num.diffs = 0 for (i in 1:num.sims) { ss = sort(sample(1:n, m)) dd = diff(ss) for (d in dd) { diffs = diffs + d num.diffs = num.diffs + 1 } } m = diffs/num.diffs # genomsnittligt antal return(c(diffs, num.diffs, m)) }

Och så körningen, där vi endast är intresserade av det sista värdet:

> round(diffs(1283*100),3) [1] 3465204.000 769800.000 4.501

Vilket alltså ger det genomsnittliga differensvärdet 4.501. Äh, vi säger 4.5.

Sammanfattning
Här ovan har vi kommit fram till bl.a. följande:

* Sannolikheten att det ska förekomma åtminstone ett konsekutivt talpar i en 7/35-dragning är cirka 77%
* Det blir i genomsnitt 1.2 konsekutiva talpar per dragning.
* Det genomsnittliga differensvärdet för en dragning är 4.5.

Som antytts så anser i alla fall jag att det är förvånande.

Upprinnelse (inspiration)
Upprinnelsen till allt detta var notisen Consecutive pairs in winning lottery numbers hos MathForge.

Se även
Se även min genomgång av födelseparadoxen,tillika med R-kod, i Sammanträffanden - anteckningar vid läsning av Diaconis och Mosteller 'Methods for Studying Coincidences'.

Jämför även med den simulering som gjordes för ett liknande problem (6/49-Lotto) i Lite om resampling, simulering, sannolikhetsproblem etc.. Sök efter "Litet problem i (engelsk) Lotto". (Det koms dock inte fram till exakt samma sak.)

Samt överhuvudtaget Simulation, probability theory etc, såsom Trisssimuleringen (Java applet).

Posted by hakank at 03:37 EM Posted to Matematik | Statistik/data-analys | Comments (7)

juli 23, 2005

Annals of Improbable Research: Namntal (Name Number) för en profession

Även om jag just nu inte kan ge några bra exempel så är det inte helt sällan jag ser t.ex. i TV-intervjuer, eftertexter till filmer eller ute i allmänna livet, efternamn som får mig att fundera kring kopplingen mellan efternamn och yrke/profession/intressen. T.ex. en elektriker som heter Ström får väl en att börja undra, liksom en murare som heter Tegel.

I den excellenta tidskriften för otroliga forskningsresultat Annals of Improbable Research (AIR, Improbable Research är kanske mest kända för sitt Ig Nobel Prize) finns det att beläsa om s.k. Name Numbers (här nedan namnat "namntal") för en forskningsgren. Ett sådant namntal anger andelen personer inom en forskningsgren som har ett efternamn kopplat till denna forskningsgren.

Ursprungsartikeln för denna forskning skrevs av Kevin Krajik och var kring geologi i The "Name Number" for Geology, and for Other Professions, Annals of Improbable Research, vol. 11, no. 2, March/April 2005 (tyvärr ej funnen å lina och därför inte läst, men borde finnas här).

Den visade att namntalet för geologi är 1.35% (dvs att 1.35% av de som forskar inom geologi har namn som är kopplat till geologin). Emedan jag inte läst originalstudien får hänvisas till sekundär källa kring dess urval (från Astronomistudien, se referens nedan, och det är min alldeles egenhändigt befetade befetning):

Krajick calculated the Name Number for geology by dividing the number of geology-related surnames for those who presented papers at the 2003 meeting of the Geological Society of America by the total number of authorial surnames for that meeting. The geology Name Number presented in Krajick's study was 117 / 8639, or 0.0135432.

Namntalet för astronomi är 0.0027143 (0.27%), dvs betydligt lägre. Detta enligt Eric Schulman och Caroline V. Cox: The Name Number for Astronomy. Namntalsnamn här är Sun, Moon, Starr etc.

Det har även gjort en liknande studie kring statvetenskap, Richard Neimi (misspelled) The Name Number(s) for Political Science (PDF) som

... examined the 4,529 names (including those occurring more than once) appearing in the on-line index to the program for the 2005 national conference of the Midwest Political Science Association.

Författaren beskriver sedan flera namntalsvärden. Det första resultatet är 1.26% (vilket alltså är lägre än geologernas). Exempel på sådana namn: King, Prince, Good.

Sedan gör dock författaren en, enligt min mening, något farlig abrovink och antar att efternamn på amerikanska presidenter och andra statsmän är giltiga namntalskopplingar och kommer således upp i hela 2.08%. (Detta är - som metod betraktad - en farlig metod eftersom det räcker med att en eller ett fåtal presidenter har väldigt vanliga namn för att namntalet ska öka drastiskt och orättvist. Såvida man inte tillåter detta inom andra professioner också, t.ex. att fysiker som heter Newton, Pascal etc är namntalsnamn. Hmmm, dessa båda är faktiskt fysiska konstanter och borde vara namntalsnamn bara av det skälet. OK, vad med de som studerar filosofi och har samma namn som kända filosofer t.ex. Wittgenstein, Sarte och Marx? Nej, jag tycker inte att det borde vara tillåtet med endast efternamnkopplingar. Eller i alla fall att de i så fall indexeras på ett speciellt sätt såsom att räknas med en faktor f, t.ex. 1/3, eller sätts i en speciell kolumn till höger om de "riktiga" kopplingarna.)

Däremot gör författaren sedan några intressanta uppföljningsstudiefrågor vilket gör att man gärna och snart förlåter dennes abrovinklande:

Further study could look into the sub-specialties of the above-mentioned authors. Do Professors Washington, Adams, and so on study the presidency? Or at least American politics? Are Gandhi and Mao Indian and Chinese specialists, respectively? Do Wiseman, Goodman, Fair and Bliss study political philosophy? Do Power, Powers, and Guerra write about International Relations? From there, one could move on to consider whether persons take on any of the characteristics of the leaders (or kinds of leaders) whose names they bear. Are Law and Lawless opposites? Is Fair fair? And what should we expect from the Nixons? And, of course, one could see whether names and particular kinds of colleges and universities are linked. Are the Popes at Catholic institutions? Does Canon’s school have an ROTC program? These and many other fascinating questions await more detailed analyses.

Vänligen bemärk att AIR har en blogg: Improbable Research What's New -- News about research that first makes people LAUGH, and then makes them THINK. This is the official blog of the Ig Nobel Prizes and of the Annals of Improbable Research (AIR).

(Parentes 1: Det mesta ovan via just denna blogg, mer specifikt blogganteckningen Name number for political scientists.)

Slutfråga
För att avsluta med en slutfråga: När kommer den första svenska tillämpningen av denna nya forskningsgren? Och man skulle också kunna undra vilken typ av efternamn som ger namntalspoäng inom området.

En ännu sistare sak: Mitt efternamn, Kjellerstrand, har inga tydliga beståndsdelar som ger namntalspoäng i någon forskningsområde alls, vilket kanske skulle förklara en del eller annat. Några möjliga men åletade kopplingar: "Kjell" skulle kunna vara "källa" och kanske kunna kopplas till filosofi (sökande efter källan till vis(s)het) eller till vatten (jag dricker rätt mycket te), och "strand" kopplas till sand som vidare kopplas till att smula sönder saker (t.ex. namn) till dess minsta beståndsdelar, såsom just denna analys. Eller helt enkelt att jag är en lat person som tycker om att sitta vid en strand (läs: på balkongen eller i en park) och läsa en bok.

(Parentes 2: Den mer eller mindre vanligtvis observante läsaren kanske observerar kategoriseringen av denna anteckning till Sammanträffanden. Och det är ingen slump.)

Posted by hakank at 08:58 FM Posted to Sammanträffanden | Statistik/data-analys | Comments (0)

juni 28, 2005

InfoVis: The Landscape Metaphor

InfoVis: The Landscape Metaphor beskriver datavisualiseringar med hjälp av landskapmetaforen.

Whoever that has followed a driving course is able to catch the meaning of a traffic signal in a fraction of a second. Nevertheless if we show the same signal to someone that hasn't been taught how to interpret it or has never seen one, it will be very difficult for him or her to get the meaning in its complete magnitude, even if staring for a long time in front of it trying to decipher it ... In the end, it seems that decoding a map, be it spatial or thematic, results quite natural for most of the human beings, requiring just very little training. This doesn't appear to be the case with more recent developments of more abstract nature, like parallel coordinates or hyperbolic geometry, just to mention a couple of examples. ... The landscape metaphor, and spatialisations in general, constitute very attractive ways of representing big amounts of data in an intuitive way. The notions of distance and height are easily understood by many people, even when they encode other variables like similarity, density, or prevalence of a certain term.

En av de tekniker som beskrivs är Islands of Music (Matlab-kod finns att tillgå):

Islands of Music are a graphical user interface to music collections based on a metaphor of geographic maps. Islands represent musical genres. Mountains and hills on these islands represent sub-genres. The islands are situated in such a way that similar genres are close together and might even be connected by a land passage while perceptually very different genres are separated by deep sea. The pieces of music from the collection are placed on the map according to their genre. To support navigation on the map the mountains and hills are labeled with words which describe rhythmic and other properties of the genres they represent.

hakank.blogg

Anteckningar från en PopSci-junkie blandat med egna mundäna försök. Samt annat.