hakank.blogg

« Agent-baserad modellering - simuleringar av emergenta fenomen | Main | Vänsterhäntas dag »

augusti 12, 2003

Kort om R och simulering - ett praktikfall

Som utlovat kommer här en liten förklaring till R-koden som finns i Ännu ett verklighetesexempel. Det är två one-liners som simulerar sannolikheten för tärningskast med lite olika förutsättningar.

Målgrupp: Denna anteckning är skriven med Mats i bakhuvudet, men alla "du":n och "man" är generella och riktar sig till alla som är intresserade, inte bara Mats.

Jag kommer här endast att förklara en liten bråkdel vad R kan. Det är mycket kompetent vad gäller statistisk analys, såväl numeriskt som grafiskt.

Först lite kort om R och hur man får reda på mer om systemet. Huvudsidan till systemet är http://www.r-project.org/. Det finns att ladda ner via CRAN (The Comprehensive R Archive Network). Binärer finns för Linux, MacOS och Windows, och skulle man inte vara nöjd med det, finns även källkoden att ladda ner och kompilera själv.

CRAN innehåller även ett stort antal Contributed Packages, av mycket skiftande syfte.

Dokumentation om R, finns via vänstermenyns Documentation. Om du får mersmak för R, rekommenderar jag att man läsa igenom FAQ-dokumentet.

Det finns en liten söksida, där man dels kan söka efter funktioner, dels kan använda den hierarkiska listan för att hitta lämpliga funktioner. Notera att dokumenten som finns via denna även sida kopieras lokalt vid installationen (i alla fall för Linux-versionen). I systemet finns även lite hjälpfunktioner, men av lite olika skäl tenderar jag att använda mitt lilla externa Perl-program för att söka, och som jag gärna lånar ut.

Så, till förklaringen:

Det första exemplet

sum(sapply(1:1000000,function(x)all(sample(1:6,6,replace=T)==1:6)))/1000000

simulerar 1000000 kast med 6 tärningar för att räkna ut sannolikheten att dessa tärningar kommer upp ordnade enligt 1,2,3,4,5,6. Notera att det var inte riktigt det Mats ville ha. Men det är något enklare att förklara, så vi börjar med det.

Vi börjar att förklara inifrån:

sample(1:6, 6, replace=T)

1:6: innebär helt enkelt listan 1,2,3,4,5,6. I R kan det även skrivas c(1,2,3,4,5,6), där c() betyder "combine".

sample är funktionen för att göra ett slumpmässigt urval av elementen i listan. Parametrarna är först listan av möjliga värden (1:6), sedan hur många värden som ska returneras varje gång, (6). Det sista replace=T betyder att det är möjligt att returnera samma värde flera gånger (av de möjliga 6), T står för True. Några exempel på möjliga returvärden: c(1,2,3,4,5,6), c(1,1,1,1,1,1), c(2,1,2,1,2,1) etc. Om man skulle användareplace=F (F=False) innebär det att man inte tillåter upprepningar av elementen, men det är inte vår förutsättning (eftersom tärningarna helt oberoende kommer upp med värdena 1:6).

Nästa steg:

all(sample(1:6,6,replace=T)==1:6)

Här görs två saker. Först kontrolleras om det erhållna värdet från sample() är listan 1:6. Detta returneras en lista av True/False (T/F). T.ex. om listan som returneras från sample() är c(1,2,3,4,3,3) kommer kontrollen mot 1:6 (==1:6) att returnera c(T,T,T,F,F,F). Eftersom vi här endast är intresserade av att alla värdena ska vara True omsluter vi detta med all(), som just kontrollerar om alla värdena är True. Är så inte fallet, som i fallet med c(1,2,3,4,3,3), returneras False.

Nästa steg är att förklara vad sapply() gör.

sapply(1:1000000,function(x)all(sample(1:6,6,replace=T)==1:6)))

sapply() används här för att skapa de 1000000 (en miljon) anrop som krävs för simuleringen. Jag föreslår definitivt att man, när man börja skriva en simulering, startar med betydligt mindre tal, t.ex. 1000, för att se om saker och ting verkar stämma. Och man kan sedan öka detta tal beroende på önkvärd precision, tid och lust.

sapply är en av de funktionella operatorerna som finns i R. Parametrarna är sapply(lista, funktion), där lista är någon form av lista, här listan av talen 1 till 1000000 och är nog inga svårigheter att förstå. Nästa parameter är den funktion som ska köras för varje element i listan.

Som funktion har vi här en "anonym funktion", dvs det är inget anrop till en skriven funktion. Sådana anonyma funktioner måste vara på formen function(parametrar), där parametrar helt enkelt är en lista på de variabler som ska använda i den anropade funktionen. I vårt fall använder vi inga parametrar, dvs vi gör egentligen inte speciellt med talen från listan, utan de används endast för att skapa ett visst antal sample-körningar. Det måste dock finnas åtminstone en parameter till function, vilket jag envetet kallar x.

Så nu har vi alltså en lista med talen 1:1000000 och en anonym funktion all(sample(1:6,6,replace=T)==1:6)), som ska köras för varje element.

Som returvärde ger sapply en lista med värdena från alla anrop som gjorts, dvs den anonyma funktionen applicerad för alla element. Vi får här en lista av True/False beroende på om sample var av korrekt slag eller inte.

Slutligen (!) har vi det fullständiga uttrycket:

sum(sapply(1:1000000,function(x)all(sample(1:6,6,replace=T)==1:6)))/1000000

som gör två saker: Först görs en summering (med sum) på listan som sapply returerar. Notera att detta är en lista med True/False, men sum på en sådan lista räknar helt enkelt hur många True som finns i listan. Dvs sum() ger här antalet sample() som returnerar talen 1:6 i den avsedda ordningen (1,2,3,4,5,6).

Detta tal (som jag fick till 23 i min körning) måste vi sedan dela med antalet körningar (1000000) för att till slut få sannolikheten.

När jag väl insåg vad Mats menade (naturligtvis efter jag publicerade min lösning), att ordningen inte skulle spela någon roll, justerade jag helt enkelt funktionen till följande:

sum(sapply(1:100000, function(x) all(sort(sample(1:6,6, replace=T))==1:6)))/100000

dvs lade endast till sort(). Detta innebär att man man först sorterar listan innan man jämför med 1:6. Annars är det exakt samma som föregående.

Några slutkommentarer:
Dessa kärva one-liners kan skrivas mer strukturerat, men, som jag skrev i inlägget, är jag van att skriva på detta sätt för mindre simuleringar. Mer komplicerade saker måste man hantera något annorlunda, fast det är samma princip. T.ex. skriver man då funktionen som en namngiven funktion i stället för en anonym, men jag vill inte gå in på sådant här.

Har du några frågor om ovanstående eller om R i övrigt, får du gärna höra av dig. Systemet är som sagt mycket kompetent.

Posted by hakank at augusti 12, 2003 09:30 EM Posted to Statistik/data-analys

Comments

Posted by: Mats Andersson at augusti 12, 2003 09:47 EM